In this article a new method is proposed for estimating the mass composition of cosmic rays in individual events with energies above inlineformula idIEq1texmath idIEq1TeX1.25 times 1019texmathinlineformula eV. It is based on a joint analysis of experimental data and simulation results obtained using the QGSjetII.04 model for muons with threshold energy inlineformula idIEq2texmath idIEq2TeXEmu   1.0 times cos theta texmathinlineformula GeV in air showers with zenith angles up to 60. The data from groundbased and underground scintillation detectors of the Yakutsk EAS array were used. Separate groups of nuclei and other primary particles were found.
['Astrophysics - High Energy Astrophysical Phenomena']
Estimation of the Composition of UltraHigh Energy Cosmic Rays Using the Muon Correlation Method Based on Yakutsk EAS Array Data
{'Estimation of the composition of ultra-high energy cosmic rays using the muon correlation method based on Yakutsk EAS array data': 'A.V.Glushkov * 1) , L. T. Ksenofontov * , K. G. Lebedev * , A. V. Saburov  \n* Yu. G. Shafer Institute of Cosmophysical Research and Aeronomy SB RAS 677027 Yakutsk, Russia  \nIn this article a new method is proposed for estimating the mass composition of cosmic rays in individual events with energies above 1 . 25 × 10 19 eV. It is based on a joint analysis of experimental data and simulation results obtained using the qgsj etii .04 model for muons with threshold energy 𝐸 𝜇 = 1 . 0 × cos 𝜃 GeV in air showers with zenith angles up to 60 degrees. The data from ground-based and underground scintillation detectors of the Yakutsk EAS array were used. Separate groups of nuclei and other primary particles were found.', '1. INTRODUCTION': 'Cosmic ray (CR) mass composition in ultra-high energy range (above ∼ 10 17 eV) so far is studied only in general terms. Information about it is often contradictory. Atmospheric depth 𝑥 max /emdash.cyr at which the maximum number of particles ( 𝑁 max ) is achieved during the extensive air shower (EAS) development /emdash.cyr is a standard parameter for obtaining information about the CR composition. Simulations demonstrate that depending on the mass number 𝐴 of primary particles, the process of air shower longitudinal development, taking into account fluctuations, leads to certain values of the mean depth ⟨︀ 𝑥 𝐴 max ⟩︀ and its dispersion 𝛿 ( 𝑥 𝐴 max ) . The ⟨ 𝑥 max ⟩ value is connected to conventional mean atomic number ⟨ 𝐴 ⟩ of primary particles via a simple ratio, which follows from the nucleon superposition principle [1]:  \n⟨ ln 𝐴 ⟩ = ⟨︀ 𝑥 𝐴 max ⟩︀ -⟨ 𝑥 𝑝 max ⟩ ⟨ 𝑥 Fe max ⟩ - ⟨ 𝑥 𝑝 max ⟩ × ln 𝐴 Fe , (1)  \nwhere values for primary protons ( 𝑝 ) and iron nuclei (Fe) were obtained using one or another hadron interaction model. In experiment the ⟨ ln 𝐴 ⟩ value is determined by substituting the ⟨︀ 𝑥 𝐴 max ⟩︀ with measured value ⟨ 𝑥 exp max ⟩ . This technique is successfully utilized at Auger [2] and Telescope Array (TA) [3] experiments for individual EAS events with registered fluorescent light emission.  \nSimulations have revealed that expression (1) is applicable to other EAS parameters that are sensitive the 𝑥 max value. This is especially true for muon component, which is actively studied in many experiments in the wide range of primary energy. In this case one has to deal with the following relations:  \n⟨ ln 𝐴 ⟩ exp = ⟨ 𝑧 exp 𝐴 ⟩ · ln 𝐴 Fe , (2)  \n⟨ 𝑧 exp 𝐴 ⟩ = ln ⟨︀ 𝜌 exp . 𝜇 ⟩︀ -ln ⟨︀ 𝜌 𝑝 𝜇 ⟩︀ ln ⟨︀ 𝜌 Fe 𝜇 ⟩︀ -ln ⟨ 𝜌 𝑝 𝜇 ⟩ , (3)  \nwhere ⟨︀ 𝜌 exp . 𝜇 ⟩︀ is the muon density registered in experiment and ⟨︀ 𝜌 𝑝 𝜇 ⟩︀ and ⟨︀ 𝜌 Fe 𝜇 ⟩︀ are densities obtained for primary protons ( 𝑝 ) and iron nuclei (Fe) using full simulation of the measurement process with a real detector. The value (3) has been heavily discussed recently due to existing fundamental disagreements [4-7].  \nAt the Yakutsk array the muon component of EAS has been registered since the very beginning of its operation in 1974. To date, a significant experimental material has been accumulated, which is episodically analyzed as the notion of the nature of CR evolves. Below, a new method is considered for estimating the CR mass composition based on these data. It relies on formulas (2) and (3) in relation to the physical picture of EAS development.', '2..1 General formulation of the problem': "The Yakutsk EAS array stands out from other similar instruments by its complex design: it simultaneously measures charged and electromagnetic particles with 2 × 2 -m 2 surface scintillation detectors (SD), muons with energy above 1 × sec 𝜃 GeV with similar ground shielded detectors (MD) with an area 20 -36 m 2 and EAS Cherenkov light emission (CLE). In the work [8] lateral distribution functions (LDFs) of SD and MD responses and zenith-angular dependencies of the corresponding densities at axis distance 𝑟 = 600 m  \nwere studied in showers with 𝐸 0 = 10 19 eV and zenith angles 𝜃 ≤ 60 ∘ . Experimentally measured values were compared to estimations obtained within frameworks of hadron interaction models qgsj et01 [9] and qgsj etii .04 [10] using the corsika code [11]. The details of the SD and MD response calculation are given in [7,8]. The whole considered data set indicates a certain agreement between experiment and theory. In [8] it is stated that the probable CR composition is close to protons with possible fraction of primary photons about 9 %. This conclusion was based on the analysis of mean LDFs in separate groups of showers with zenith-angular directions ⟨ cos 𝜃 ⟩ = 0 . 95 , 0 . 90 , 0 . 85 , 0 . 80 , 0 . 75 , 0 . 65 and 0 . 55 . Here we present the results of a further study of the CR mass composition with energy above 1 . 25 × 10 19 eV in individual events. Let's reiterate on some major points of the procedure for processing of registered events adopted at the Yakutsk array, which we will need later.", '2..2 Energy estimation': 'CLE contains information about approximately 80 % of primary energy 𝐸 dispersed by a shower in the atmosphere and provides the possibility to estimate the 𝐸 value with the use of calorimetric method [1216]. The EAS energy was determined from following relations [12]:  \n𝐸 = 𝐸 1 × 𝑆 600 (0 ∘ ) 𝐵 [eV], (4) 𝑆 600 (0 ∘ ) = 𝑆 600 ( 𝜃 ) × × exp [︂ (sec 𝜃 -1) × 1020 𝜆 ]︂ [m -2 ], (5) 𝜆 = 400 ± 45 [g/cm 2 ], (6)  \nwhere 𝐸 1 = (4 . 1 ± 1 . 4) × 10 17 eV, 𝐵 = 0 . 97 ± 0 . 04 . Numerically the 𝐸 1 proportional coefficient is equal to the energy of a vertical shower with density of the SD response measured at axis distance 600 m 𝑆 600 (0 ∘ ) = 1 m -2 . The experiment essentially measures primary energy in units of some reference EAS event. Later the parameters of the relation (4) changed slightly: 𝐸 1 = (4 . 8 ± 1 . 6) × 10 17 eV, 𝐵 = 1 . 00 ± 0 . 02 and the expression for attenuation length (6) took the following form [13,14]:  \n𝜆 =(450 ± 44) + (32 ± 15) × × log 10 𝑆 600 (0 ∘ ) [g/cm 2 ]. (7)  \nIn [15, 16] we have reconsidered the energy calibration according to results of simulations performed  \nwith corsika code. After that, the relations (4), (5) and (7) were finally adopted with parameters 𝐸 1 = (3 . 76 ± 0 . 30) × 10 17 eV and 𝐵 = 1 . 02 ± 0 . 02 [16]. The error in formula (4) was mainly conditioned by the precision of absolute calibration of CLE detectors and incorrect estimation of the atmospheric transparency [12].', '2..3 Primary events reconstruction': "Arrival direction of a shower is reconstructed from relative delays of SD's firing using the approximation of flat shower front. The axis is located with the use of the Greisen-Linsley LDF approximation with parameters obtained at the Yakutsk array during the initial period of its operation [17]:  \n𝑆 ( 𝑟, 𝜃 ) = 𝑆 600 ( 𝜃 ) · 600 𝑟 · (︂ 600 + 𝑟 M 𝑟 + 𝑟 M )︂ 𝛽 ( 𝜃 ) -1 , (8)  \nwhere 𝑆 600 ( 𝜃 ) is the all-particle response density measured by scintillation SDs at axis distance 𝑟 = 600 m and 𝑟 M is the Moliere radius. It depends on atmospheric parameters 𝑇 [K] and 𝑃 [mb] [17,18]:  \n𝑟 M = 7 . 5 × 10 4 𝑃 × 𝑇 273 . (9)  \nValues of 𝑇 and 𝑃 are constantly monitored and recorded in the primary data bank at every triggering of the array. The annual mean value of the Moliere radius for Yakutsk is ⟨ 𝑟 M ⟩ ≈ 70 m. The structural parameter 𝛽 was determined in [17]:  \n𝛽 ( 𝜃 ) = 1 . 38 + 2 . 16 × cos 𝜃 + +0 . 15 × log 10 𝑆 600 ( 𝜃 ) . (10)  \nLater it was established that the LDF approximation (8) poorly describes experimental data for showers with 𝐸 ≥ 10 19 eV in the axis distance range 𝑟 ≃ 30 -2000 m, and a modified LDF approximation was introduced [19]:  \n𝑆 ( 𝑟, 𝜃 ) = 𝑆 600 ( 𝜃 ) · (︂ 600 𝑟 )︂ 2 · (︂ 600 + 𝑟 M 𝑟 + 𝑟 M )︂ 𝛽 ( 𝜃 ) -2 × × (︂ 600 + 𝑟 1 𝑟 + 𝑟 1 )︂ 10 , (11)  \nwhere 𝑟 1 = 10 4 m. The 𝛽 ( 𝜃 ) parameter in this case virtually does not change with the energy but depends on the zenith angle. The arrival direction, axis coordinates and the 𝑆 600 ( 𝜃 ) value are determined  \nvia 𝜒 2 -minimizations during the preliminary processing of registered events after a diurnal cycle of the array operation.  \nThe LDF of muon component in this work is approximated in the axis distance range 𝑟 ≈ 300 -2000 m with the function:  \n𝑆 𝜇 ( 𝑟, 𝜃 ) = 𝑆 𝜇, 600 ( 𝜃 ) · (︂ 600 𝑟 )︂ 0 . 75 × × (︂ 600 + 𝑟 0 𝑟 + 𝑟 0 )︂ 𝑏 𝜇 ( 𝜃 ) -0 . 75 · (︂ 600 + 𝑟 1 𝑟 + 𝑟 1 )︂ 9 , (12)  \nwhere 𝑟 0 = 280 m and 𝑟 1 = 2000 m. Absolute and daily calibrations of SDs and MDs were described in [7].", '2..4 The muon correlation method': "Here we describe a method for estimating the CR composition in individual EAS events from muon component, which was called 'the muon correlation method' . It consists in comparing the measured density of the MD response 𝑆 exp 𝜇 ( 𝐸,𝜃, 𝑟 ) in a shower with energy 𝐸 and zenith angle 𝜃 at axis distance 𝑟 with the expected value 𝑆 sim 𝜇 ( 𝐸 * , 𝜃, 𝑟 ) obtained within the framework of the qgsj etii .04 model for primary proton with fixed energy 𝐸 * , the same zenith angle at the same axis distance. Simulations indicate that for registered events with energies 𝐸 ≥ 10 19 . 1 eV and 𝜃 ≤ 60 ∘ , which were selected for further analysis, it is possible to determine the expected muon density from the relation:  \n𝑆 sim 𝜇 ( 𝐸,𝜃, 𝑟 ) = 𝑆 sim 𝜇 ( 𝐸 * , 𝜃, 𝑟 ) · 𝑆 exp 600 ( 𝐸,𝜃 ) 𝑆 sim 600 ( 𝐸 * , 𝜃 ) , (13)  \nwhere 𝐸 * = 3 . 16 × 10 19 eV. The 𝑆 sim 600 ( 𝐸 * , 𝜃 ) value reflects the zenith-angular dependency of the SD response density at 𝑟 = 600 m. For a vertical shower it is connected via expression (4) to the energy 𝐸 * , which is included in the simulation, and within ≈ 6 % agrees with primary energy estimation obtained at the Yakutsk array [16]. The 𝑆 exp 600 ( 𝐸,𝜃 ) = 𝑆 600 ( 𝜃 ) parameter is estimated experimentally and is included in expression (4). Further we consider the correlation between densities 𝑆 sim 𝜇 ( 𝐸,𝜃, 𝑟 ) at 𝑟 = 600 and 1000 m and the corresponding experimentally derived values.  \nThe sequence of calculations for determining the 𝑆 sim 𝜇 ( 𝐸,𝜃, 𝑟 ) value with relation (13) is demonstrated in Fig. 1 for a particular registered shower with energy 𝐸 = 10 19 . 73 eV, which was estimated according to formulas (4), (5), (7), and zenith angle 𝜃 = 48 ∘ (sec 𝜃 = 1 . 49 ± 0 . 03) . For this event the value  \nFig. 1. Zenith-angular dependencies of responses in SD (empty circles) and MD with 1 . 0 × sec 𝜃 GeV threshold at 600 m (empty squares, crosses) and 1000 m (dark squares, x-es) from the axis in showers with energy 𝐸 = 10 19 . 5 eV obtained within the framework of the qgsj etii .04 model for primary protons ( 𝑝 ) anf iron nuclei (Fe). Simulation results obtained using the eposlhc model (dotted lines) were obtained for primary protons.  \n<!-- image -->  \nlog 10 [ 𝑆 exp 600 ( 𝐸, 48 ∘ )] = 1 . 77 ± 0 . 02 was obtained (dark circle '1'). The estimated energy of this event isn't used further and is only given for information. The value of parameter log 10 [ 𝑆 sim 600 ( 𝐸 * , 48 ∘ )] = 1 . 54 for relation (13) follows from the zenith-angular dependency displayed with empty circles. The length of the arrow '1' equals to the normalization shift:  \nΔlog 10 [ 𝑆 600 (48 ∘ )] = = log 10 [ 𝑆 exp 600 ( 𝐸, 48 ∘ )] -log 10 [ 𝑆 sim 600 ( 𝐸 * , 48 ∘ )] = = 1 . 77 -1 . 54 = 0 . 23 ,  \nwhich is to be added to the calculated value log 10 [ 𝑆 sim 𝜇, 600 ( 𝐸 * , 48 ∘ )] = 1 . 08 (curve with squares at 𝑟 = 600 m) in order to obtain the expected density of the MD response:  \nlog 10 [ 𝑆 sim 𝜇, 600 ( 𝐸, 48 ∘ )] = = 1 . 08 + 0 . 23 = 1 . 31 ± 0 . 02 . (14)  \nIt appears that, within measurement errors, the value (14) is close to the experimentally obtained density (dark square '1'):  \nlog 10 [ 𝑆 exp 𝜇, 600 ( 𝐸, 48 ∘ )] = 1 . 30 ± 0 . 04 . (15)  \nNote that this result does not depend on the chosen hadron interaction model. To make it sure, let's consider, for example, zenith-angular characteristics of SD and MD responses obtained with the above mentioned technique using the epos-lhc model [20]. It is evident from Fig. 1 that both models give curves with identical shapes and differing in absolute values by factor 1.1. Let's determine the parameter log 10 [ 𝑆 sim 𝜇, 600 ( 𝐸 * , 48 ∘ )] = 1 . 60 for relation (13). The length of arrow '1' in this case is  \nΔlog 10 [ 𝑆 600 (48 ∘ )] = = log 10 [ 𝑆 exp 600 ( 𝐸, 48 ∘ )] -log 10 [ 𝑆 sim 600 ( 𝐸 * , 48 ∘ )] = = 1 . 77 -1 . 59 = 0 . 18 ,  \nwhich is to be added to the new calculated value log 10 [ 𝑆 sim 𝜇, 600 ( 𝐸, 48 ∘ )] = 1 . 12 (dotted line at 𝑟 = 600 m) in order to obtain the expected MD response density:  \nlog 10 [ 𝑆 sim 𝜇, 600 ( 𝐸, 48 ∘ )] = = 1 . 12 + 0 . 18 = 1 . 30 ± 0 . 02 . (16)  \nIt coincides with the value (15). From the considered example of two different hadron interaction models, with a certain degree of confidence one can assume that the considered shower was initiated by primary proton.  \nThe equality of values (14) and (16) is not a random coincidence, but a natural characteristic of EAS development. In [21] it was shown that estimations of the CR mass compositions obtained from muon fraction with the use of expressions (2) and (3) within frameworks of different models were virtually identical. This feature of the behaviour of muon fraction is utilized in expression (13), where it is clearly seen in the identically transformed expression:  \n𝑆 sim 𝜇 ( 𝐸,𝜃, 𝑟 ) = 𝑆 exp 600 ( 𝐸,𝜃 ) · 𝑆 sim 𝜇 ( 𝐸 * , 𝜃, 𝑟 ) /𝜀 𝑆 sim 600 ( 𝐸 * , 𝜃 ) /𝜀 (17)  \nwith 𝜀 ≈ 1 . 1 for epos-lhc . This coefficient varies within (8 -10) % in different hadron interaction models. From a physical standpoint, a formal substitution 𝜀 → 𝐸 * means normalization of all responses in (17) by primary energy.", '2..5 The giant event': "On May 7 1989 a unique event was registered at the Yakutsk array, with energy 𝐸 ≈ 1 . 12 × 10 20 eV and zenith angle 𝜃 = 58 . 7 ∘ ± 1 . 3 ∘ , and to this day it remains the most powerful event in its data bank.  \nFig. 2. The unique EAS event with energy 𝐸 ≈ 1 . 12 × 10 20 eV on the array map: digits /emdash.cyr summary SD responses; empty circles /emdash.cyr stations that were not operating during the observation; rectangles - MDs with 1 . 0 × sec 𝜃 GeV threshold; triangles /emdash.cyr CLE detectors; cross and arrow indicate shower axis and arrival direction.  \n<!-- image -->  \nFig. 3. Lateral distribution of particles in the Arian event registered in experiment by SDs (circles, triangles) and by MDs with 1 . 0 × sec 𝜃 GeV threshold (dark squares). Solid and dotted lines represent theoretical LDFs of SD and MD responses calculated within the framework of the qgsj etii .04 model for primary protons with energy 10 20 eV and 𝜃 = 60 ∘ . Open squares and corresponding curve represent approximation (12) of experimental data.  \n<!-- image -->  \n10  \nThe shower was dubbed Arian . In Fig. 2 the shower footprint is shown on the array map with its axis, which was located with Δ 𝑟 = 26 m accuracy. The fired detector stations are shown with filled circles, their size is proportional to the logarithm of the number of measured SD responses. This shower covered the whole array and provided the opportunity to test its  \noperation. In Fig. 3 readings of all fired SD and MD are shown. Empty circles indicate measurements with six additional stations concentrated around the array center with 250 m spacing. Each station contained a single 2-m 2 scintillation detector. Also during this period, a compact cluster of 0 . 25 -m 2 scintillation detectors spaced by 50 m operated in the central area of the array. Readings of these detectors are shown with black triangles. All the data provided a precise measurement of the value log 10 [ 𝑆 exp 600 ( 𝐸, 59 ∘ )] = 1 . 735 ± 0 . 013 (dark circle '2' in Fig. 1). Solid curve indicates best fit of the SD data with approximation (11). It happened to be very close, in terms of both absolute value and shape, to the LDF obtained using the qgsj etii .04 model for primary protons with energy 10 20 eV and 𝜃 = 60 ∘ . Empty squares represent approximation (12) of the experimentally measured lateral distribution of muons with energies 𝐸 𝜇 ≥ 1 . 92 GeV. According to the 𝜒 2 -test it satisfies the experimental data with parameters 𝑏 𝜇 ≈ 0 . 28 and  \nlog 10 [ 𝑆 exp 𝜇, 600 ( 𝐸, 59 ∘ )] = 1 . 81 ± 0 . 04 . (18)  \nThe density (18) is represented in Fig. 1 with dark square '2'. Let's determine the normalization shift for this shower according to the procedure described in Section 2.4:  \nΔlog 10 [ 𝑆 600 (59 ∘ )] = = log 10 [ 𝑆 exp 600 ( 𝐸, 59 ∘ )] -log 10 [ 𝑆 sim 600 ( 𝐸 * , 59 ∘ )] = = 1 . 735 -1 . 225 = 0 . 5 , (19)  \nwhich is shown with arrow '2' in Fig. 1 and calculate the expected density of muon response:  \nlog 10 [ 𝑆 sim 𝜇, 600 ( 𝐸, 59 ∘ , 𝑝 )] = = 1 . 04 + 0 . 5 = 1 . 54 ± 0 . 02 , (20)  \nwhich turned out to be not only lower than the expected value from primary protons by Δ 𝜇 = 1 . 81 -1 . 54 = 0 . 27 , but also lower than the expected value from iron nuclei:  \nlog 10 [ 𝑆 sim 𝜇, 600 ( 𝐸, 59 ∘ , Fe )] = = 1 . 16 + 0 . 5 = 1 . 66 ± 0 . 02 (21)  \nby Δ 𝜇 = 1 . 81 -1 . 66 = 0 . 15 . Similar result follows from muon density at axis distance 𝑟 = 1000 m, where the expected value for iron nuclei is:  \nlog 10 [ 𝑆 sim 𝜇, 1000 ( 𝐸, 59 ∘ , Fe )] = = 0 . 51 ± 0 . 5 = 1 . 01 ± 0 . 02 , (22)  \nand the measured value is  \nlog 10 [ 𝑆 exp 𝜇, 1000 ( 𝐸, 59 ∘ )] = 1 . 15 ± 0 . 04 , (23)  \nshown in Fig. 1 with dark square.  \nThe relation (13) is convenient for implementing this method because it eliminates the need to perform heavy and time consuming Monte Carlo calculations of 𝑆 sim 𝜇 ( 𝐸,𝜃, 𝑟 ) in each individual shower. In this case we have performed them only for showers with 𝐸 * = 10 19 . 5 eV and zenith angle values separated by step Δcos 𝜃 = 0 . 1 . This approach is fundamentally different from that used at the TA and Auger arrays, where each registered EAS event is simulated with Monte Carlo method ∼ 10 6 times to estimate its main parameters. The previously performed computational analysis of the EAS development have demonstrated that estimations of the CR mass composition obtained with the use of the relation (13) are quite correct. This is evident from the above-mentioned example of the giant shower where experimentally measured parameters were 𝜃 = 58 . 7 ∘ ± 1 . 3 ∘ and log 10 𝑆 exp 600 ( 𝐸, 59 ∘ ) = 1 . 735 ± 0 . 013 . Densities of the MD response calculated for this event with Monte-Carlo method using the qgsj etii .04 model for primary proton are  \nlog 10 [ 𝑆 sim 𝜇, 600 ( 𝐸, 59 ∘ , 𝑝 )] = 1 . 65 ± 0 . 02 , (24)  \nand  \nlog 10 [ 𝑆 sim 𝜇, 1000 ( 𝐸, 59 ∘ , 𝑝 )] = 0 . 89 ± 0 . 02 , (25)  \nwhich within errors are close to the expected values from the relation (13) and Fig. 1.", '3..1 Correlation dependencies': "The geometry of the Yakutsk array during the observational period 1974.01.01 - 1990.06.23 is shown in Fig. 2. In 1976 two new 36-m 2 MDs were deployed at ≈ 500 m from the array center. In 1986 another three 20-m 2 MDs were deployed. During the renovation in 1990-1992, the outermost SD stations were dismantled and moved to a central circle with a radius 2 km. The current analysis of the MD data includes only showers with axes lying inside this circle.  \nIn Fig. 4 a sample of 127 showers with 𝐸 ≥ 10 19 . 1 eV and 𝜃 ≤ 60 ∘ is shown. It includes events with muons registered at axis distances 600 m (dark circles) and 1000 m (open circles). Showers with at least two fired MDs were selected for the analysis. A detector was considered fired even if it didn't register a single particle  \nFig. 4. Correlations of the response densities recorded by MDs with 1 × sec 𝜃 GeV threshold at axis distances 600 m (dark circles) and 1000 m (open circles) in 127 showers with energies 𝐸 ≥ 10 19 . 1 eV and 𝜃 ≤ 60 ∘ . Horizontal axis /emdash.cyr densities determined in experiment, vertical axis /emdash.cyr those calculated from simulation results performed within the framework of qgsj etii .04 model for primary protons (solid line), iron nuclei (dashed line). Dotted line highlights EASs with abnormally high muon content.  \n<!-- image -->  \nbut was in standard accepting mode. Lines represent the expected correlations between measured and expected muon densities in EASs initiated by different primary particles. Solid line corresponds to the correlation for primary protons. Dashed line indicates correlation of muon densities in EASs originating from iron nuclei. It was obtained by averaging the differences of muon densities shown in Fig. 1 in all interval of zenith angles. One can see that it is shifted to the right along the horizontal axis from the proton correlation line by the value  \n⟨ Δ 𝑝 -Fe ⟩ = log 10 𝑆 sim 𝜇 ( 𝐸,𝜃, 𝑟, Fe ) 𝑆 sim 𝜇 ( 𝐸,𝜃, 𝑟, 𝑝 ) ≈ 0 . 15 . (26)  \nThe dotted line in Fig. 4 indicatively highlights a few showers (further marked with the symbol 'X') with abnormally high muon content. It is shifted to the right along the horizontal axis from the proton correlation line by the value  \n⟨ Δ 𝑝 -X ⟩ = 2 ×⟨ 𝑝 -Fe ⟩ ≈ 0 . 30 . (27)  \nThe considered sample includes showers with at least 2 MDs fired within the axis distance range 800 -1600 m. Their readings were approximated with LDF (12) and  \nthe densities 𝑆 exp 𝜇, 1000 ( 𝐸,𝜃 ) were determined. In this case, the steepness parameter 𝛽 𝜇 ( 𝜃 ) was taken equal to the mean value following from the previously obtained experimental data [22], as it was impossible to calculate it in the majority of individual events. About 27% of all events had one MD with readings (or more, in rare cases) at axis distances from 300 m to 800 m. In such a case, the 𝑆 exp 𝜇, 600 ( 𝐸,𝜃 ) density was determined in addition to 𝑆 exp 𝜇, 1000 ( 𝐸,𝜃 ) . It was the most optimal approach to data selection which had its impact on the distribution of some events presented in Fig. 4.  \nFig. 5. Deviations of individual pairs of muon densities from correlation lines displayed in Fig. 4 at axis distance 1000 m.  \n<!-- image -->  \n10  \nDeviations of individual pairs of muon densities from the correlation lines presented in Fig. 4 are shown in Fig. 5:  \n𝜂 = 𝑆 exp 𝜇 /𝑆 sim 𝜇 ( 𝑝 ) (28)  \nHere zero-value on horizontal axis corresponds to the position of proton peak according to the qgsj etii .04 model. Labels 'Fe' and 'X' indicate shifts (26) and (27) relative to this peak. Events marked with the symbol 'D' refer to events with observed muon deficit. They are positioned at peaks with mean values in the region of their location at 68% confidence level. For all considered groups of events the following values were obtained for deviations (28):  \n𝜂 ( 𝑝 ) = 0 . 99 +0 . 08 -0 . 08 , (29)  \n𝜂 ( Fe ) = 1 . 44 +0 . 18 -0 . 12 , (30)  \n𝜂 ( X ) = 1 . 95 +2 . 05 -0 . 33 , (31)  \n𝜂 ( D ) = 0 . 09 +0 . 05 -0 . 02 . (32)  \nThus, the distribution of peaks presented in Fig. 5 can be interpreted as a probable composition of primary particles.", '3..2 CR composition': "The differentiated data on the CR composition depending on the EAS energy are given in Table 1. It is evident that the fraction of protons in this sample is nearly constant and amounts to ≈ 69 % except for the last row. The fraction of showers with muon deficit is also roughly constant and is about 10%. It does not contradict our previous estimation [23,24]. The fraction of iron nuclei gradually increases from 0 to 24%. In [23] this fraction was estimated to be about 36% from the sample of 33 events. The difference between the previous results and Table 1 probably arises from different sizes of the considered samples. It is also possible that different methods of estimation of the CR composition are also reflected here.  \nThere are 5 showers worth noting in the column 'X' with the observed abnormally high muon content. A detailed analysis of the obtained data have confirmed that these were the real densities recorded by MDs. From lateral distribution of particles in the most powerful event presented in Fig. 3 it is evident that the signals recorded by SDs and MDs were close in value and probably were produced by muons with energies 𝐸 𝜇 ≥ 1 . 92 GeV. At axis distances 𝑟 ≥ 400 m the contribution from electro-magnetic component is almost absent. So far it is hard to interpret this fact in any way.  \nThe positions of peaks remained virtually unchanged relative to each other. This indicates that the muon correlation method is not sensitive to the choice of hadron interaction model.  \nLet's estimate the parameter ⟨ 𝑧 exp 𝐴 ⟩ (3) from the data in Table 1. For this we will use the expression:  \nln ⟨︀ 𝜌 exp. 𝜇 ⟩︀ -ln 𝜌 𝑝 𝜇 = 2 . 3 × log 10 (︂ 𝑆 exp 𝜇 𝑆 sim 𝜇 ( 𝑝 ) )︂ = = -0 . 127 ± 0 . 035 , (33)  \nwhere the right side can be determined from the data in Fig. 4 and Fig. 5 by summation of all 127 events. For this sample with average energy ⟨ 𝐸 ⟩ ≈ 1 . 82 × 10 19 eV and ⟨ sec 𝜃 ⟩ ≈ 1 . 33 we will obtain the following value:  \n⟨ 𝑧 exp 𝐴 ⟩ = -0 . 127 ln 𝜌 Fe 𝜇 -ln 𝜌 𝑝 𝜇 = -0 . 03 ± 0 . 02 , (34)  \nwhich within errors agrees with our earlier estimations of this parameter obtained from mean LFDs [6-8]. Formally, one can conclude from (34) and formula (2)  \nthat in the considered energy range CRs consist mainly of protons. But in fact this is not the case. Estimations of nuclear composition of primary particles derived from the EAS muon content, on average, can be significantly distorted due to certain fraction of muon-poor and muon-rich EAS events presented in the data.", '4. CONCLUSION': "A new algorithm for estimation of the CR mass composition in individual EAS events is described above. It is based on a comparison of the measured MD responses and those calculated using the qgsj etii .04 hadron interaction model for primary particles with a given mass. In this particular case, protons with zenith angles 𝜃 ≤ 60 ∘ were considered. The chosen model accurately describes the development of all EAS components. It has proven to be highly effective during the investigation of the so-called 'muon puzzle' [7,8,21]. Calculations of SD and MD responses for this work were performed without resorting to the Monte Carlo method for individual events. For this purpose, zenith-angular dependences of these values were obtained beforehand (Fig. 1). This approach significantly simplified the analysis of the air showers muon content without losing the quality of the obtained results.  \nComparison of experimentally determined and calculated responses for all EASs was made for a given value of CR energy 𝐸 * = 3 . 16 × 10 19 eV. The measured MD responses in individual showers with energy 𝐸 were normalized by the corresponding calculated values according to relation (13). The essence of the method is demonstrated in Fig. 4, where correlations of the MD response densities at axis distances 𝑟 = 600 m and 1000 m are shown (such picture is observed at any axis distance). We selected events with densities determined at 1000 m(as in [23,24]), which give maximum statistics with a good quality of muon data. The relation between the experimentally derived muon response and the expected from the qgsj etii .04 model is shown in Fig. 5. This distribution has several pronounced peaks. The first one (29) undoubtedly refers to protons. The second peak (30) conventionally could be attributed to iron nuclei. Third one (31) is formed by several events with abnormally high muon content, including the largest air shower registered at the Yakutsk array (see Fig. 2 and Fig. 3). Its origin is yet to be discovered. The fourth peak (32) refers to abnormally muon-poor EAS events. Probably it has direct connection to primary gamma photons. Some characteristics of events in these peaks are presented in Table 1.  \nTable 1. The number of EAS events in energy bins with a width Δlog 10 𝐸 = 0 . 1 .  \nFurther plans are to continue studying the CR mass composition by applying this method to the Yakutsk EAS array data in the lower energy range.", 'ACKNOWLEDGEMENTS': "This work was made using the data obtained at The Unique Scientific Facility 'The D. D. Krasilnikov Yakutsk Complex EAS Array' (YEASA) ( https: //ckp-rf.ru/catalog/usu/73611/ ). Authors express their gratitude to the staff of the Separate structural unit YEASA of ShICRA SB RAS.", 'FUNDING': 'This work was made within the framework of the state assignment # 122011800084-7.', 'CONFLICT OF INTEREST': "The authors of this work declare that they have no conflict of interest.  \n- 1. J. R. Horandel, J. Phys.: Conf. Ser. 47 , 41 (2006).\n- 2. A. Aab, P. Abreu, M.Aglietta et al., Nucl. Instr. Methods A 798 , 172 (2015), arXiv:1502.01323 [astroph.IM].\n- 3. T. Abu-Zayyad, R. Aida, M. Allen et al., ApJ Lett. 768 , L1 (2013), arXiv:1205.5067 [astro-ph.HE].\n- 4. A. Aab et al. (Pierre Auger Collaboration), Phys. Rev. Lett. 117 , 192001 (2016), arXiv:1610.08509 [hep-ex].\n- 5. R. U. Abbasi et al. (Telescope Array Collaboration), Phys. Rev. D 98 , 022002 (2018), arXiv:1804.03877 [astro-ph.HE].\n- 6. H. P. Dembinski, J. C. Arteaga-Vel'azquez, L. Cazon et al., EPJ Web Conf. 210 , 02004, arXiv:1902.08124 [astroph.HE].\n- 7. A. V. Glushkov, A. V. Saburov, L. T. Ksenofontov et al, Phys. Atom Nucl. 87 , 25 (2024), arXiv:2306.17039 [astro-ph.HE].\n- 8. A. V. Glushkov, K.G.Lebedev, A.V.Saburov, JETP Lett. 117 , 257 (2023), arXiv:2304.09924 [astro-ph.HE].\n- 9. N. N. Kalmykov, S. S. Ostapchenko, A. I. Pavlov. Nucl. Phys. B /emdash.cyr Proc. Suppl. 52 , 17 (1997).\n- 10. S. Ostapchenko, Phys. Rev. D 83 , 014018 (2011), arXiv:1010.1869 [hep-ph].\n- 11. D. Heck, J. Knapp, J. N. Capdevielle et al. CORSIKA: AMonte Carlo Code to Simulate Extensive Air Showers . Forshungszentrum Karlsruhe, FZKA 6019. 90 p. (1988).\n- 12. A. V. Glushkov. Lateral distribution and full flux of Cherenkov light emission in EASs with primary energy 𝐸 0 ≥ 10 17 eV . PhD thesis, SINP MSU 1982 (in Russian).\n- 13. A. V. Glushkov, M. N. Dyakonov, T. A. Egorov et al, Izv. AN SSSR: Ser. Fiz. 55 , 713 (1991) (in Russian); B. N. Afanasiev, M. N. Dyakonov, T. A. Egorov et al. in Proc. of the Tokyo Workshop on Techniques for the Study of Extremely High Energy Cosmic Rays , edited by M.Nagano, p. 35 (1993).\n- 14. A. V. Glushkov, A.A.Ivanov, S. P. Knurenko et al., in Proc. of the 28th ICRC , vol. 1, p. 393, edited by T. Kajita, Y. Asaoka, A. Kawachi et al. (2003).\n- 15. A. V. Glushkov, M. I. Pravdin and A. Sabourov, Phys. Rev. D 90 , 012005 (2014), arXiv:1408.6302 [astroph.HE].\n- 16. A. V. Glushkov, M. I. Pravdin and A. V. Saburov, Phys. Atom. Nucl. 81 , 575 (2018), arXiv:2301.09654 [astroph.HE].\n- 17. L. I. Kaganov. Lateral distribution of charged particles in extensive air showers with primary energy above 𝐸 0 ≥ 10 17 eV . PhD thesis, SINP MSU 1981 (in Russian).\n- 18. M. I. Pravdin. Spectra obtained from particle fluxes of extensive air showers with energy above 10 17 eV . PhD thesis, SINP MSU 1985 (in Russian).  \n- 19. A. V. Saburov. Lateral distribution of particles in EASs with energy above 10 17 eV according to the data of Yakutsk Array . PhD thesis, INR RAS 2018.\n- 20. T. Pierog, Iu. Karpenko, J. M. Katzy et al. Phys. Rev. C 92 , 034906 (2015), arXiv:1306.0121 [hep-ph].\n- 21. A. V. Glushkov, A.V.Sabourov, L.T.Ksenofontov et al., JETP Lett. 117 , 645 (2023), arXiv:2304.13095 [astro-ph.HE].\n- 22. A. V. Glushkov, I. T. Makarov, E. S. Nikiforova et al., Astropart. Phys. 4 , 15 (1995); A. V. Glushkov, M.I. Pravdin, I. E. Sleptsov et al., Phys. Atom Nucl. 65 , 1313 (2002).\n- 23. A. V. Glushkov, I. T. Makarov, M.I. Pravdin et al., JETP Lett. 87 , 190 (2008), arXiv:0710.5508 [astro-ph].\n- 24. A. V. Glushkov, D. S. Gorbunov, I. T. Makarov et al., JETP Lett. 85 , 131 (2007), arXiv:astro-ph/0701245."}