Ultra-high energy cosmic rays : 40 years retrospective of continuous observations at the Yakutsk array : Part 2 . Mass composition of cosmic rays at ultra high energies

In the paper, we describe methods for the analysis and present results for the mass composition of cosmic rays, obtained by using these techniques over a large time span. The data were obtained at the Small Cherenkov array over a 20 – year period of continuous observation and 40 – years of observations at the main Yakutsk array. Our experimental data indicate a change in the mass composition in the energy range 1016–1018 eV and is confirmed by independent results obtained by other


Introduction
In this review, we used the work of the Yakutsk group in which was introduced an estimation of the mass composition of cosmic rays at ultrahigh energies [1][2][3][4][5].In these studies we analyzed the characteristics of longitudinal and lateral development of EAS, reconstructed according to observations at the Yakutsk array.This primarily refers to measurements of Cherenkov light of EAS and muons with a threshold energy ≥ 1 GeV.These components, in accordance with the calculations, were considered the most sensitive characteristics of the shower to the atomic weight of the primary particle.The results were obtained using different models of hadron interactions [6][7][8].One can assume that all presented results on the mass composition are rather indicative, because they are not direct but indirect measurements of the mass composition and depend on many factors: 1) the accuracy of the measurement of the basic characteristics (they are all different for different arrays), 2) the conditions of registration, the techniques used, mathematical processing and selection of events, and finally 3) from a large uncertainty in the choice of a single model of hadron interactions to describe the development of extensive air showers at ultrahigh and especially in the region of huge energies.
However, it is necessary to estimate the CR mass composition in the region of ultra-high energies and compare these results with direct measurements, the so called normal composition, which were obtained at high energies from satellite and balloon measurements.This leads to a refinement of our knowledge of the nature of a e-mail: knurenko@ikfia.sbras.rub e-mail: igor.petrov@ikfia.sbras.rucosmic radiation and a better understanding of the physics of EAS development in the field of the highest energies.

The method of joint analysis of the average characteristics of the longitudinal development of EAS and their fluctuations: X max , σ (X max ), dE/dX max
In paper [9] it was suggested that composition of primary particles consists of a mixture of protons and iron nuclei.The analysis also used the superposition hypothesis, where it was assumed that the collapse of the primary nucleus did not occur at the top of the atmosphere, but at a depth corresponding to a run for the collision of nuclei in the air.Therefore, the average depth of the shower maximum for this superposition is modified by the value of the path and is lower than X max calculated from the diffusion equations of the nuclear cascade process.In the method the hydrodynamic model with n ch ∼ E 1/3 was used [6].
Let the distribution of the maximum depth of showers from primary nuclei have an exponential form with the first moment equal to their run for the nuclear interaction.Then, the average depth of the shower maximum for the sum of two exponential functions will be where η -the fraction of protons in the primary cosmic radiation, X p and X Fe -the depth of maximum development of the primary proton and iron nuclei by the chosen model of the EAS.For the dispersion we have the X max expression

EPJ Web of Conferences
Here λ p λ Fe are respectively the path for the nuclear interaction of the proton and iron nuclei.The β multiplier takes into account the increase of the dispersion X max due to fluctuations of the inelasticity coefficient and is taken to be (1 −1 < k >) −1 , where ¡k¿ is the mean value of the inelasticity coefficient of the leading particle.If we assume that the ratio λ Fe (E)/λ p (E) is constant and known, then Eqs. ( 1) and ( 2), in the framework of a two-component composition, can determine the proportion of protons in the primary radiation and the cross section of the proton -nucleus of an atom from the experimental values X max and D(X max ).Technically, the above simplification can be extended to the multi-component composition of the primary particles, but the accuracy of formulae of type 1 and 2 for X max and D (X max ) will be slightly worse.
We also found an indication for a gradual increase of the protons in the energy region 3 • 10 17 -3 • 10 18 eV.

Method of comparing the asymmetry of distribution at different X max for different fixed energies
This method does not depend upon models of air shower development.Since the distribution of X max at fixed energy is formed by nuclei of different masses its shape will therefore reflect their contribution to a statistical contribution of X max , g/cm 2 .This is understandable, since showers produced by particles of different masses have either rapid development, for example the iron nucleus, or slow development as happens if we consider protons.In paper [2] the ideology of asymmetry of X max distributions at different energies was used.Its essence is as follows.The value of the effective cross section for inelastic collisions of protons with air nuclei on the distribution of heights of the maxima of the EAS in the energy range of the primary particles 10 17 -10 18 eV and 10 18 -10 19 eV was studied.Right-hand side of such distributions is determined mainly by the effective cross section for inelastic collisions of protons with air nuclei.The height of the maximum development of the shower was determined by the spatial distribution of the Cherenkov light at distance ranges of 100-600 m from the shower axis.
Next, let us assume that at high energies only protons are present and ration the distribution at lower energies by the proton (deeper than 700 g/cm 2 ), by simply subtracting the estimated fraction of nuclei in the primary radiation with energy ∼ 10 17 -10 18 eV.
Thus, in the framework of QGSJET 01 an indication was found that the mass composition of the PCR in the transition from the energy range (5-30)• 10 17 eV to (3-10)•10 18 eV changes.For E 0 • 3•10 18 eV primary cosmic radiation consists of ∼70% of protons and helium nuclei, the proportion of other nuclei does not exceed ∼30% (see Fig. 1).

Multicomponent analysis
In order to interpret the experimental data of the Yakutsk array we used CORSIKA code (v.6.0.QGSJET model) to generate a database of X max and ρ e (600).Simulations were made for five primaries (P, He, C, Si, Fe) and three energies 10 17 , 10 18 , 10 19 eV.For each energy we simulated 100 showers in the standard atmosphere.In this paper, we used two-dimensional probability density F(X max , ρ e (600)), preliminary standardized experimental data of the entire array data (X max , ρ e (600)) for a given energy.At numerical implementation of this method instead of (X max , ρ e (600)) and used variables τ and ρ e : where σ is the standard deviation of the value.Standardization performed on pooled data (X max , ρ e (600)) for all groups of nuclei and each energy 10 17 , 10 18 , 10 19 eV.Distribution on X max and ρ e (600) separately and joint distributions for τ and ρ are described respectively by dimensional F (X max ) F (ρ(600)) and twodimensional f(τ , ρ) logarithmically normal distribution.For each given energy and different types of primary nuclei, including for nuclei, combined in groups P + He, C, Si + Fe, probability distribution density f(τ , ρ) were plotted (Fig. 2).The intersection of f(τ , ρ) layers gives lines m1 and m2, which optimally separates nuclei into 3 groups: (P + He), C and (Si + Fe) respectively.
Figure 2 shows the result of a multi-component analysis of the data binding (X max ρ e (600)) of the Yakutsk array.A cloud of points in such a representation reflects standardized values, and lines represent areas that are directly associated with the mass number of the primary particle.In this case, the line m1 and lines m2 optimally separate nuclei into groups (P + He), C and (Si + Fe).
Analysis has shown that the proportion of nuclei (P + He) increases from 50% to 53%, and the proportion of carbon nuclei from 23% to 31%.At the same time the proportion of nuclei of heavy chemical elements decreases from 27% to 16% when increasing energy from 2.4•10 17 to 4.8•10 18 eV.

Proportion of muons analysis method depending on the length of the track of the particles in the atmosphere
In paper [5] we considered the dependence ρ µ /ρ s the length of the track of the particles after the maximum of EAS λ = X 0 / cos θ -X max where X 0 = 1020 g/cm 2 for Yakutsk.Here X max was determined from measurements of the Cherenkov light and ρ µ and ρ s by measuring large EAS.Next, the experiment was compared with model calculations of QGSJETII 03 and EPOS (see Fig. 3).It is known that X max of showers greatly differs depending on the primary nucleus, therefore, this fact can be used to analyze the mass composition of cosmic rays, for example, by fixing the parameter λ and analyzing fluctuations in the ratio ρ µ /ρ s .This method is somewhat similar to that proposed by Christiansen in 1981 [10].With sufficient precision of measurements of each parameter (better than 5%) in the distribution single peaks from different nuclei can be allocated.
Comparison of the distribution of the muon proportion with calculation indicates a mixed composition at energies above 10 18 eV.Large fluctuations do not allow allocating separate groups of nuclei with a good precision nor evaluating the percentage of each group.But, the use of a pure response of muon detectors leads to the conclusion that the MC at energies 10 18 -10 19 eV is light [11].

Evaluation of the mass composition at average depth of maximum of EAS development. Interpolation method
In papers [9,12], the dependence of X max on energy in the range ∼ 10 15 to 5 • 10 19 eV was considered.The MC of PCR was evaluated by this formula: where a i is the relative proportion of nuclei with mass number A i .
In each case experimental data was compared with QGSJET 03 calculations for proton and iron in the frame of the superposition model: ln A = ((P ex p. − P p ))/(P Fe − P p )) • ln A Fe (8) where P i -the parameter that characterizes the longitudinal development of air showers X max .
Figure 4 shows the dependence of X max on energy (dots) derived from experiment and from simulation (lines) of this characteristic calculated using the QGSJET 03 and SIBYLL model for proton and iron nuclei.Figure 5 shows Yakutsk array results of MC PCR derived by the method described above.The data was obtained in the frame of QGSJETII 03 model and dual component MC (proton iron).The value < ln(A) > was determined by the interpolation method in each case.

Conclusion
a) For more than 40 years the Yakutsk array has continuously recorded ultrahigh energy air showers.We obtain information about all the main components of  the shower: electrons, photons, hadrons and muons.All these data from different times were used to estimate the mass composition of cosmic rays using different methods.This follows from numerous publications in journals and proceedings of scientific conferences.According to the data shown in Fig. 6 [13] the mass composition is not uniform over a wide energy range, and has a peak at (0.8-2) 10 17 eV.
Figure 7 shows the latest results on the MC obtained using the interpolation method (see Sect. b) At Yakutsk we measured the energy spectrum of CR [14,15] and evaluated the MC over a wide energy range [16].If we compare the energy scale with the studied CR spectrum and obtained results of the cosmic ray MC, than we observe a matching of energy intervals, where the change in the shape of the spectrum and the change in the value of A =< ln A > are the same.Most likely, these two results are related and caused by the same astrophysical processes.
c) For the boundary of the transition from galactic to metagalactic cosmic rays recently developed nonlinear kinetic theories of CR acceleration in supernova remnants has allowed not only to achieve agreement in the shape of the CR spectrum up to energies ∼ 10 17 with experimental data, but also to choose a class of SNR, which is responsible for the MC of particles similar to those observed in satellite, balloon and ground experiments [9].This is confirmed by the results of the calculation of work [16], which are shown in Fig. 7 (lines).Figure 7 shows a comparison of MC obtained at different arrays, with MC generated in the sources, which are supernova remnants.There is not only a satisfactory agreement of experimental data with calculation in the energy range 10 15 -10 19 eV, but also an indication that the sharp change in the MC at ∼ 2 • 10 17 eV may be associated with the boundary of the transition from galactic CR to metagalactic CR.In this case the mass composition of CRs at energies above ∼ 2 • 10 17 eV should be represented primarily by protons, which is consistent with the mass composition obtained by the Yakutsk EAS array.

Figure 3 .
Figure 3.A of the ρ µ /ρ s relation normalized to the track length 500 g/cm 2 .On the left -according to models QGSJet II(FLUKA), on the right -according to EPOS(UrQMD).

Figure 4 .
Figure 4. Dependence of X max on energy.Lines are calculated values for proton and iron nuclei.

Figure 5 .
Figure 5.The mass composition of the highest energy Cosmic rays is obtained at Yakutsk.Model QGSJETII-03.

Figure 5
Figure 5 shows that the nature of the dependence of < ln A > with increasing energy changes reaching a maximum in the energy range (5-30)•10 16 eV.This means, that the MC of PCR changes after the first kink in the spectrum at ∼ 3 • 10 15 eV, requiring heavier particles at (3-30)•10 16 eV and then, starting at 3 • 10 17 eV, becomes much lighter.

Figure 6 .
Figure 6.Estimation of MC CR with different methods using different characteristics of air showers.

Figure 7 .
Figure 7. Mass composition of CRs from measurements of different EAS arrays.Lines show MC calculated from paper[16] in the case of near SNR.
2.6) and the QGSJET 04 model.The figure also shows results obtained at other compact and large arrays.From Fig. 7 it follows 04003-p.4ISVHECRI 2014 that the MC is undergoing change in the energy range (8-20) •10 16 eV and (8-20) •10 18 eV.Most likely this is due to the nature of the formation of cosmic rays in the sources and their distribution in galactic and intergalactic space.