*Article* **Measurement of Energy Spectrum and Elemental Composition of PeV Cosmic Rays: Open Problems and Prospects**

#### **Giuseppe Di Sciascio**

INFN—Roma Tor Vergata, Department of Physics, University of Roma Tor Vergata, Viale della Ricerca Scientifica 1, I-00133 Roma, Italy; disciascio@roma2.infn.it

**Abstract:** Cosmic rays represent one of the most important energy transformation processes of the universe. They bring information about the surrounding universe, our galaxy, and very probably also the extragalactic space, at least at the highest observed energies. More than one century after their discovery, we have no definitive models yet about the origin, acceleration and propagation processes of the radiation. The main reason is that there are still significant discrepancies among the results obtained by different experiments located at ground level, probably due to unknown systematic uncertainties affecting the measurements. In this document, we will focus on the detection of galactic cosmic rays from ground with air shower arrays up to 1018 eV. The aim of this paper is to discuss the conflicting results in the 10<sup>15</sup> eV energy range and the perspectives to clarify the origin of the so-called *'knee'* in the all-particle energy spectrum, crucial to give a solid basis for models up to the end of the cosmic ray spectrum. We will provide elements useful to understand the basic techniques used in reconstructing primary particle characteristics (energy, mass, and arrival direction) from the ground, and to show why indirect measurements are difficult and results are still conflicting.

**Keywords:** cosmic ray physics; multi-messenger astrophysics; extensive air showers

**Citation:** Di Sciascio, G.

Measurement of Energy Spectrum and Elemental Composition of PeV Cosmic Rays: Open Problems and Prospects. *Appl. Sci.* **2022**, *12*, 705. https://doi.org/10.3390/ app12020705

Academic Editors: Roberta Sparvoli and Matteo Martucci

Received: 8 November 2021 Accepted: 23 December 2021 Published: 11 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

Cosmic rays (CRs) are the most outstanding example of accelerated particles and represent about 1% of the total mass of the universe [1]. The riddle of the origin of this radiation has been unsolved for more than a century. The study of CRs is based on two complementary approaches [2]:


In Figure 1, the primary CR all-particle energy spectrum (namely, the number of nuclei as a function of total energy) is shown. The spectrum exceeds 10<sup>20</sup> eV, showing a few basic characteristics [2]:


**Figure 1.** All–particle energy spectrum of primary cosmic rays, the flux is multiplied by E3. Results from direct and indirect experiments updated to the year 2021 are shown.

Despite the differences in flux, emphasized by multiplying the differential spectrum by E3, all the measurements of the all-particle energy spectrum are in fair agreement when taking into account the statistical, systematic and energy scale uncertainties. Nevertheless, uncertainties affecting flux measurements could be underestimated for a number of reasons discussed in this document. In a conservative approach, the spread of different results provides a more realistic estimate of the uncertainty.

All the observed features are believed to carry fundamental information that sheds light on the key questions of the origin, acceleration and propagation of CRs. However, from the all-particle results alone, it is not possible to understand the origin of different features. All models concerning sources, acceleration and propagation of the primary flux, differ considerably for what concerns expected elemental composition as a function of the energy. A measurement of the chemical composition is therefore crucial to disentangle between different hypotheses.

The main structure is the *"knee"* observed for the first time by R.W. Williams in 1948 in the experiment which first located individual shower cores from symmetry of the fired detectors [4,5]. The knee as a feature connected to the end of the Galactic CR flux was first suggested in 1959 by Kulikov and Khristiansen [6]. They speculated that particles above 1016 eV may have a *"metagalactic origin"*. Consequently, the observed spectrum is a superposition of the spectra of particles of galactic and metagalactic origin. In 1962, Miura and Hasegawa [7] reported the first observation of two spectral kinks (in both N*<sup>e</sup>* and N*<sup>μ</sup>* spectra) correlating them to a steepening of the primary energy spectrum.

All experiments observed the knee at about 4 × 1015 eV but a general consensus about the chemical component responsible for such a feature does not exist yet because experimental results are still conflicting, as will be discussed in Section 5. Determining elemental composition in the knee energy region is crucial to understand where Galactic CR spectrum ends and to give a solid basis to CR models up to the highest observed energies. The maximum energy at which the various nuclei are accelerated should be subject to a rigidity cutoff, as proposed originally by Peters [8]. Protons will cutoff first, followed by other nuclei according to the relation

$$E\_{\text{max}}(Z) = Z \times E\_{\text{max}}(Z=1) \tag{1}$$

If the dominant primary mass of the knee is light (protons and helium), then, according to this scheme, the Galactic CR spectrum is expected to end around 1017 eV with iron. The sum of the fluxes of all elements, with their individual knees at energies proportional to the nuclear charge, makes up the CR all-particle spectrum shown in Figure 1. With increasing energies, not only does the spectrum become steeper due to such cutoffs, but also heavier. In this scenario, the knee would represent the end of the spectrum of CR accelerated by SNRs in the galaxy.

Indeed, it is widely believed that the bulk of CRs up to about 1017 eV are galactic, produced and accelerated by the shock waves of SuperNova Remnants (SNR) expanding shells [9], and that the transition to extra-galactic CRs occurs somewhere between 10<sup>17</sup> and 1019 eV. The experimental results, however, do not demonstrate the capability of SNRs to produce the power needed to sustain the population of galactic CRs and to accelerate particles up to the knee, and beyond. Indeed, to accelerate protons up to the PeV energy domain, a significant amplification of the magnetic field at the shock is required, but this process is problematic [10].

Unlike neutrinos that are produced only in hadronic interactions of CRs, the question whether the observed *γ*-rays are produced by the decay of *π*<sup>0</sup> from CR interactions (*'hadronic'* mechanism), or by a population of relativistic electrons via Inverse Compton scattering or bremsstrahlung (*'leptonic'* mechanism), still needs a conclusive answer. In a hadronic interaction, the secondary photons have, on average, an energy factor of 10, lower than the primary proton. Therefore, the quest for CR sources to be able to accelerate particles up to the PeV range in a multi-messenger approach requires the observation of the *γ*-ray sky above 100 TeV. However, the first results reported by the LHAASO experiment [11,12], that is, the observation of a number of gamma sources emitting photons beyond 500 TeV, show that SNRs are likely not the main sources of PeV CRs in our galaxy. In fact, none of the 12 observed ultra-high energy gamma sources can be clearly described with hadronic mechanisms operating in SNRs. We note that the highest photon emission at 1.4 PeV comes from a system of massive stars in the Cygnus Region, the so-called *'Cygnus Cocoon'*, a possible factory of fresh CRs, as suggested by other experiments [13,14].

In this note, we will focus on galactic CRs in the PeV energy range detected from ground with air shower arrays. This is not a place for a complete review of CR physics and models (for which we recommend, for instance, [2,15–19] and the references therein), but only to provide elements useful to understand the main techniques used in reconstructing primary particle characteristics from the ground with particle arrays, and to show why indirect measurements are difficult and the results are still conflicting.

In the next section, we will introduce the detection techniques. In Section 3 we will describe the main characteristics of Extensive Air Showers to understand how different observables measured by arrays are related to the properties of the primary CRs. In Section 4, we will discuss the general scheme of the air shower array analysis. In Section 5 the experimental results in the 1014–1018 eV energy range are summarized. The prospects for new measurements in the knee region are introduced in Section 6.

#### **2. Detection Techniques**

We can divide the experimental methods adopted to measure fluxes and elemental composition of CRs into two categories: *'direct'* and *'indirect'* measurements. Generally speaking, for all particle types:


The direct measurements, in principle, detect and directly identify the primary particles with detectors outside the atmosphere (on board of stratospheric balloons or satellites), since the atmosphere behaves as a shield (see below). Since the CR flux rapidly decreases with increasing energy and the size of detectors is constrained by the weight that can be carried in flight, their *'aperture'* (i.e., the acceptance measured in m2·sr) is small and determines a maximum energy (of the order of a few hundred TeV/nucleon), at which a statistically significant detection is possible. In fact, the number of detected events is given by the CR flux times the detection area times the total observation time. Therefore, the detection area limits the smallest measurable flux. In addition, the limited volume of the

detectors makes the containment of showers induced by high-energy nuclei difficult, thus limiting the energy resolution of the instruments in direct measurements.

At higher energies, the flux is so low (about 1 particle/m2/year around 1015 eV) that the only chance is to have earth-based detectors of large area, operating for long times. In that case, the atmosphere is considered as a target, and we study the primary properties in an *'indirect'* way, through the measurement of secondary particles produced in the interaction of the primary particle with the nuclei of the atmosphere, the so-called *'Extensive Air Shower'* (EAS).

Approaching the hundred TeV energy region, even in space-borne experiments, the energy assignment is indirect since it is generally based on the energy deposition of particles produced in the interaction of primaries in the detector itself. The reconstruction of the total energy is then obtained by comparison with some model prediction, and therefore, at least in that region, the boundary line between 'direct' and 'indirect' experiments is more uncertain. In fact, important results obtained by 'direct' methods are conflicting due to some still unknown systematic uncertainties probably related to the interaction model used to assign the energy. A neutrino energy-dependent component must be estimated via Monte Carlo simulations, an evaluation which adds some additional model dependency for 'indirect' measurements.

At the ground, the study of CRs is based on the reconstruction and interpretation of EAS observables in the different components, electromagnetic (e.m.), muonic and hadronic, Cherenkov photons, nitrogen fluorescence, radio emission. Therefore, different detectors must be used to detect different observables.

Two different approaches are exploited:


Shower arrays are made by a large number of detectors (scintillators, Resistive Plate Chambers (RPCs) or water Cherenkov tanks, for example) distributed in a regular grid over very large areas, of an order of 104–105 m<sup>2</sup> (see Figure 2). The shower *"size"*, the total number of charged particles, and the shower arrival direction are the two key parameters reconstructed by all arrays. The majority of EAS arrays do not distinguish between the charged particles. From the measurement of the particle densities on the fired detectors of the array it is possible to determine the shower core position, that is, the point where the shower axis intersects the detection plane, and, via a Lateral Density Function (LDF), reconstruct the size of the shower. The LDF is of phenomenological nature, determined via Monte Carlo simulations for the particular experimental set-up [20]. The direction of the incoming primary particle is reconstructed with a *'time of flight'* method making use of the relative times at which the individual detection units are fired by the shower front [20].

**Figure 2.** Example of a typical air shower array (Tibet AS*γ* experiment located at the YangBaJing Cosmic Ray Observatory in Tibet (P.R. China) 4300 m asl).

On general grounds, the instrumented area A determines the rate of high energy events recorded, that is, the maximum energy via limited statistics. The grid distance *d* determines the low energy threshold (small energy showers are lost in the gap between detectors) and the quality of the shower sampling. The particular kind of detector (scintillator, RPC, water tank) determines the detail of measurement (efficiency, resolution, energy threshold, quality) and impact on the cost per detector *Cd*. In principle, best physics requires large area *A*, small distance *d* and high quality of the sampling. However, the cost of an array increases with *Cd* · *<sup>A</sup>*/*d*2, therefore a compromise is always needed. This is one of the reason why the typical total sensitive area of a classical array is less than 1% of the total enclosed area. This results in a high degree of uncertainty in the reconstruction due to sampling fluctuations which add to the shower fluctuations.

The experiments devoted to study the PeV energy range have been operated at different altitudes, ranging from the extreme altitude (5200 m asl) of BASJE-MAS [21] to the sea level of KASCADE [22–24].

In Tables 1 and 2, the characteristics of air shower arrays operated in the last two decades to study Galactic CR physics from ground are summarized. The atmospheric depths of the arrays, the main detectors used, the energy range investigated, the sensitive areas of e.m. and muon detectors, the instrumented areas and the coverage (i.e., the ratio between sensitive and instrumented areas) are reported. The depth in atmosphere is crucial to fix the energy threshold, the energy resolution, the impact of shower-to-shower fluctuations, then the sensitivity to elemental composition.


**Table 1.** Characteristics of different air shower arrays.

**Table 2.** Characteristics of different muon detectors operated in some shower arrays.


Generally speaking, near the depth of the maximum of the shower development, the number of secondary charged particles is almost independent of the mass of the primary particle, and the shower fluctuations are at minimum. For the knee energy region, this depth corresponds to ≈5000 m asl. Therefore, these extreme altitudes are suitable to have good energy resolution, to reconstruct the primary energy in a mass-independent way and to study the shower core region in great detail, where the hadronic component feeds the e.m. one deep in the atmosphere. As demonstrated by the ARGO-YBJ [30] and Tibet AS*γ* [31] experiments, observables related to the shower core properties are almost independent on the details of hadronic interaction models. At high altitudes, due the low energy threshold (≈TeV), it is possible to cross-check the fluxes with direct measurements on a wide energy range (ARGO-YBJ in the 5–250 TeV range). This cross-calibration is important due to the conflicting results obtained not only by ground-based detectors, but also by direct experiments. In addition, the absolute energy scale can be calibrated at a level of 10%, exploiting the so-called *"Moon Shadow"* technique [20].

On the other hand, experiments located deep in the atmosphere enhance the differences in the longitudinal development of EAS of different primary masses, as the shower is sampled well beyond its maximum. Therefore, the ratio N*e*/N*<sup>μ</sup>* is, in principle, more suitable for elemental composition studies. However, shower fluctuations are much larger, making it difficult to interpret the data. In addition, the reconstruction of the energy is typically strongly model-dependent because 'a priori' assumptions on the primary composition are needed, and the calibration of the absolute energy scale is one of the major open issues.

The great variety of layouts, observables, and reconstruction procedures to infer the elemental composition is at the origin, in part, of the conflicting results reported by different ground-based experiments. Arrays focused on the investigation of the knee region operated so far are also characterized by a limited size of the instrumented area. They collected limited statistics above 10<sup>16</sup> eV, and were, therefore, unable to give a conclusive answer to the origin of the knee. The poor sensitivity to elemental composition, due to the small statistics, prevents discrimination against different mass groups, and only general trends can be investigated in terms of the evolution of ln *A* or of *"light"* and *"heavy"* components with energy.

#### **3. Extensive Air Showers: The Heitler-Matthews Model**

A general idea of the main characteristics of EAS and of how different nuclei produce showers with different properties can be obtained from some relatively simple arguments, as suggested by Heitler [32] and Matthews [33]. This toy model is useful to show how different observables depend on the primary mass and energy, and why certain techniques have historically been used to study elemental composition or to reconstruct the energy spectrum. Nevertheless, detailed Monte Carlo simulations must be used to describe quantitatively all the characteristics of these random processes, with particular care to the role of shower fluctuations.

In a nutshell, the collision of a primary CR with a nucleus of the atmosphere produces one large nuclear fragment and many charged and neutral pions (with a smaller number of kaons) (Figure 3) [34]. A significant fraction of the total energy is carried away by a single *"leading"* particle. This energy is unavailable immediately for new particle production. Roughly speaking, half of the energy of the primary particle is transferred to the nuclear fragment and the other half is taken by the pions (and kaons). The fraction of energy transferred to the new shower particles is referred as *inelasticity*. Accurate description of the leading particles is crucial because these high-energy nucleons feed energy deeper into the atmosphere. Approximately equal number of positive, negative and neutral pions are produced. The e.m. component, the most intense of an EAS, is produced by the photons coming from the decay of the neutral pions. At each interaction before the charged pions decay, nearly a third of the hadronic component energy is released into the e.m. one.

As the number of particles increases, the energy per particle decreases. They will also scatter, losing energy, and many will range-out. Thus, the number of particles (or, with less ambiguities in the definition, the quantity of energy transferred to secondaries and eventually released into the atmosphere) will reach a maximum at some depth X*max* which is a function of energy, of the nature of the primary particle and of the details of the interactions of the particles in the cascade. After that, the energy/particle is so degraded (will be below some *"critical energy"*) that energy losses dominate over particle multiplication process, and the shower "size" will decrease as a function of depth: it grows 'old'. Once the pions have reached an energy which is low enough, they will decay into muons and neutrinos (*π*<sup>+</sup> <sup>→</sup> *<sup>μ</sup>*+*νμ* or *<sup>π</sup>*<sup>−</sup> <sup>→</sup> *<sup>μ</sup>*−*ν*¯*μ*). The resulting muons propagate unimpeded to the ground. The muon cascade grows and maximizes, but the decay is slower as a consequence of the relative stability of the muon and small energy losses by ionization and pair production.

**Figure 3.** Schematic evolution of cascades initiated by a CR particle. At each step, roughly 1/3 of the energy is transferred from the hadronic cascade to the e.m. one. Figure taken from [34].

These are the most common processes, but not at all the only ones. As an example, successive hadronic interactions of the primary CR, interactions/decays of kaons and muon decays, multiple scattering and production angles must also be considered (see, for example, the the Ref. [16]). Only detailed simulations with Monte Carlo methods are able to describe all the characteristics of these random processes.

Historically, one of the main problems in analyzing data from shower arrays was related to the fact that each experiment used its own simulation of shower development and detectors. This made difficult the comparison of the results and the understanding of their differences. Starting in the 1990s, all experiments began to use the same Monte Carlo simulation code CORSIKA [35], a framework containing different hadronic interaction models to describe the shower development in the atmosphere, and the software GEANT [36] to simulate the detectors operated in the arrays. Over the years, other simulation codes have been developed, in particular to describe the development of showers at ultra-high energies, such as AIRES [37]. The main characteristics of hadronic interactions that are relevant for EAS physics are: cross-sections (p–air, *π*–air, N–air), inelasticity of the collisions, multiplicity/composition of secondaries, transverse momentum distribution, fraction of diffractive dissociation.

New data coming from the LHC (at an energy *Elab*∼10<sup>15</sup> eV) allowed to improve the models even if some points remain critical. In fact, the situation is much worse than it may appear from energy considerations. Measurements at colliders are limited to an angular region that excludes the beam pipe (the so-called *'central region'*), and therefore a very large majority of the high-energy particles that are emitted at small angles (in the so-called *'forward region'*) are unobservable. In EAS physics the forward region is the most relevant because the high-energy particles feed energy in the shower down in the atmosphere. Therefore, models tuned to accelerator measurement in the central region are extrapolated to describe the interactions of CRs.

Nevertheless, this simple toy model predicts the basic features of EAS development. In the following, the e.m. and hadronic processes will be described separately in more detail.

#### *3.1. Electromagnetic Showers*

The main features of an e.m. shower profiles can be described within the simple Heitler's toy model of particle cascades [32]. Let us suppose that a particle (electron, positron or photon) with energy E0 splits its energy equally into two particles after traveling a radiation length X0 in air, and let this process be repeated by the secondaries (see Figure 4).

**Figure 4.** Schematic view of an e.m. cascades (**a**), and of a hadronic shower (**b**). In the hadron shower, dashed lines show *π*<sup>0</sup> which do not re-interact but decay, producing e.m. sub-showers.

Let *X* describe the depth in the atmosphere and define the depth at which the average CR starts interactions with the atmosphere to be *X* = 0 g/cm2. After *n* radiation lengths, we obtain a particle cascade which has evolved into *N* = 2*<sup>n</sup>* particles of equal energy *E* = *E*0/*N*. Multiplication stops when the energies of the particles are too low for pair production or bremssthralung. This energy is the critical energy *εem <sup>c</sup>* in the air (≈80 MeV, below which the collisional energy losses are dominant).

The maximum number of particles *Nmax* is reached at this moment, when all particles have the same energy *εem <sup>c</sup>* , *E*<sup>0</sup> = *εem <sup>c</sup>* · *Nmax*. The depth *Xmax* at which the shower reaches the maximum size is *Xmax* = *nmax* · *X*0, where *nmax* is the number of radiation lengths required for the primary energy to be reduced to *εem c* .

Since *Nmax* = 2*nmax* , we have

$$m\_{\max} = \ln\left(\frac{E\_0}{\varepsilon\_c^{em}}\right) \cdot \frac{1}{\ln 2} \tag{2}$$

so that

$$X\_{\text{max}}^{\text{em}} = \frac{X\_0}{\ln 2} \cdot \ln \left(\frac{E\_0}{\varepsilon\_c^{\text{em}}}\right). \tag{3}$$

Finally, it is interesting to estimate the *elongation rate* Λ, that is, the rate of increase of *Xmax* with the primary energy. From the relation (3), we have

$$
\Lambda^{em} = \frac{d\,X\_{\text{max}}}{d\,\log\_{10}E\_0} = 2.3 \cdot X\_0 = 85 \text{ g/cm}^2 \qquad \text{per decade of energy.} \tag{4}
$$

This simple model predicts two basic features of e.m. shower development:


#### *3.2. Hadronic Showers*

Air showers initiated by protons have been modeled by different authors (see, for example, the the Ref. [38–40]) following the Matthews approach [33], similar to the Heitler one. The main differences with the e.m. cascades are


According to the Ref. [33], constant values *Nch* = 10, corresponding to an energy of about 100 GeV, and *ε<sup>π</sup> <sup>c</sup>* = 20 GeV are adopted in the following.

Protons travel one interaction length and interact producing *Ntot* pions, all having equal energies, *Nch* are charged and *Ntot* <sup>3</sup> <sup>=</sup> <sup>1</sup> <sup>2</sup> · *Nch* neutral, which immediately decay into photons, initiating e.m. showers. As for the e.m. cascade, we assume equal division of energy during particle production.

In turn, the charged pions can decay in muons and neutrinos and hence, as long as their decay length remains larger than their interaction length, they will re-interact rather than decay. This happens for *γcτπ* > *λπ* <sup>−</sup>air/*ρ*air, with the Lorentz factor *γ* = *Eπ*/*mπ*, the charged pion lifetime *τπ* 26 ns and *λπ* <sup>−</sup>air 1.5*λ<sup>p</sup>* <sup>−</sup>air 120 g cm−<sup>2</sup> (since the *<sup>π</sup><sup>p</sup>* cross-section is about 2/3 the *pp* cross-section). This implies that pions will re-interact as long as their energy satisfies *<sup>E</sup>* > *<sup>E</sup>*<sup>d</sup> 100 GeV(10−<sup>4</sup> g cm−3/*ρ*air). Hence, at the heights above 10 km, where the initial development of the shower takes place, *π*± will re-interact for energies greater than ∼20–30 GeV [40].

After *n* interactions, the *N<sup>π</sup>* = (*Nch*)*<sup>n</sup>* charged pions produced carry a total energy of ( 2 <sup>3</sup> )*<sup>n</sup>* · *<sup>E</sup>*0. The energy per charged pion after *<sup>n</sup>* interactions is then *<sup>E</sup>π*<sup>±</sup> <sup>=</sup> *<sup>E</sup>*<sup>0</sup> (3/2*Nch*)*<sup>n</sup>* . The remainder of the primary energy goes into the e.m. component from *π*<sup>0</sup> decays

$$E\_{\rm em} \simeq E\_0 \left[ 1 - \left( \frac{2}{3} \right)^n \right]. \tag{5}$$

After only six interactions, about 90% of the initial energy is transferred to the e.m. component of the shower, with the remaining 10% being essentially the muons and neutrinos from the charged pion decays. As a consequence, most of the energy of an air shower can be observed in its e.m. component. This is the so-called *calorimetric energy* which allows to estimate the primary energy with good accuracy to detectors able to observe the longitudinal air shower development.

Assuming that at *ε<sup>π</sup> <sup>c</sup>* , all pions decay, the number of muons is *N<sup>μ</sup>* = *Nπ*<sup>±</sup> = (*Nch*)*nc* , where *nc* is the number of interaction lengths required for the charged pion's interaction length to exceed its decay length

$$m\_{\varepsilon} = \frac{\ln(E\_0/\varepsilon\_{\varepsilon}^{\pi})}{\ln(\frac{3}{2}N\_{ch})} = 0.85 \lg\left(\frac{E\_0}{\varepsilon\_{\varepsilon}^{\pi}}\right). \tag{6}$$

Thus, the total energy is divided into two channels, hadronic and electromagnetic

$$E\_0 = E\_{cm} + E\_h = \varepsilon\_c^{cm} \cdot N\_c + \varepsilon\_c^{\pi} \cdot N\_\mu. \tag{7}$$

This equation represents *energy conservation*, apart from a fraction of a few percent of the primary energy spent in the neutrino component. The relative magnitude of the contribution from *Nμ* and *Ne* does not depend on the details of the model, but only on the respective critical energies, the energy scales at which e.m. and hadronic multiplication ceases. An important conclusion of this description of the hadronic cascades is that the energy is given by a linear combination of muon and electron sizes. This result is insensitive to fluctuations in the division of energy between the hadronic and e.m. channels and independent on the mass of the primary particle.

The number of muons is given by

$$\ln N\_{\text{fl}} = \ln N\_{\pi^{\pm}} = n\_{\text{c}} \ln N\_{\text{cl}} = \frac{\ln(E\_0 / \varepsilon\_{\text{c}}^{\pi})}{\ln(3/2 N\_{\text{cl}})} \cdot \ln \left(N\_{\text{cl}}\right) = \beta \cdot \ln \left(\frac{E\_0}{\varepsilon\_{\text{c}}^{\pi}}\right). \tag{8}$$

Following [33], we can estimate *β* = ln (*Nch*) ln (3/2*Nch*) <sup>=</sup> 0.85 for *<sup>E</sup>*<sup>0</sup> in the range 1014–1017 eV, obtaining

$$N\_{\mu} = \left(\frac{E\_0}{\varepsilon\_c^{\pi}}\right)^{\beta} = \left(\frac{E\_0}{\varepsilon\_c^{\pi}}\right)^{0.85} \sim 9900 \left(\frac{E\_0}{10^{15} \text{eV}}\right)^{0.85}.\tag{9}$$

Including inelasticity in the Heitler model [33] changes the parameter *β*

$$\beta = \frac{\ln\left(N\_{\rm ch}\right)}{\ln\left(3/2N\_{\rm ch}\right)} \to \frac{\ln\left[1+N\_{\rm ch}\right]}{\ln\left[\left(1+N\_{\rm ch}\right)/\left(1-\frac{1}{3}\kappa\right)\right]} \approx 1 - \frac{\kappa}{3\ln\left(N\_{\rm ch}\right)} = 1 - 0.14\kappa. \tag{10}$$

The elasticity for the most energetic meson in pion–air interactions yields (1 − *κ*) between 0.26 and 0.32, resulting in *β* = 0.90.

The electronic size can be calculated by inserting the expression (9) for the muon size in the energy conservation relation (7)

$$\frac{E\_{\rm em}}{E\_0} = \frac{E\_0 - N\_{\mu} \varepsilon\_{\mathcal{E}}^{\pi}}{E\_0} = 1 - \left(\frac{E\_0}{\varepsilon\_{\mathcal{E}}^{\pi}}\right)^{\beta - 1}.\tag{11}$$

The e.m. fraction is 66% at *E*<sup>0</sup> = 1015 eV, increasing to 83% at 10<sup>18</sup> eV for protoninduced showers.

Therefore, the number of electrons at a maximum shower for proton-induced showers is

$$N\_{\varepsilon} = \frac{E\_{\varepsilon m}}{\varepsilon\_{\varepsilon}^{cm}} = \frac{E\_0}{\varepsilon\_{\varepsilon}^{cm}} - \frac{\varepsilon\_{\varepsilon}^{\pi}}{\varepsilon\_{\varepsilon}^{cm}} \left(\frac{E\_0}{\varepsilon\_{\varepsilon}^{\pi}}\right)^{\beta} \approx \frac{E\_0}{\varepsilon\_{\varepsilon}^{\epsilon m}} = N\_{\varepsilon|\_{\max}}^p. \tag{12}$$

The approximation is justified at high energies when the fraction of energy transferred to muons is small [39].

In the framework of the *superposition model*, each nucleus is taken to be equal to *A* individual single nucleons, each with energy *E*0/*A* and each acting independently. The shower resulting from the interaction of the primary nucleus *A* can be treated as the sum of *A* proton-induced independent showers all starting at the same point. Thus, while a proton creates one shower with energy *E*0, an iron nucleus of the same total energy is expected to create the equivalent of 56 proton showers, each with reduced energy (*E*0/56). The average properties of showers are well reproduced by this model, though the fluctuations are clearly underestimated and can be studied only with detailed Monte Carlo simulations of the intra-nuclear cascade. The superposition of *A* independent showers naturally explains why the shower-to-shower fluctuations are smaller for shower initiated by nuclei as compared to proton showers.

By substituting the lower primary energy (*E*0/*A*) into the previous expressions and summing *A* such showers, we obtain the following relations for the number of electrons and muons in a shower induced by a nucleus *A*:

$$N\_{c|\_{max}}^A = A \left(\frac{E\_0/A}{\varepsilon\_c^{cm}}\right) = N\_{c|\_{max}}^p \tag{13}$$

$$N^{A}\_{\mu|\_{\rm tur}} = \left(\frac{E\_0}{\varepsilon\_c^{\pi}}\right)^{\beta} A^{1-\beta} = N^{p}\_{\mu|\_{\rm tur}} A^{1-\beta} \approx 1.69 \cdot 10^4 \cdot A^{0.10} \left(\frac{E\_0}{1 \text{ PeV}}\right)^{0.90} \cdot \tag{14}$$

From these relations valid *at shower maximum* follows:


A large number of ground-based arrays studying the knee energy region are located deep in the atmosphere and do not sample the number of electrons at shower maximum. Therefore, the experimental situation is not ideal because the size, used to recover the energy of the primary particle, is mass-dependent, as discussed in Section 4. Only experiments located at extreme altitude (above 4000 m asl) observe the electrons in the shower maximum region for near-vertical showers with an energy in the PeV range.

Deeper in the atmosphere, arrays measure only the attenuated size

$$N\_{\varepsilon|\_{graul}} \approx N\_{\varepsilon|\_{max}} \cdot \exp\left(-\frac{\Delta X}{\Lambda}\right) \tag{15}$$

where <sup>Δ</sup>*<sup>X</sup>* is the distance of the shower maximum from the ground and <sup>Λ</sup> ≈ 60 g/cm<sup>2</sup> is the attenuation length of the electron size after the shower maximum. Since heavy nuclei reach the maximum of longitudinal development at smaller depths than light ones, on the ground we have a larger electron number for air showers initiated by light particles. This implies that, due to the steeply falling CR spectrum, showers of equal ln *Ne* are enriched in light elements.

#### *3.3. Longitudinal Development*

The longitudinal development of a hadronic shower is dominated by the parallel e.m. sub-showers produced by the *π*<sup>0</sup> decays in the first interaction, at an atmospheric depth *<sup>X</sup>*<sup>∗</sup> <sup>=</sup> *<sup>λ</sup>p*−*air* · ln <sup>2</sup> <sup>≈</sup> 55 g/cm2. In a good approximation following cascades can be neglected. The energy of the single photon is *<sup>E</sup><sup>γ</sup>* <sup>=</sup> *<sup>E</sup>π*<sup>0</sup> <sup>2</sup> <sup>=</sup> *<sup>E</sup>*<sup>0</sup> <sup>3</sup> <sup>2</sup> *Nch* <sup>=</sup> *<sup>E</sup>*<sup>0</sup> <sup>3</sup>*Nch* .

From Equation (3), we have

$$\begin{split} X\_{\text{max}}^p &= X^\* + X\_0 \cdot \ln \left( \frac{E\_0}{3N\_{\text{cl}} \cdot \mathbf{e}\_c^{\text{cm}}} \right) \\ &= X^\* + X\_{\text{max}}^{\text{cm}} - X\_0 \cdot \ln(3N\_{\text{cl}}) \quad \text{g/cm}^2 \end{split}$$

where *Xem max* is the atmospheric depth of the maximum of *γ*-induced showers with *E*<sup>0</sup> primary energy and *Nch* is the multiplicity of charged pions in the first interaction. The elongation rate for showers induced by protons is then

$$
\Lambda^p = \Lambda^\gamma + \frac{d}{d\log\_{10} E\_0} \left[ X^\* - X\_0 \cdot \ln(\text{3 } N\_{\text{ch}}) \right] = 58 \quad \text{g/cm}^2 \text{ per decade}, \tag{16}
$$

reduced from the elongation rate for purely e.m. showers. This estimation verifies Linsley's elongation rate theorem [41], which points out that e.m. showers represent an upper limit to the elongation rate of the hadronic showers. The shower maximum is expected to be influenced by the elasticity of the first interaction, (1 − *κ*) = *Elead*/*E*0, where *Elead* is the energy of the leading particle. For interactions with (1 − *κ*) > 0.5 most of the primary energy will be transferred deeper into the atmosphere and correspondingly the shower maximum will be deeper.

The extrapolation to a primary particle with mass *A* with the superposition model yields

$$X\_{\max}^{A} = X\_{\max}^{p} - X\_0 \cdot \ln A \tag{17}$$

Detectors able to observe the longitudinal air shower development can estimate the primary energy with good accuracy measuring the so-called *calorimetric energy*, that is, the energy of the e.m. component. With this estimator of the energy of the primary particle, the orthogonal variable sensitive to its primary mass is the depth of the shower maximum in terms of the number of particles, *Xmax*.

Therefore:


Despite the simplicity and the approximations of the toy model, the main characteristics of the EAS development are quite well reproduced. Obviously, a detailed description of the cascade, in particular for what concern the role of fluctuations, can be provided only by detailed Monte Carlo simulations.

#### *3.4. Energy and Mass*

The relevance of muon measurements to the question of the primary composition has been first remarked by the Institute for Nuclear Studies (INS) group in Tokyo [42]. They were the first group to point out the key information that the mass of the primary particle could be derived from a study of plots of muon versus electron number.

Due to the intuitive relation between shower to shower fluctuations and primary mass, the study of fluctuations in the muon number distributions was historically the first method employed to study the primary CR mass composition [42–44]. The narrowing of the distribution of *Nμ*/*Ne* was considered to be due partly to the change in the composition of primary particles with energy [43].

On general grounds, the elemental composition can be investigated if the total size and the muon component depend differently from the primary energy. If we assume that their dependences from the energy of a primary proton *E*<sup>0</sup> can be described as

$$\mathcal{N}\_{\mathfrak{C}} \propto E\_0^{\beta\_{\mathfrak{F}}} \, \prime \qquad \qquad \mathcal{N}\_{\mathfrak{V}} \propto E\_0^{\beta\_{\mathfrak{V}}} \, \prime \tag{18}$$

for a nucleus *A*, we have

$$N\_c \propto \left(\frac{E\_0}{A}\right)^{\beta\_\varepsilon} \cdot A\_\prime \qquad\qquad\qquad N\_\mu \propto \left(\frac{E\_0}{A}\right)^{\beta\_\mu} \cdot A\_\prime \tag{19}$$

with a mass-number dependency of the type 1 − *β<sup>e</sup>* and 1 − *βμ*, respectively. The relation *Nμ*/*Ne* can be easily deduced

$$N\_{\mu} \propto N\_{\mathfrak{e}}^{\mathbb{S}\_{\mu}/\mathbb{R}\_{\mathfrak{e}}} A^{1 - (\mathbb{S}\_{\mu}/\mathbb{R}\_{\mathfrak{e}})}.\tag{20}$$

In 1962, Linsley, Scarsi, and Rossi working at the MIT Volcano Ranch Station observed, for the first time, a muon/electron correlation: *<sup>N</sup>μ*∼*A*1−*<sup>α</sup>* · (*Ne*)*α*, thus establishing that the muon size is a mass-sensitive observable [45].

The Equation (13) can be transformed to obtain the energy *E*<sup>0</sup> to be introduced in the relation (14) to obtain *βμ*/*βe*∼0.86. The muon size for a given mass *A* as a function of the total size *Ne* is then

$$N\_{\mu} \propto N\_{c}^{0.86} A^{0.14} . \tag{21}$$

In a similar way, we can obtain the muon size as a function of the total size for a given primary energy *E*0. The Equation (13) is transformed to obtain the mass *A* which in turn is introduced in the relation (14)

$$N\_{\mu} \propto \left(\frac{E\_0}{1\text{ PeV}}\right)^{3.17} N\_c^{-2.17}.\tag{22}$$

In experiments with ground-based arrays the reconstructed number of muons and electrons are plotted in a ln *N<sup>μ</sup>* − ln *Ne* plane to recover the energy and mass of the primary particle. This diagram, when combined with detailed shower simulations, proved to be a powerful tool for extracting information on primary mass.

Therefore, it is interesting to study the electron-to-muon ratio at shower maximum

$$\frac{N\_{\text{f}}}{N\_{\text{fl}}} \approx 35.1 \cdot \left(\frac{E\_0}{A}\right)^{0.15}.\tag{2.3}$$

with the energy in PeV. This ratio depends on the energy per nucleon *E*0/*A* of the primary particle, thus showing that *Ne*/*N<sup>μ</sup>* can be used to infer the mass of the primary particle if the energy is measured with a different, independent observable.

We can use the relation (23) to investigate the sensitivity of EAS arrays to the primary mass A [38,46]

$$\log\left(\frac{N\_{\varepsilon}}{N\_{\mu}}\right) = 1.54 + 0.15 \cdot \lg\left(\frac{E\_0}{1 \text{ PeV}}\right) - 0.065 \cdot \ln A = C - 0.065 \cdot \ln A \tag{24}$$

If the energy is reconstructed from another independent observable, the mass of the primary CR can be determined by measuring the ratio *Ne*/*Nμ*. Therefore, the relative error on the electron-to-muon ratio is

$$\frac{\Delta(N\_{\rm c}/N\_{\rm \mu})}{N\_{\rm c}/N\_{\rm \mu}} \sim 0.15 \left[ \frac{\Delta E\_0}{E\_0} + \frac{\Delta A}{A} \right] \sim 0.15 \left[ \frac{\Delta A}{A} \right] \tag{25}$$

with the consequence that to measure the elemental composition with a resolution of one unit in ln *A* the relative error on *Ne*/*N<sup>μ</sup>* must be ≈15%. A resolution of one unit in ln *A* in principle allows to reconstruct 4 (or 5 ?) different mass groups: p, He, CNO, MgSi (?) and Fe. The large shower-to-shower fluctuations often only allow one to trace the light and heavy components or the parameter ln *A* with energy. Similarly, from the relation (17) follows that the position of the shower maximum must be determined with a resolution of about one radiation length X0∼37 g/cm2 to have a resolution in ln *<sup>A</sup>* of one unit.

#### **4. Reconstruction of the Energy and Mass of the Primary Particle**

The crucial point in air shower observations with EAS arrays is the reconstruction of the primary particle properties (especially energy and mass number) from the measured quantities. In fact, analysis of shower data consists in the disentanglement of a threefold problem involving primary energy, primary mass and modelling of hadronic interactions (for a discussion about hadronic interactions in CR physics see, as an example, refs. [47,48] and references therein). An intrinsic ambiguity affects the interpretation of data. Different combinations of the two following elements can produce similar showers. As an example, a *"short"* shower can be produced by a large cross-section, high inelasticity or heavy primary mass. On the contrary, a *"long"* shower, penetrating deeper in the atmosphere, can be produced by small cross-section, low inelasticity or light mass.

(1) shower development, mainly governed by the inelasticity and by the inelastic cross section (2) elemental composition of the primary flux, that we don't know and want to measure

Strictly speaking, when operating with shower arrays there are no observables directly related to the mass of the primary particle, and its measurement is very indirect. We note, however, that the Cherenkov light emitted by a primary heavy nucleus high up in the atmosphere (the so-called *"direct Cherenkov light"*) is directly related to the charge (and therefore to the mass) of the primary particle [49,50]. Since this light is proportional to *Z*2, heavy nuclei are more suited for detection. Charge resolution is about 10% for Z > 10. The main limitation is that it can only be used over a small energy range for each atomic charge Z.

The majority of experiments with shower arrays can therefore apply only 'statistical' methods according to a classical scheme:


Therefore, a typical data analysis consists in finding a combination of primary energy spectrum, elemental composition and hadronic interaction characteristics to obtain a consistent description of the experimental results. Clearly, this is not a measurement, but only a consistency check of some trial models. In case of discrepancy, it is difficult to identify the origin; in case of agreement, is the parameter combination unique?

Due to the reduced resolution in the measurement of the primary mass (see Section 3), the majority of shower arrays displayed the results only as a function of the total energy per particle with the so-called *"all-particle"* energy spectrum, that is, as a function of the total energy per nucleus, and not per nucleon. Any tentative to infere informations about elemental composition are limited, at most, to study the evolution of the *"light"* ("protonlike") or *"heavy"* ("iron-like") components as a function of the energy, with results which critically depend on Monte Carlo predictions.

In the last two decades, a number of multi-component experiments have started to measure, with high statistics, at the same time, different shower observables, on an event-by-event basis. This fact allowed to exploit sophisticated analysis techniques to infer the characteristics of the primary particle by measuring the correlation between different components (for a review see, for example, the Refs. [39,51] and references therein).

In a nutshell,

• *How to obtain the energy spectrum in shower arrays?*

This is the first step in the analysis of CR data. We measure the spectrum in one observable and make a conversion to the energy spectrum. The observable typically used is the shower size because, as discussed in Section 3.2, the number of electrons at shower maximum is nearly independent on the primary mass: *N<sup>A</sup> <sup>e</sup>*|*max* <sup>≈</sup> *<sup>N</sup><sup>p</sup> <sup>e</sup>*|*max* . However, surface detectors are usually located deep in the atmosphere and do not measure the number of electrons at shower maximum. Beyond the maximum, the number of electrons is a mass-sensitive parameter, with a larger electron number for air showers initiated by light primaries, according to a relation of the type

$$N\_{\mathcal{E}}(E, A) = \mathfrak{a}(A) \cdot E^{\mathfrak{f}(A)} \tag{26}$$

where the parameters *α* and *β* depend on the primary mass *A*. This implies a degeneracy in the reconstruction procedure because to recover the primary energy from the size spectrum we must assume a given elemental composition to be measured. If the

composition changes in the investigated energy range, the relationship between the measured electron size and inferred energy will also vary. The number of electrons in the core region has been used in some experiment, as well as the particle density at a suitable given distance from the shower axis, in some large arrays (see, for example, the Refs. [11,52]). In both cases this densities, according to Monte Carlo simulations, are nearly independent of the primary mass.

• *How do we measure elemental composition at ground?*

The inelastic cross-section *σFe*−*Air inel* of iron at 1 PeV is about six times larger than for protons of equal energy. Hence, nuclei develop showers higher in atmosphere (smaller *Xmax*) than protons, dissipating their energy much faster. Due to the shorter interaction length and the smaller energy per nucleon and because of the reduced attenuation of the muon component, nuclei-induced showers contain less particles in the e.m. component deep in the atmosphere, but they carry more muons than a proton shower of the same energy. This is the basis of the electron-muon correlation method. Therefore, as discussed in Section 3.4, the measurement of electron and muon contents simultaneously (with their fluctuations) has become the first and most commonly employed technique to infer the CR elemental composition with arrays. However, intrinsic shower to shower fluctuations limit mass resolution to a few mass groups (see Section 3.4) and electron and muon numbers are not independent. In addition, the muon component is heavy dependent on the details of the hadronic interactions and the results strongly depend on the particular model used to interpret the data.

The other common technique, below 10<sup>18</sup> eV, involves the observation of the Cherenkov light and the study of its shape. In fact, the characteristics of the photon distribution depend on the depth of the shower maximum, therefore on the mass of the primary particle. The overall Cherenkov intensity provides a calorimetric measurement of the CR energy. Cherenkov light has been measured, for instance, in hybrid experiments by ARGO-YBJ and Tunka apparatus.

The KASCADE multi-component array was the first experiment that claimed the measurement of the energy spectra of 5 different mass groups (p, He, CNO, MgSi, Fe) through a complex unfolding of the *Ne*/*N<sup>μ</sup>* diagram [22–24]. In the last two decades, other multi-component experiments measured a number of observables that, in principle, are mass-sensitive: steepness of the lateral distribution, characteristics of shower core region, distribution of the relative arrival times and angles of incidence of the muon component, characteristics of the lateral distribution of high energy muons (the so-called *"muon bundles"*) measured underground, pulse shape and lateral distribution of the air Cherenkov light, depth X*max* of the shower maximum (see, for example, ref. [15,17]). But the study of the muon component has remained the most used technique.

#### **5. Elemental Composition in the 1014 to 1018 eV Region**

Several experimental results associate the knee with the bending of the light component (p and He), and are compatible with a rigidity-dependent cut-off [22–24,53–55]. However, the flux of the different components vary significantly depending on the interaction model used to interpret the data [22–24]. On the contrary, other results (in particular those obtained by arrays located at high altitudes) seem to indicate that the knee of the all-particle energy spectrum is due to heavier nuclei and that the light component cuts off well below 1 PeV [21,26,31,52,53,56].

In this section the measurements of the light component energy spectrum in the 10<sup>14</sup> to 1018 eV region will be presented by using different plots to point out the conflicting results between the experiments.

In Figure 5 the energy spectra of the light component as measured by Tibet AS*γ* [25,31] and ARGO-YBJ [52] are shown. Both experiments are located in the YangBaJing Cosmic Ray Laboratory in Tibet (China) at 4300 m a.s.l. and did not exploit a measurement of the muon component to determine the elemental composition of the primary CR flux.

**Figure 5.** Energy spectra of the light (p+He) component as measured by Tibet AS*γ* [25,31] and ARGO–YBJ [52] experiments with different techniques and analyses, compared with results obtained in direct observations by CREAM [57] and NUCLEON [58].

The Tibet AS*γ* Collaboration reconstructed the energy spectrum studying the shower core region with a burst detector as well as with emulsion chambers. The ARGO-YBJ experiment measured the CR energy spectra exploiting completely different and independent approaches [52]:


All the results are in excellent agreement. In the ARGO-YBJ experiment the selection of (p+He)-originated showers is performed not by means of an unfolding procedure after the measurement of electronic and muonic sizes, but on an event-by-event basis exploiting showers topology, that is, the lateral distribution of charged secondary particles. This approach is made possible by the full coverage of the central carpet, the high segmentation of the read-out and the high altitude location of the experiment that retains the characteristics of showers lateral distribution in the core region. The contamination of nuclei heavier than helium is estimated smaller than 15% at 1 PeV in all analyses.

In Figure 5 the direct measurements reported by CREAM [57] and NUCLEON [58] are also shown. ARGO-YBJ is the only experiment that traced the (p+He) spectrum across the knee starting from an energy so low (≈TeV) to overlap with direct measurements and to cross-calibrate the fluxes on a wide energy range (5–250 TeV). These results show that, when indirect measurements are capable of selecting almost pure beams, their findings are in fair agreement with direct ones and confirm that current simulation models provide a satisfactory description of the EAS development in the atmosphere. The cross-calibration of fluxes in this energy range, where the boundary line between 'direct' and 'indirect' measurements is uncertain, is very important. The low energy threshold allowed also a calibration of the absolute energy scale at a level of 10% exploiting the *Moon Shadow* technique in the 1–30 TeV/Z range [20].

As can be seen from the figure, the observations of Tibet AS*γ* and ARGO-YBJ are in good agreement each other showing that the knee of the (p+He) energy spectrum is at ≈500–700 TeV, well below the energy of knee in the all-particle spectrum. Similar conclusions have been obtained by the BASJE-MAX experiment located at 5200 m asl [21] and by EAS-TOP at 2000 m asl [53] and by CASA-MIA at 1450 m asl [26].

In Figure 6 the energy spectra of the light component reconstructed by the KASCADE experiment [22–24] with two different hadronic interaction models are added for comparison. The energy threshold is about 1 PeV and the experiment was located at sea level. KASCADE did use of a complex unfolding procedure to recover the elemental composition from the *Ne*−*N<sup>μ</sup>* diagram in terms of 5 mass groups (p, He, CNO, MgSi, Fe). As can be seen from the figure, both the spectra are at variance with the results obtained by Tibet AS*γ* and ARGO-YBJ, suggesting that the knee of the CR all-particle spectrum at a few PeV is due to the bending of the light component.

**Figure 6.** The energy spectra shown in Figure 5 compared with the results obtained by the KASCADE experiment by using two different interaction models to interpret data [22–24].

All measurements of the light component up to about 10<sup>18</sup> eV are summarized in the Figure 7 and compared with the parametrization provided by Horandel [62]. Roughly speaking, we can separate the experiments that measured the (p+He) energy spectrum in the PeV range in 2 different groups


Measurements exploiting the longitudinal profile of the showers with Cherenkov detectors are conflicting too. The results obtained with the ARGO-YBJ hybrid detector (carpet and Cherenkov telescope) are in agreement with those of the carpet only, whereas the observations of Cherenkov light by Tunka-133 are consistent with KASCADE and KASCADE-Grande findings.

**Figure 7.** The energy spectra shown in Figure 6 compared with the results obtained by HAWC in the 10–100 TeV region [63] and by KASCADE–Grande [27] and IceTop/IceCube [29] combined above the PeV. The parametrization of the light component provided by Hörandel [62] is also shown.

In Figure 8 the all-particle energy spectrum measured by several experiments in the energy region 1016–1018 eV is shown. As can be seen, the spectrum cannot be fitted by a single power law. We observe a spectral hardening at ∼<sup>2</sup> × 1016 eV and a steepening at ∼1017 eV. This result was first pointed out by KASCADE-Grande experiment [64,65], then more firmly assessed with higher statistics and precision by Tunka-133 and IceTop-73 [28,29], in particular for the feature at ∼<sup>2</sup> × 1016 eV.

The *light* (p+He) and *heavy* (C-Fe group) components measured by KASCADE-Grande are also shown. A knee is observed in the heavy component of CRs at *E* = 1016.92±0.04 eV, which coincides within the uncertainties with the change of the slope in the all-particle energy spectrum around 1017 eV. The spectral index changes from −2.76 ± 0.02 below the knee to −3.24 ± 0.05 above. At slightly higher energies (*<sup>E</sup>* = 1017.08±0.09 eV), the light component shows a hardening of the slope, with the spectral index changing from −3.25 ± 0.06 below this ankle-like structure to −2.79 ± 0.09 above. The positions of the changes of the slope as well as the intensities of the different components depend on the interaction model adopted to interpret the data. The knee in the heavy component seems visible also in the all-particle spectrum, as it is the dominant component.

The results obtained by Tunka-133, measuring the Cherenkov light deep in the atmosphere, suggest that the mass composition becomes heavier in the energy range <sup>10</sup>16–3 × <sup>10</sup><sup>16</sup> eV, then stays heavy till 10<sup>17</sup> eV, where the composition starts becoming lighter. IceTop/IceCube, exploiting the high-energy muons underground, indicates an increase of ln *<sup>A</sup>* in the energy range 1016–3 <sup>×</sup> <sup>10</sup><sup>17</sup> eV [65]. However, as discussed in the previous sections, the average logarithmic mass of CR ln *A* is used to describe the evolution of the composition as a function of energy when the mass resolution of the experiments do not allow a discrimination between different mass groups.

**Figure 8.** The all–particle energy range in the 'transition' region measured by different experiments. The *light* and *heavy* components measured by KASCADE–Grande are also shown.

In Figure 8 recent measurements of the all-particle energy spectrum down to about 100 PeV by the Pierre Auger Observatory are also reported. These observations suggest that the second knee is not a sharp feature, but a softening that extends in the interval 100–200 PeV [66]. Preliminary results based on the distribution of the depth of the shower maximum are consistent with a spectrum dominated by heavy nuclei in the 10<sup>17</sup> eV range becoming lighter with increasing energy [67].

The findings of Tunka-133, IceTop/IceCube and Auger are qualitatively in agreement with KASCADE-Grande. Despite the large uncertainty in the absolute composition, a common general trend is reported, composition gets heavier through the knee region and becomes lighter approaching the ankle.

In conclusion, the observations of the different ground-based arrays show two conflicting results regarding the maximum energy at which the light component is accelerated in CR sources. This disagreement is summarized in Figure 9 where the ARGO-YBJ results are compared with the KASCADE-Grande light and heavy spectra. These results cannot be reconciled and show the existence of a still unknown systematic uncertainty that, as discussed in previous sections, could be due to the different array characteristics (altitude, coverage), the observables used (muons or shower core characteristics), and the dependence on hadronic interaction models. These are certainly among the major sources of systematic errors that affect the interpretation of shower array measurements (for a recent discussion see, for instance, refs. [68,69]), although recent re-analyses of the KASCADE-Grande data with the latest versions of the post-LHC codes confirm previous results [70].

Important information could be deduced, in principle, by the measurement of the flux of atmospheric neutrinos, sensitive to the spectrum of parent CRs. In particular, the tail of the spectrum of atmospheric neutrinos is mainly shaped by the parent protons rather than by heavier element. As a consequence, we expect different predictions for the flux of atmospheric neutrinos according to ARGO-YBJ and KASCADE proton energy spectra, predictions that, in principle, can be checked at energies *E<sup>ν</sup>* ≥ 100 TeV if the atmospheric neutrinos could be properly identified. Unfortunately, in this energy region the total neutrino flux detected by IceCube departs from the existing predictions for atmospheric neutrinos suggesting the onset of an astrophysical component. The origin of such neutrinos is still unknown and current experimental uncertainties do not allow to draw clear conclusions [71].

**Figure 9.** The energy spectra of the light component measured by ARGO–YBJ compared to the *light* and *heavy* components measured by KASCADE–Grande.

#### **6. What's Next**

The experimental situation in the 100 TeV–100 PeV energy region must be clarified to solve the longstanding problem of the origin of the knee and to give solid foundations to CR models up to the highest observed energies. A new experiment, able to measure, at the right altitude and with high statistics, the elemental composition exploiting the techniques used so far in different apparatus, is mandatory to investigate the unknown uncertainties affecting the results so far obtained by shower arrays.

The only experiment that meets these requirement is LHAASO, a new multi-component array developed starting from the experience of the high altitude experiment ARGO-YBJ. The apparatus is located at high altitude (4410 m asl, 600 g/cm2) in the Daochen site, Sichuan province, P.R. China. LHAASO is expected to measure the energy spectrum, the elemental composition and the anisotropy of CRs in the energy range between 10<sup>12</sup> and 1017 eV [11,12,72,73]. The experiment is constituted by a 1 km<sup>2</sup> dense array of plastic scintillators and muon detectors. At the center of the array a 300 × 300 m<sup>2</sup> water Cherenkov facility will allow the detection of TeV showers. An array of 18 wide field of view Cherenkov telescopes will image the longitudinal profile of events. Neutron monitors will study the hadronic component in the core of air showers. LHAASO will study CR physics with different detectors and techniques starting from the TeV range, thus overlapping direct measurements in a wide interval. In Tables 1 and 2 the characteristics of the LHAASO-KM2A array are compared with other experiments. As can be seen, LHAASO will operate with a coverage of ∼0.5% over a 1 km2 area. The sensitive area of muon detectors is unprecedented (more than 40,000 m2), about 17 times larger than the CASA-MIA experiment, with a coverage of about 5% over 1 km2. For the first time the *Ne*/*N<sup>μ</sup>* correlation will be studied at high altitude with high statistics. This suite of independent instruments will also allow a deep investigation of the characteristics of the hadronic interaction models. The capability of hybrid measurements with Cherenkov telescopes operated in combination with a shower array have been demonstrated by the ARGO-YBJ measurement of the light component energy spectrum.

In addition, LHAASO will act simultaneously as a wide aperture (∼2 sr), continuouslyoperated gamma-ray telescope in the energy range between 10<sup>11</sup> and 1015 eV. The first results obtained during the first year of data taking with only a portion of the apparatus opened for the first time the PeV sky to observations, showing that the Northern hemisphere contains a lot of galactic PeVatrons.

Other projects under way to investigate, with a much higher energy threshold, the high energy tail of the galactic spectrum and the transition region are HiSCORE [74] and GRAND [75].

#### **7. Conclusions**

The results obtained by different experiments in the 10<sup>14</sup> to 10<sup>18</sup> eV region can be summarized as follows:

	- 1. All experiments observe an all-particle knee at ≈<sup>4</sup> × <sup>10</sup><sup>15</sup> eV.
	- 2. The absolute fluxes are in good agreement with each other and with the direct measurements.
	- 3. The elemental composition is conflicting. Experiments located at high and extreme altitude (BASJE-MAS, Tibet AS*γ*, ARGO-YBJ, EAS-TOP and CASA-MIA) reported evidence that the knee of the (p+He) component is below 1 PeV and that the composition at the all-particle knee energy is dominated by heavier nuclei. Experiments located deeper in atmosphere (KASCADE, KASCADE-Grande, IceTop/IceCube, Tunka-133) reported evidence that the proton knee is at the same energy of the all-particle knee.
	- 4. A 10−3–10−<sup>4</sup> Large Scale Anisotropy (LSA) amplitude is found at TeV energies [76].
	- 5. A 10−<sup>4</sup> Medium Scale Anisotropy (MSA) amplitude is observed at TeV energies [76].
	- 1. The all-particle energy spectrum measured by different experiments are in good agreement within the systematics and with the measurements of UHE experiments.
	- 2. A concave region is observed above 2 × 1016 eV with a steepening at ∼1017 eV.
	- 3. The dipole component of the LSA is smaller than 10<sup>−</sup>2.

The observed features in the all-particle energy spectrum seem to be consistent with the bending of different components in a rigidity-based scenario. However, rigidity models can be


The first PeVatrons observed in the northern hemisphere by the LHAASO experiment show that SNRs are probably not the main sources of PeV CRs in our galaxy. The observation of sources emitting photons above the PeV in the North suggests the need of a wide field of view instrument in the Southern Hemisphere to monitor the Inner Galaxy and the Galactic Center looking for super-PeVatrons (SWGO [77], STACEX [78]) to operate with CTA-South [79].

In the coming years, the LHAASO experiment is expected to be able to measure the energy spectra of different mass groups up to 1017 eV and to determine the energy of the proton knee, thus clarifying the origin of the knee in the all-particle spectrum. The suite of independent instruments that will be operated will also allow a deep study of the characteristics of the hadronic interaction models and to investigate the uncertainties related to the main techniques used to recover the elemental composition.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest

#### **References**


#### MDPI

St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Applied Sciences* Editorial Office E-mail: applsci@mdpi.com www.mdpi.com/journal/applsci

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com

ISBN 978-3-0365-3858-7