{'fraction_non_alphanumeric': 0.03460726585934896, 'fraction_numerical': 0.019997778024663928, 'mean_word_length': 4.404533173221255, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The KASCADE-Grande observatory and the composition of very high-energy cosmic rays

KASCADE-Grande is an air-shower observatory devoted to the detection of cosmic rays with energies in the range of 1016 to 1018 eV. This energy region is of particular interest for the cosmic ray astrophysics, since it is the place where some models predict the existence of a transition from galactic to extragalactic origin of cosmic rays and the presence of a break in the flux of its heavy component. The detection of these features requires detailed and simultaneous measurements of the energy and composition of cosmic rays with sufficient statistics. These kinds of studies are possible for the first time in KASCADE-Grande due to the accurate measurements of several air-shower observables, i.e., the number of charged particles, electrons and muons in the shower, using the different detector systems of the observatory. In this contribution, a detailed look into the composition of 1016 — 1018 eV cosmic rays with KASCADE-Grande is presented.


Introduction
A long standing problem in astroparticle physics has been the identification of the sources and the acceleration mechanism of galactic cosmic rays. Possible scenarios range from supernova remnants (SNR) [1], superbubbles [2] and magnetars [3] to cannonballs from supernovae [4] and blast waves from an early gamma-ray burst in our galaxy [5], being the SNR model the accepted paradigm. Modern γ-ray observations seems to support the SNR and the superbubble scenarios [6], however, the question whether these are the accelerators of the very high-energy cosmic rays (E > 10 2 TeV) is still open. Undoubtedly, the start of neutrino astronomy with ICECUBE [7], the arrival of new gamma-ray detectors, such as HAWC [8] and CTA [9], together with detailed measurements of galactic cosmic rays will pave the way to settle down the question as to the origin and acceleration of these energetic particles.
Nowadays, there are different models to explain how cosmic rays can reach high-energies at different sources. However, there exist some common characteristics shared among them: they involve astrophysical plasma motions and magnetic fields and predict the existence of individual knee-like structures or breaks in the corresponding spectra of cosmic nuclei, which are associated with the efficiency of the accelerator. Such structures lead to the formation of a knee at around (3 − 5) · 10 15 eV in the all-particle energy spectrum, which is a prominent feature discovered in the cosmic ray flux caused by a steepening of the spectrum [10]. According to the theory of cosmic rays, the positions of the knees are distributed inside the energy interval from E = 10 15 to 10 18 eV and, depending of the model, they scale with the electric charge, Z, of the ion (E Z ∝ Z × B × R, with B: the magnitude of the magnetic field and R: the size of the confinement region) [1,5,11], or the atomic mass, A, of the nucleus being accelerated [4]. Further differences among the models can also give rise to distinct predictions on the relative abundances of each element in the primary flux and the existence of fine structures on the overall cosmic ray energy spectrum between 10 15 and 10 18 eV.
The spectrum observed at Earth is not only shaped by the acceleration mechanisms working at the sources but also by the propagation process of cosmic rays inside the galaxy, which is also a matter that is not completely understood (see [12,13] and references therein). An important ingredient of the propagation models of cosmic rays is the interstellar magnetic field, which has a deep influence on the paths of galactic cosmic rays and their spectrum. In general, the particular effects of the intergalactic magnetic fields on the shape of the all-particle spectrum depend on the model under consideration. Models can be found in which the knee-like feature in the total energy spectrum of cosmic rays is explained as a natural consequence of the loss of efficiency of the galactic magnetic field at confining energetic charged particles in the Milky Way. In this context, the leakage of cosmic rays from the galaxy would imply also the presence of kinks in the spectra of the different kinds of cosmic ray nuclei. The kinks would be characterized by a rigidity dependent position.
Some clues to elucidate among the aforementioned scenarios are to be found in the study of the properties of cosmic rays, i.e., the analysis of their energies, arrival directions and composition. First measurements performed by the KASCADE experiment 1 in the energy range E = 10 15 − 10 17 eV revealed in fact the existence of individual knees in the spectra of the light and intermediate mass groups of cosmic rays (Z ≤ 14) [14], which was confirmed by the EAS-TOP [15], Tibet [16] and GAMMA [17] observations. The knees are formed because the fluxes of the corresponding primaries begin to decrease with the energy. Interestingly, KASCADE discovered that the observed depletions, particularly, in the energy spectra of the light mass groups (Z < 6) were responsible for the origin of the knee feature of the all-particle cosmic ray flux [14]. From the results of the experiment a pattern emerged, in which the positions of the individual knees are shifted to higher energies as the atomic/mass number of the primaries increases. However, these data were not enough to rule out a Zor Adependence of the individual knees due to the uncertainties in the measured locations of the breaks 2 . A convincing argument in favor of one or another type of dependence has to come from the measurement of the knee of the heavy mass group (for which Fe is the representative element) in cosmic rays. That structure was missing in the KASCADE data, presumably due to the lack of statistics at the energies of interest. Assuming E H ≈ 3 · 10 15 eV [14] for the position of the knee of the hydrogen mass group and supposing that the atomic or mass dependence holds for the location of the individual knees, the break for the iron nuclei is expected at energies in the interval E F e = (Z · E H , A · E H ) ≈ (7.8 · 10 16 eV, 1.7 · 10 17 eV), outside the sensitivity range of the KASCADE experiment. To explore this energy regime in search for the so called iron-knee, KASCADE was upgraded to KASCADE-Grande [18] using the components of the EAS-TOP observatory. During the period from 2003 until the date of its shut-down, in 2012, measurements of cosmic ray events were performed with the KASCADE-Grande experiment in the interval 10 16 − 10 18 eV. So far several dedicated and careful analyses have been applied to the data concerning the study of the chemical composition and the all-particle energy spectrum of very high-energy cosmic rays. The effort has result in very interesting results [19,20,21,22]. In this contribution, we will take a fast look into the KASCADE-Grande instrument and will sum up the details of some of the aforementioned analyses and their corresponding results.

Detection systems of KASCADE-Grande
Above ∼ 1 PeV, cosmic rays are studied by means of indirect techniques as a consequence of the low flux at the high-energy regime. The techniques are based on the detection of the extensive air showers (EAS) that cosmic rays induce in the Earth's atmosphere. In order to derive significant information about composition and primary energy the simultaneous observation of several EAS parameters [23,24] is imperative, e.g., the total number of electrons/muons in the shower, the density of particles in the shower front, the hadron content, the muon production height and the depth of the shower maximum. The measurements are done combining different kinds of particle detectors such as plastic scintillator detectors, tracking detectors, fluorescence and Cherenkov telescopes, water Cherenkov detectors and radio antennas among others [23,24]. The devices are arranged in arrays and are deployed over large surfaces on the ground in order to achieve big collecting areas. However, measurements are not only bounded to the ground. Fluorescence and Cherenkov telescopes, can be put in orbit around the Earth in order to study extensive air showers from space [25].
In case of KASCADE-Grande, we speak of a ground-based EAS array, which consisted of different detection systems. This instrument was designed to measure and separate the muon and electromagnetic components of the EAS [18]. The KASCADE-Grande detector (located at 49.1 • N, 8.4 • E, 110 m a.s.l.) was the successor of the KASCADE experiment [26] and incorporated the original electromagnetic detectors and muon devices of KASCADE to a bigger system of detectors, called Grande, composed of a 700× 700 m 2 array with 37× 10 m 2 scintillator stations regularly spaced by an average distance of 137 m [18]. A smaller array, named Piccolo was also added and was used to coordinate the KASCADE and Grande triggers. The main characteristics of the KASCADE-Grande detectors are displayed in table 1 and the layout of the experiment is presented in figure 1.
The Grande array was used to sample the density of charged particles of the shower front at ground level and to measure the particle arrival times of the EAS. The core position, the 2 Part of these uncertainties are associated with the hadronic interaction models.  number of charged particles (N ch ) and the arrival direction of the shower were extracted from the Grande data through an iterative fit and a careful modeling of the EAS front. To reconstruct the arrival direction of the shower from the arrival times a χ 2 fit was applied assuming a curved shower front as suggested by CORSIKA/QGSJET II simulations. On the other hand, to obtain the core position and the shower size from the density of charged particles of the EAS at ground, a modified NKG lateral distribution function was fitted to the data by means of a maximumlikelihood procedure [18]. An important component of this experiment was the KASCADE muon array, composed by 192 × 3.2 m 2 shielded scintillator detectors, which are sensitive to muons with threshold energy above 230 MeV for vertical incidence. With the Grande information and the measurements of the shielded array of the lateral distribution of muons in the shower front, the muon size (N µ ) was reconstructed event-by-event at KASCADE-Grande. The procedure involved the maximum log-likelihood technique along with a Lagutin-Raikin distribution function [18]. Once the muon content was known, the number of electrons, N e , in the EAS was estimated subtracting N µ from N ch .
Buried in a tunnel under several layers of soil, iron and concrete, the muon tracking detector (MTD) was installed close to the center of the KASCADE array. The MTD was composed of streamer tubes grouped in modules. It was used to reconstruct the individual tracks of penetrating particles of a section of the EAS and to measure the muon pseudorapidities [27]. By applying a triangulation procedure on the reconstructed tracks the mean muon production height of each event was also derived [28].
As it can be seen, different EAS parameters are at disposal for the research of high-energy cosmic rays at the KASCADE-Grande data base. However, due to the lack of space, only those analyses based on the muon and charged particle numbers will be considered for the present work.

Accuracy of EAS reconstruction
Systematic uncertainties for the core, N ch and arrival direction of the EAS are studied directly [18] by comparing the results of the Grande and KASCADE reconstructions, which work independently. Comparisons were performed for a subset of data with cores located inside a common area for both detectors and shower sizes in the interval log 10 N ch = 5.8−7.2. Accuracies of the EAS core positions and arrival directions are found to be of the order of 5 m and 0.7 • , respectively. Meanwhile, for the total number of charged particles a resolution ≤ 15% is achieved. Those values are in full agreement with expectations from Monte Carlo simulations.
To study the accuracy of the reconstruction of the muon size, Monte Carlo simulations have to be used. The reason is that shielded detectors were only available for KASCADE. In this case, the results show that, for N ch > 10 4.7 , the N µ systematic uncertainties are ≤ 25%, decreasing for high-energies. They exhibit a dependence with the core position and the arrival direction [19]. Since N µ is reconstructed with less accuracy than the charged particle number a muon correction function is applied to the data. This is possible because the behavior of muon uncertainties are well understood. The function is built from Monte Carlo simulations (using the high-energy hadronic interaction model QGSJET II-02 [29]) for a mixed composition assumption (H, He, C, Si, and Fe in equal abundances) and is parameterized as a function of the muon size, the EAS arrival direction and the distance of the EAS core to the KASCADE center.

Analyses and results
The available shower parameters to study the properties of cosmic rays with the KASCADE-Grande experiment are the muon and electron contents, the number of charged particles, the muon production height, the muon pseudorapidities and the muon/charged particle density at the shower front and at a given distance from the core. These parameters can be used in different ways to get insight into the energy and composition of cosmic rays (see, e.g., [19]- [22] and [30]- [33]) and, even more, to test the available high-energy hadronic interaction models [33]- [36]. In the next subsections, the main concern will be on the all-particle energy spectrum and the energy spectra of mass groups of cosmic ray particles in the sensitive range of KASCADE-Grande (E = 10 16 − 10 18 eV). In particular, the analysis techniques and the results based on the study of the EAS observables N ch and N µ will be summarized [19,20,21,22,30]. These analyses were applied over a subset of selected data that passed several quality cuts. Since unfortunately these cosmic-ray analyses must rely on MC simulations to interpret the EAS events several sets of simulated data had to be built. They included both the development of the EAS and its interaction with the detectors. These MC data were analyzed in the same way as the experimental ones. The production and development of the showers were simulated with CORSIKA [37] and the low and high-energy hadronic interaction models FLUKA [38] and QGSJET II [29]. Events were sampled from a differential spectrum described by a power law distribution (Φ ∝ E γ ) with spectral index γ = −2, which was later re-weighted to simulate a γ = −3 spectrum. The EAS were generated for the zenith angle interval θ = 0 • − 40 • using several types of primaries (H, He, C, Si and Fe). Finally, smaller sets of simulations were also created using EPOS 1.99 [39], SIBYLL 2.1 [40] and QGSJET II-04 [41] to check the sensitivity of the results to the high-energy hadronic interaction model. Left panel: Reconstructed all-particle energy spectrum, multiplied by E 3 , as obtained from the KASCADE-Grande N ch and N µ data and the different methods described in the present paper. The coloured band represents the preliminary systematic uncertainty for the spectrum reconstructed through the N ch − N µ method. Right panel: Residual plot after fitting the reconstructed all-particle energy spectrum derived from the N ch − N µ method with a power law function in the interval E = (2 − 8) · 10 16 eV. The systematic error band is also shown (dotted lines) [19].
4.1. The all-particle energy spectrum [19]. The reconstruction method in this case makes use of the correlation between N ch and N µ to reduce the uncertainties introduced in the result by the composition of the primaries. For the analysis, data with θ < 40 • is divided in different zenith angle intervals. Only EAS events with cores located in a central area of KASCADE-Grande were selected. The assignment of energy is done event-by-event according to a MC derived formula, which depends on the magnitudes of N ch and the parameter k = k(N ch , N µ ). The latter is defined as: These functions are calibrated for each zenith angle interval according to the predictions of the MC program CORSIKA and the FLUKA/QGJSET II hadronic interaction models for pure proton and iron nuclei (an E −3 energy spectrum is assumed). From this procedure, the energy spectrum for each zenith angle interval is derived. Before adding up all the spectra to obtain the final result, each flux is corrected for the effect due to the bin-to-bin migration of events. That is done building a response matrix from MC simulations and applying an unfolding procedure. The reconstructed all-particle energy spectrum is shown at the left panel of figure 2. It is based on data measured by KASCADE-Grande during the period from 2003 to 2009 with total exposure ≈ 2 · 10 13 m 2 s sr after applying selection cuts. On the other hand, a residual plot of the spectrum is presented on the right panel of figure 2. It reveals the presence of two significant structures in the all-particle energy spectrum, specifically, a concavity around 10 16 eV and a weak knee-like feature at 8.3 · 10 16 eV, where the spectrum changes from a power law behavior with spectral index γ = −2.92 ± 0.02 to another one with α = −3.39 ± 0.07. Both structures are significant with respect to the systematic and statistical uncertainties. A statistical analysis reveals, for example, that the spectrum above log 10 (E/eV) = 10 16 Figure 3: Left panel: Residual plot of the experimental all-particle energy spectrum reconstructed with the N ch − N µ method and using different hadronic interaction models for the energy assignment to the data [42,43]. Right panel: Experimental spectrum for the heavy mass group of cosmic rays derived from the N ch − N µ technique using four hadronic interaction models. The error bars represent the statistical uncertainties and the bands, the systematic ones. Spectral indexes obtained from a power-law fit to the energy spectra before and after the corresponding knees are indicated on the corners of the plot (left and right, respectively) [43]. the all-particle energy spectrum are not an artifact of the hadronic interaction model employed for the energy calibration [42,43]. In all cases, the features remain stable. Just a shift in the energy scale is observed among the spectra as seen on the left panel of fig. 3. On the other hand, recent results of the TUNKA [44], GAMMA [45] and ICETOP [46] experiments have also shown the presence of the concavity in their spectra around the energy scale reported by KASCADE-Grande. In addition, the TUNKA [44] and ICETOP [46] collaborations have reported the existence of the small break in the spectrum close to 10 17 eV. A change from a light mass group (Z = 1 − 2) to the intermediate one (Z = 6 − 12) in the all-particle energy spectrum could explain the concave behavior of the spectrum [47]. However, a transition between two different kinds of galactic sources could also have the same effect [48]. On the other hand, the origin of the observed break is also intriguing. That feature is related with the structure known as the iron-knee, as it will be seen later on in the section devoted to the composition studies.
4.1.2. N ch and N µ techniques [19]. The N ch and N µ observables can be employed individually for the task of reconstructing the energy spectrum. However, the procedure has the disadvantage that it is more sensitive to the primary composition, which is reflected in the form of a wide systematic error band in the energy spectrum. Such analyses are important for several reasons, for instance, in order to study the influence of systematic uncertainties, to cross-check the reconstruction methods and to investigate the validity of the used hadronic interaction model. The technique goes as follows: Applying the constant intensity cut method on EAS with θ < 40 • , the attenuation curves for the N ch and N µ parameters in the atmosphere are obtained. These curves are used to calculate the equivalent muon/charged particle number of the showers at a Right panel: All-particle energy spectrum and energy spectra (without unfolding) of the heavy (electron poor) and the light and medium (electron rich) mass groups derived by applying the k technique to the KASCADE-Grande data. Error bands are obtained by considering the changes in the spectra introduced by varying the frontier lines inside their respective interval of uncertainties [19].
Since the composition is a priori unknown, fits are performed for two extreme cases: pure protons and iron nuclei. To end the procedure, data derived from each observable under a given composition assumption are combined to construct the associated energy spectrum. The reconstructed spectra for the data from the period 2003−2009 are shown in the left panel of figure  2 (migration effects are taken into account). From the widths of the bands generated by the H/Fe solutions, it is clear that the N ch variable is more sensitive to the primary composition than the N µ case. For the N ch case, the region of the expected solution covers a bigger area than that obtained from the N µ analysis technique. However, it is seen some room for a common solution, where in fact the N ch − N µ spectrum is found. That speaks in favor of the internal consistency of the results when interpreted in the framework of the QGSJET II hadronic interaction model.

4.2.
The energy spectra of mass groups 4.2.1. k parameter technique [20]. One way to take a look into the composition of cosmic rays from the KASCADE-Grande data is by using the k parameter, already described in the last section. This quantity is sensitive to the composition of the primary particles (H, He, C, Si, Fe) and can be employed to separate the data into a light/medium and a heavy mass groups, at least. To achieve the separation, a plot with the mean value of k for different primaries as a function of the reconstructed energy is created from MC simulations for every zenith angle interval in which data is divided. Here, energy is estimated through formula (1). An example, of such a graph can be seen in the left panel of figure 4. It happens that, in general, the mean k value grows with the mass of the particle. That is not surprising, since the k parameter is defined in such a way that it returns the biggest values for iron nuclei and the lowest, for protons. With the above plots, data is divided into mass groups by defining a frontier line,k Si−C , which is obtained from a fit to the average between the mean k values of Si and C. If k >k Si−C , events are classified as part of the heavy component (also called electron poor sample). In other case, events are included in the light and medium group (electron rich sample). The energy spectra for   figure 4 (right panel) along with the all-particle energy spectrum obtained by adding the individual spectra. In this plot, the knee around 10 17 eV in the all-particle energy spectrum is seen again. But now, it is more clear what is its origin. The answer must be track down to the existence of a knee in the energy spectrum of the heavy component, whose relative abundance around 10 17 eV is bigger than that of the other mass groups. A fit performed to the spectrum of the heavy component using a broken power-law function with a smooth knee [49], indicates a change ∆γ = γ 2 − γ 1 = −0.48 of the spectral index and locates the knee at log 10 (E knee /eV) = 16.92 ± 0.04. The significance that this spectrum is not described by a single power-law is of 3.5σ. Applying the same fit to the all-particle energy spectrum results in ∆γ = −0.29 and log 10 (E knee /eV) = 16.92 ± 0.10. It is worth to mention that the presence of the knee-feature at the spectrum of the heavy component is independent of the high-energy hadronic interaction model (although relative abundances are strongly affected). That was checked out by reproducing the analysis within the framework of EPOS 1.99, SIBYLL 2.1 and QGSJET II-04 (see right panel of figure 4) [42,43]. In these analyses, the significance levels for the existence of the heavy knee in the data goes from 2.4, for QGSJET II-04, up to 7.4, in case of SIBYLL 2.1. Turning now the attention to the light/medium spectrum, the data seems to reveal that this component is present even up to 10 18 eV. That is an important issue. Some astrophysical models predict the existence of an extragalactic component dominated by protons in the cosmic ray flux, which would take over the galactic contribution in the energy regime of 10 17 − 10 19 eV [50]. Therefore, somewhere in this energy range, part of it accessible to KASCADE-Grande, clues of a possible galactic-extragalactic transition may be buried. One way to look for them is by examining the light/medium spectrum intently. For this purpose, a larger data sample was selected from a bigger data acquisition period 2003 − 2012 and for an extended fiducial area as described in [21]. The steps from the k-technique were also applied here. But additionally, to Error bars represents the RMS of the Y distributions. Right panel: Light and heavy energy spectra derived from the KASCADE-Grande data using the Y ratio as a mass estimator. The all-particle energy spectrum derived by adding the individual mass group fluxes is also shown [51].
separate the light from the medium contribution a new cut was introduced by defining another frontier linek C−He , which is obtained by averaging the mean k values of C and He at a given energy (see figure figure 5, left side). The resulting energy spectra for the light component is shown in figure 5 (right side). An ankle-like feature in the light component of cosmic rays can be clearly observed. A double power-law fit with formula (3) indicates that the hardening of the spectrum of the light component is characterized by a change of spectral index from γ = −3.25 ± 0.05 to γ = −2.79 ± 0.08 and occurs at log 10 (E ankle /eV) = 17.08 ± 0.08, slightly above the position of the heavy-knee. Analyses conducted in [21,43] point out that the changes in the spectrum of heavy and light primaries around 10 17 eV are not connected by a bias introduced by the classification procedure or reconstruction techniques. The appearance of the ankle-like feature is striking and may be indicative of an early transition from galactic to extragalactic cosmic rays. [20,51,52]. This technique makes use of the ratio between the equivalent muon and charged numbers at a zenith angle of reference (

Y -cut technique
which is also a sensitive quantity to the composition of cosmic ray particles. In equation (4) the equivalent number of particles of an EAS is estimated with the CIC method, as explained in the subsection 4.1.2. However, the energy of the event is calculated using the expression (1). As it was the case for the k parameter, plots for the Y ratio as a function of the reconstructed energy are produced for several elements in order to be used in the classification of events (c.f figure  6). The difference with the k parameter is that for the plots of the Y ratio, data from all zenith angles are considered to build up the graph. It is worth to mention that for this analysis a smaller data sample was selected (θ < 30 • ) in order to have events with lower systematic uncertainties (errors tend to increase with the value of the zenith angle). As seen in figure 6, for a fixed energy, the Y ratio grows with the mass of the primary nuclei. The effect can be understood as a result of the fact that physics of hadronic interactions favors the production of more muons and less electrons for heavy nuclei. From the same figure, it is also noted that the mean Y distributions are almost energy independent above the threshold (log 10 (E/GeV) ≥ 7.4). Taking advantage of the observed behavior of Y a cut on this parameter is applied at Y cut = 0.84 to divide data into two different sets: the electron rich (Y < 0.84) and electron poor (Y ≥ 0.84) groups. Using figure 6, it is clear that the first set covers the mean Y values of light mass elements (here H and He), while the second set contains the expected region for heavy mass nuclei (Si and Fe). In this way, the electron rich and electron poor samples become representative of the light and heavy mass groups (the fraction of misclassified events inside each group is less than 20% in the region of maximum efficiency).
The results for the energy spectra obtained from the 2003−2009 data sample are presented in the right panel of figure 6. Plots were not corrected for migration effects. From this figure, the knee-like structure of the heavy component of cosmic rays is immediately recognized. Applying again a fit with a broken power-law spectrum to the spectrum of the heavy mass group a change in the spectral index ∆γ = −0.47 at log 10 (E knee /GeV) = 7.77 ± 0.06 is obtained. The same analysis for the all-particle energy spectrum gives the following values: ∆γ = −0.29 at log 10 (E knee /GeV) = 7.92 ± 0.10. In this analysis, the light component of cosmic rays is also present. Relative abundances, of course, change with the position of Y cut , however, the presence of the knee-like future in the heavy mass group is independent of the cut, i.e. it is inherent to the cosmic ray data [52]. [30]. In order to avoid the influence of increasing systematic uncertainties with θ in the results, the present analysis is restricted to vertical showers (θ < 24 • ). From experimental data a bidimensional histogram of N µ /N ch vs log 10 (N ch ) is constructed. Similarly, two-dimensional histograms for simulated data organized into three mass groups (light, intermediate and heavy) are created. The light MC data set is composed solely of protons, the intermediate one, of He (50%) and C (50%), while the heavy one, of Si (50%) and Fe (50%). Then, inside a single bin log 10 N ch,j , each experimental N µ /N ch distribution, n exp j , is fitted with a function, n M C j , formed by the linear combination of those simulated distributions derived for the three mass groups (see, for example, left panel of figure 7). Here where i runs over all N µ /N ch -bins and ℓ, over each mass component. The fitting parameters are the α ℓ j coefficients, which represent the relative abundances of the different mass groups inside each log 10 N ch,j channel. Hence, the parameters must be constraint to the interval 0 ≤ α ℓ j ≤ 1 and satisfy ℓ α ℓ j = 1. The coefficients are extracted by minimizing the χ 2 functions: with σ i defined as the statistical error on the corresponding n M C ji entry. With the relative abundances, the number of events per mass group for each charged particle number interval is found. Finally, N ch is used to estimate the primary energy by using a power-law expression (E = ζN ǫ ch ) derived from MC simulations accordingly to each mass group. The energy spectra obtained from this method for data from the period 2003 − 2009 (exposure ≈ 8 · 10 11 m 2 s sr) are presented in the right panel of figure 8. Again, in this plot, the outstanding knee-like feature in the spectrum of the heavy mass group is visible. The change in the spectral index and the position of the knee-like feature are ∆γ = −0.62 and log 10 (E knee /eV) = 16.79 ± 0.04, respectively. This analysis also confirms the non-vanishing  [22]. The basis for this analysis is the two-dimensional histogram log 10 N ch vs log 10 N µ . The purpose is to unfold from that distribution the cosmic ray fluxes, dJ/d log 10 E, for five different mass groups (represented by H, He, C, Si and Fe). The problem is formulated in terms of the following coupled equations: n exp ji = ℓ dJ ℓ d log 10 E P ℓ (log 10 N ch,j , log 10 N µ,i | log 10 E) d log 10 E cos θdA dΩ dt, where P ℓ represents the conditional probability that an EAS with a given combination of charged particle and muon numbers, (log 10 N ch,j , log 10 N µ,i ), will be observed given that the primary cosmic ray belongs to the mass group ℓ and has an energy E, n exp ji is the number of entries in the bin (log 10 N ch,j , log 10 N µ,i ) of the experimental distribution, t is the observation time, Ω is the solid angle of observation, θ is the arrival zenith angle of the EAS and A is the fiducial area. P ℓ must be known apriori for the unfolding procedure. These probabilities comprise all relevant information regarding the properties of the observables of the shower, e.g. EAS fluctuations, and the efficiency as well as the reconstruction accuracy of the instrument/experiment. They are calculated using MC simulations (based on FLUKA 2002.4 and QGSJET II-02). To unfold the bidimensional distributions, Gold's unfolding algorithm [53] is applied. Solutions are checked out using also the Bayes' algorithm [54] and a method based on the principle of reduced crossentropy [55]. The corresponding bidimensional distribution used in this analysis (with θ < A knee-like feature in the spectrum of the Fe component is visible in figure 8. Remarkably, the kink of the heavy component obtained from this analysis is more pronounced 3 . A double power-law fit applied to the iron spectrum -see equation (3) -reveals a change of spectral index equivalent to ∆γ = −1.08 at the position log 10 (E knee /eV) = 16.9 ± 0.1. Moreover, the analysis shows that the iron primaries dominate in the high-energy regime (under the framework of QGSJET II). On the other hand, the all-particle energy spectrum does not exhibit the small knee around 10 17 eV. That could be the result of both the lower statistics of the subsample used for this analysis and the uncertainties that affect the unfolding of the five different mass groups, which are propagated to the all-particle spectrum.
With the position of the iron-knee at around 63 − 100 PeV and that for the proton-knee between 3 PeV and 5 PeV (according to the newest analysis of the KASCADE data, see first reference of [22] and figure 9) a rigidity-dependent picture for the position of the knees of the elemental primary cosmic rays seems to be preferred.

Conclusions
In this contribution, the composition of cosmic rays in the energy interval 10 16 −10 18 eV has been explored with the KASCADE-Grande experiment. Composition results from several analyses based on the study of the N ch and N µ observables of EAS were presented. The main conclusions of this research are that the spectrum of the heavy component of cosmic rays exhibits a kneelike feature around 10 17 eV, which is responsible for a subtle break in the all-particle energy spectrum around the same energy value. The spectrum of the light component of cosmic rays presents a hardening around 10 17 eV. Two mass groups, light and heavy (represented by H and Fe+Si nuclei, respectively), are not enough to describe the N µ /N ch distributions, 3 Here, it is important to mention one important difference with the composition analyses already described and it is that the separation of the mass groups by using the unfolding technique is based on the properties of each of the five elemental groups. In case of the other techniques, mass groups with a different degree of mixing of the elemental groups are employed. This difference could explain why the iron-knee is more pronounced in the unfolding analysis. QGSJet-II-02 / FLUKA Figure 9: The all-particle energy spectrum [19,21] and spectra of individual mass groups obtained with KASCADE [22] and KASCADE-Grande data [21]. Results are compared with the energy spectra measured with other experiments (figure taken from [56]).
an intermediate mass group is needed (at least for a QGSJET II-03 based analysis). From the different components, the heavy one is dominant inside the framework of the QGSJET II hadronic interaction model. Finally, relative abundances at very high-energies depend on the high-energy hadronic interaction model. Nevertheless, the main features observed in the all-particle and the heavy and light mass group spectra remain stable.