Towards a better estimation of energy and species of primary cosmic rays in the knee region with the Tibet hybrid experiment: Utilization of refined EAS lateral distributions

A new hybrid experiment (YAC-II+Tibet-III+MD) located at Yangbajing has started to improve the capability of explicit measurements of the cosmic-ray components (P, He, CNO, Fe, etc) at the knee energy region since 2014. Considering the different features of air shower development for different primary cosmic-ray nuclei, using a full Monte Carlo simulation, we modified the lateral distribution functions for proton-induced, helium-induced and iron-induced air showers, respectively. The results show that the air shower size of different nuclei obtained by the modified Nishimura-Kamata-Greisen function is consistent with the true shower size within 5% systematic errors. Furthermore, we studied the Extensive Air Shower age parameter determination from the lateral distribution functions of charged particles near the air shower core with the Tibet hybrid experiment. The derived age parameter can be well used to estimate the particle type of the incident primary cosmic rays.


Introduction
The spectra of cosmic rays (CR) can be well described by a power law dj/dE ∝ E −γ over a very broad energy range, separated with some "spectral features". The features can be softenings or hardenings of the spectrum, and appear as "knee-like" or "ankle-like" in the usual log-log graphic representation of the spectrum [1,2]. The most prominent and well known feature in the all particle spectrum is the "knee" at E 4 PeV. That is, the power index suddenly steepens from approximately -2.7 to -3.1 at around 4 PeV, resulting in a distinctive "knee" shape in the spectrum [1], and the corresponding energy range is the so called "knee region". The identification and study of its origin is of great importance to understand the astrophysical mechanisms of acceleration and propagation [3,4]. However, there is still controversy to its origin. For this explanation, many authors have proposed various models such as a change of acceleration mechanisms at the sources of cosmic rays (supernova remnants, pulsars, etc.), an assumption of nearby sources emitting high energy cosmic rays, a change of cosmic ray propagation in the Galaxy (diffusion, drift, escape from the Galactic disk), and some unknown new processes in the atmosphere during air-shower development [5][6][7][8]. Hence, in order to resolve the origin of the knee, it is critical to measure the primary chemical composition or mass group at energies of 50 TeV-10 16 eV, especially, to measure the primary energy spectra of individual * e-mail: yingzhang@ihep.ac.cn components and determine a break energy of the spectral index for individual components [3,9,10].
A new hybrid experiment (YAC-II+Tibet-III+MD) located at Yangbajing with the advantageous geographical conditions (4300 m above sea level (a.s.l.), atmospheric depth: 606 g/m 2 ), has started to improve the capability of explicit measurements of the cosmic-ray components (P, He, CNO, Fe, etc) at energies from 10 14 to 10 16 eV. Considering the different features of air shower (AS) development for different primary cosmic-ray nuclei, in this paper, we firstly modified the lateral distribution functions (LDFs) for proton-induced, helium-induced and iron-induced air showers, respectively. Then, based on the modified Nishimura-Kamata-Greisen (NKG) function, we studied Extensive Air Shower (EAS) age determination from the LDFs of charged particles near the air shower core with the Tibet hybrid experiment to better estimate the particle type of the incident primary cosmic rays.

The new hybrid experiment
In 2014, an important upgrade of the Tibet ASγ project was completed, and we then started a new hybrid experiment consisting of an air-shower-core array (YAC-II), a high-density air-shower array (Tibet-III) and a large underground water-Cherenkov muon-detector array (MD). The Tibet hybrid experiment is not only a simple combination of three type of detectors, but also hybrid measurements of electromagnetic and hadronic components of the air shower. For an air shower event, Tibet-III provides the arrival direction (θ) and the air shower size (Ne f it ) which are interrelated to primary energy, YAC-II is used to record high energy core events, both of the two arrays measure the electromagnetic components generated by air showers, while MD is used to measure the high-energy (above 1 GeV) µ component.
The Tibet-III array (∼50000 m 2 ) at present consists of 761 fast timing (FT) counters with 28 density (D) counters around them. In the inner 36,900 m 2 , the FT counters are deployed at 7.5 m lattice intervals, and 249 counters among the 761 FT counters are equipped with a density-PMT (D-PMT) too.
The YAC-II array is located near the center of the Tibet-III air-shower array. YAC-II consists of 124 detectors, covering a total area of ∼500 m 2 . Each detector of YAC-II comprises a plastic scintillator with the size of 80 cm × 50 cm × 1 cm and a lead plate with a thickness of 3.5 cm (6.3 radiation lengths) being put on the scintillator. Each YAC detector is viewed by two photomultipliers (PMT) of high-gain (HAMAMATSU: R4125) and lowgain (HAMAMATSU: R5325) to cover the wide dynamic range from 1 MIP (Minimum Ionization Particle) to 10 6 MIPs.
The MD array now consists of 5 pools set up 2.5 m underground, each with 16 cells, covering a total area of ∼4500 m 2 . Each cell of the MD array is composed of a concrete water tank 7.2 m wide × 7.2 m long × 1.5 m deep, equipped with two downward-facing 20-inch-in-diameter PMTs (HAMAMATSU: R3600) on the ceiling, and the inside of each cell is painted with white epoxy resin for waterproof and efficient reflection of the water-Cherenkov light.

Simulation
A full Monte Carlo (MC) simulation has been carried out on the development of extensive air showers in the atmosphere and the response in the Tibet hybrid experiment (YAC-II+Tibet-III+MD) array. The simulation code CORSIKA including QGSJETII and SIBYLL2.1 hadronic interaction models is used to generate air-shower (AS) events with primary energy above 50 TeV [11]. Considering the main purpose of this work is to study EAS age parameter to separate iron nuclei, only the He-poor primary composition model is used as the input energy spectra mentioned in the paper [1]. Primaries isotropically incident at the top of the atmosphere within the zenith angles from 0 to 60 degrees are injected into the atmosphere. The minimum primary energy of this simulation is set at 1 TeV. Secondary particles are traced to the altitude of 4300 m till 1 MeV. The air-shower events are randomly dropped onto the YAC-II detectors array plane, 25 m wider in each side of the YAC-II array. In order to treat the MC events in the same way as experimental data analysis, simulated air-shower events were input to the detector with the same detector configuration as the (YAC-II+Tibet-III+MD) array. All detector responses are calculated using Geant4 (version 9.5) [12].

Analysis
As the particle lateral density distribution of extensive air showers is the key quantity for most of the ground-based cosmic-ray experiments, the LDFs of EAS play an important role in the analysis of the air-shower events. It can be used to derive air-shower size and age parameter (s ), etc, which, in turn, can be used to estimate the energy and particle type of the incident primary cosmic rays. Usually, we can use the following Nishimura-Kamata-Greisen function Here B denotes the beta function, r is the distance from air shower axis, N e and s are free parameters representing the air shower size and the age of the air shower (s corsika ), respectively. In the original NKG function, a(s) = s − 2, b(s) = s − 4.5 are used and r m (Moliere unit) = 130 m. However, the original NKG function is not completely suitable for the Tibet-III array because there is a 5 mm thick lead plate placed above each scintillation counter. In a previous paper, we modified the NKG function which can be well fitted to the lateral distribution but for allparticle shower particles under a lead plate of 5 mm thickness, where r m =30 m. The details of this work have been described in [1]. In this paper, considering the different feature of AS development for different primary cosmicray nuclei, the lateral distribution of air showers for proton, helium or iron primaries may be different from that for the all particle spectrum, and we try to develop LDFs for proton, helium and iron components, respectively, for the current new Tibet ASγ hybrid experiment. The function of a(s) and b(s) was developed as follows. The MC events were classified by s corsika and the averaged lateral distribution (normalized by total number of electrons) at a given s corsika was fitted by Eq. (1) (r m =30 m) to obtain 'a' and 'b' as constants. Then the dependence of 'a' and 'b' on s corsika was obtained as a(s) and b(s). The functions a(s) and b(s) are determined for the case of pure protons, pure helium and pure iron, respectively. Fig. 1(a) shows a comparison of a(s) for protons, helium and iron nuclei. Fig. 1(b) is the same one for b(s). Proton, helium and iron show different behaviors when s corsika >1.2. Although our result shows dependencies of a and b on s corsika of different nuclei that differ from each other, it is confirmed that the lateral distribution of the shower particles is well reproduced by our modified NKG function. When reconstructing the air shower size (Ne f it ), the age parameter, is also obtained by fitting the lateral density distribution of air showers. As can be seen in Eq. (1), the reconstructed age parameter has the nature of the longitudinal age s corsika although it was obtained from single observation level.
Based on the MC simulation, the correlation between the true shower size (Ne true ) and the estimated shower size (Ne f it ) (Proton, Helium and Iron, respectively) is demonstrated in Fig. 2. Here, the true shower size (red lines) means the number of particles calculated for a carpet array, while the estimated shower size (black open circles) is for the real Tibet-III array using the modified NKG function induced by different nuclei mentioned above. As seen in these figures, a good correlation (∼5%) is obtained between the true shower size and the estimated shower size. However, there are ∼10% (for smaller zenith angles) and ∼20% (for larger zenith angles) lower estimate of N e by the original NKG function (blue triangles) than Ne true . The resolution of Ne f it using the modified NKG function is ∼9% and ∼5% at around 300 TeV and 1 PeV, respectively, based on the QGSJETII+He-poor model with 1.0<secθ≤1.1.

Results and Discussion
In this work, we studied the separation of the mass composition by the (YAC-II+Tibet-III) and (YAC-II+Tibet-III+MD) arrays, respectively, by using a feed-forward artificial neural network (ANN) method, whose applicability to our experiment was well confirmed by the MC simulation [9]. Here, a three-layered feed-forward network was applied to select iron-induced events from others.

Selection of iron-induced events by (YAC-II+Tibet-III) array
Firstly, we investigated the separate ability of iron by (YAC-II+Tibet-III) array, both of the two arrays measure the electromagnetic components generated by air showers.
The following 10 parameters are input to the ANN with 50 hidden nodes and 1 output unit: (10) s . To train the ANN in separating iron from other nuclei, the input patterns for irons and others are set to 1 and 0, respectively. In this paper, the ANN output value T is used to separate the primary nuclei, and we define a critical value of T c to calculate the purity and the selection efficiency of the selected nucleus events. The purity (p) and selection efficiency ( ) of iron are calculated using Eq. (3).
The events with T ≥ T c are regarded as iron group. We accomplished the ANN training of iron under two hadronic interaction models (QGSJETII, SIBYLL2.1). The average purity and selection efficiency over the whole energy range is: 80.0%, 29.4% for (QGSJETII + Hepoor), and 79.0%, 27.9% for (SIBYLL2.1 + He-poor) at T c ≥ 0.7. Fig. 3(a) shows the ANN test results for iron based on (QGSJETII + He-poor) model while the events with T≥ 0.7 are regarded as iron-like events. One can see that the ANN method works well to separate the iron group from other nuclei almost independently of the used hadronic interaction models. However, without the reconstructed s parameter, the average purity and selection efficiency is: 69.3%, 29.5%, for (QGSJETII+ He-poor), and 67.9%, 28.2% for (SIBYLL2.1 + He-poor) when we set almost the same efficiency, which indicate that s plays an important role in the iron separation from other nuclei.

Selection of iron-induced events by the (YAC-II+Tibet-III+MD) array
Then we added reconstructed parameters from the MD array as input to the ANN. In the same way, the following 14 parameters are input to the ANN with 70 hidden nodes and 1 output unit:   Fig. 3(b) is the ANN test result for iron by the (YAC-II+Tibet-III+MD) array based on the QGSJETII+He-poor model while the events with T≥ 0.7 are regarded as the iron-like events. The average purity and selection efficiency over whole the energy range is about 86.2%, 65.7% for (QGSJETII + He-poor) and 85.2%, 62.3% for (SIBYLL2.1 + He-poor) at T c ≥ 0.7, respectively. We can see that the MD array is very helpful to separate the primary iron group from other nuclei. Considering the uncertainties of µ multiplicities among current hadronic interaction models, we will use (YAC-II+Tibet-III) array and (YAC-II+Tibet-III+MD) array to obtain primary iron spectrum in the analysis of experimental data in the future, respectively.

Summary
In this paper, using a full Monte Carlo simulation, we modified the lateral distribution functions (LDFs) for proton-induced, helium-induced and iron-induced air showers, respectively. The results show that the air shower size (Ne f it ) obtained by the modified NKG function is consistent with the true shower size (Ne true ) within 5% systematic errors. Based on the modified NKG function, we derived an age parameter which can be well used to estimate the particle type of the incident primary cosmic rays. The result shows that the constructed parameters work well to separate iron nuclei around and beyond the knee.