A Monte Carlo study to check muon numbers in hadronic interaction models observed by the Tibet hybrid experiment (YAC-II + Tibet-III + MD)

The interpretation of the EAS data relies on the employment of hadronic interaction models, which are subject to theoretical and experimental uncertainties that may hamper composition studies of cosmic rays. To check the reliability of such models at the energies relevant for EAS studies, the predictions of the models can be compared with data from air shower observatories. In this regard, the study of the number of muons becomes extremely useful, since they are sensitive to the hadronic interactions that occur in the early phases of the EAS development. In the paper, we propose a Monte Carlo study to test the number of muons in hadronic interaction models by the hybrid experiment (YAC-II + Tibet-III + MD). For an air-shower event, the Tibet air-shower array (Tibet-III) provides the arrival direction and the air-shower size which are interrelated to primary energy, the Yangbajing Air shower Core detector (YAC-II) array measures the high energy electromagnetic particles in the very forward region so as to obtain the characteristic parameters of air-shower cores, at the same time, the underground MDs record the number of high-energy muons above 1 GeV. Since we can select proton events with high accuracy by YAC-II almost independently of the hadronic interaction models, the accompanying number of muons induced by proton events can be fed out. With the unique settle of YAC-II, our results show that the description of muon numbers in diﬀerent hadronic interaction models can be well systematics-checked to avoid the ambiguity of the primary cosmic-ray mass composition around the knee energy region by the Tibet hybrid experiment (YAC-II + Tibet-III + MD).


A
: The interpretation of the EAS data relies on the employment of hadronic interaction models, which are subject to theoretical and experimental uncertainties that may hamper composition studies of cosmic rays. To check the reliability of such models at the energies relevant for EAS studies, the predictions of the models can be compared with data from air shower observatories. In this regard, the study of the number of muons becomes extremely useful, since they are sensitive to the hadronic interactions that occur in the early phases of the EAS development. In the paper, we propose a Monte Carlo study to test the number of muons in hadronic interaction models by the hybrid experiment (YAC-II + Tibet-III + MD). For an air-shower event, the Tibet air-shower array (Tibet-III) provides the arrival direction and the air-shower size which are interrelated to primary energy, the Yangbajing Air shower Core detector (YAC-II) array measures the high energy electromagnetic particles in the very forward region so as to obtain the characteristic parameters of air-shower cores, at the same time, the underground MDs record the number of high-energy muons above 1 GeV. Since we can select proton events with high accuracy by YAC-II almost independently of the hadronic interaction models, the accompanying number of muons induced by proton events can be fed out. With the unique settle of YAC-II, our results show that the description of muon numbers in different hadronic interaction models can be well systematics-checked to avoid the ambiguity of the primary cosmic-ray mass composition around the knee energy region by the Tibet hybrid experiment (YAC-II + Tibet-III + MD).

K
: Performance of High Energy Physics Detectors; Simulation methods and programs

Introduction
The all-particle energy spectrum of primary cosmic rays is roughly described by a power-law / ∝ − over many orders of magnitude, with changing sharply from 2.7 to 3.1 at around a few PeV. Such a structure of the all-particle energy spectrum is called the 'knee', and the corresponding energy range is so called 'knee region' [1,2]. The special structure of the power law spectrum is considered to be closely related to the origin, acceleration and propagation mechanisms of cosmic rays in the Galaxy. Although the existence of the knee has been well established experimentally, there are still controversial arguments on its origin due to the lack of detailed information about the chemical composition around the knee [3,4]. The best ways to study the chemical composition are direct measurements of primary cosmic rays by balloon flights or satellites, but the energy range with sufficient event statistics is limited to 10 14 eV because of limited exposure time and small effective area. Therefore, the study of chemical composition around the knee region depends on ground-based experiments. However, ground-based experiments can only observe secondary cosmic ray particles, and rely on the Monte Carlo simulation to reconstruct primary cosmic ray information (such as cosmic ray composition, energy and direction, etc.). Thus the uncertainty of primary composition models and hadronic interaction models used in Monte Carlo simulation made the precise measurement of chemical composition at the knee region difficult (the so-called 'ambiguity' problem).
To overcome these problems, the Tibet AS experiment has recently upgraded the Tibet Emulsion Chamber (Tibet-EC) experiment and developed a new low threshold core detector named -1 -YAC (Yangbajing Air shower Core detector), which is capable of observing the high-energy core events with high statistics [5,6]. The new Tibet hybrid experiment (YAC-II + Tibet-III + MD) consists  of the YAC-II detector array (YAC-II), the Tibet air-shower array (Tibet-III) and an underground water-Cherenkov Muon-detector array (MD). The three types of detectors play different roles in the process of studying the energy spectrum of primary cosmic rays. For an air-shower event, the Tibet-III provides the arrival direction and the air-shower size which are interrelated to primary energy, the YAC-II can record the electromagnetic showers induced by high energy electrons and photons in the EAS forward region (later referred to as "high-energy core events"), which are very sensitive to identifying the main components of primary cosmic rays. At the same time, the underground MDs record high-energy muons above 1 GeV. A primary proton-induced Air Shower (AS) produces much fewer muons than a primary iron-induced AS does. The Tibet MD array thus enables us to separate different primary cosmic-ray nuclei by means of counting the muon numbers ( ) in an AS. However, the description of muon numbers in different hadronic interaction models was built based on some theoretical scenario that was checked by the hadron collider experiments. Indeed, because of the limitation of the beam pipe, especially at higher energy, very poor information on hadronic interactions was obtained. Therefore, some extrapolation was applied during the modeling of the hadronic interactions. This made the description of muon numbers in different hadronic interaction models have sizeable uncertainties theoretically as well as experimentally [7].
In the paper, we develop a method to check muon numbers in hadronic interaction models by a Monte Carlo study. In fact, it is very difficult to check the number of muons in the hadronic interaction models. At present, the accelerator experiment cannot check the number of muons in the hadronic interaction models above 100 TeV. The ground-based experiments are also difficult to do it. The reason is that the ground-based experiments do not only rely on the hadronic interaction models, but are also limited by the unknown measurement of the primary cosmic-ray component around the knee energy region. Once the primary cosmic-ray composition models are assumed in the Monte Carlo simulation, the uncertainty of the adoption of different primary cosmic-ray composition models is introduced. Therefore, if we check the number of muons in the hadronic interaction models, we should not only overcome the dependence of the hadronic interaction models, but also overcome the interference caused by the uncertainty of the primary cosmic-ray composition model around the knee energy region.
In this work, we use the unique advantages of the YAC-II experiment to solve this problem. The method is introduced as follows: 1) Firstly, we use the characteristic parameters observed by (YAC-II + Tibet-III) which are extremely sensitive to proton events, and using Artificial Neural Network (ANN), we can select proton events with high precision from MC simulation data and experimental data. 2) It is worthwhile to note that there is no serious difference among the current interaction models on the particle production in the forward region at the knee energy region since all the models are well tuned using recent accelerator data including LHC [8,9]. Therefore, we can separate primary proton events from other nuclei almost independently of the used hadronic interaction models by using the Tibet hybrid experiment (YAC-II + Tibet-III) [10]. This is the key to successfully check the muon numbers in this work. 3) At the same time, its accompanying is recorded by the Tibet-MD experiment. 4) Then, as long as we compare the differences in the number of muons associated with these selected proton events in the MC simulation and the Tibet-MD experiment, we can check which hadronic interaction model has the closest number of muons to the Tibet-MD experiment, -2 -and then the number of muons in the hadronic interaction model is correct. Since we have selected pure proton events with high accuracy in both MC simulation and experimental data, this method can avoid the interference caused by the uncertainty of the primary cosmic-ray composition models around the knee region assumed in MC simulation. 5) The key to the success of this work is that the YAC-II array plays an extremely important role, because YAC-II can well observe the EAS forward region, which not only overcomes the dependence of the ground-based experimental observation results on the hadronic interaction models, but also strictly limits the systematic error caused by the uncertainty of the primary cosmic-ray composition models at the knee region.
In this paper, we first introduced the Tibet hybrid experiment in section 2. Then, the Monte Carlo simulations including various combinations of hadronic interaction models and primary cosmic-ray composition models are described in section 3. Section 4 is the analysis including reconstructed multiple parameters observed by the (YAC-II + Tibet-III + MD) array. In section 5, we demonstrate the results and discussion to check muon numbers in hadronic interaction models observed by the Tibet hybrid experiment. Finally, the work and the results are summarised in section 6.

Experiment
The new Tibet hybrid experiment (YAC-II + Tibet-III + MD) has been operated in Tibet, China, since March, 2014. The current Tibet-III array, covering an area of 65,700 m 2 , consists of 569 Fast-Timing (FT) scintillation detectors [2] and 28 Density (D) detectors [2], as shown in figure 1. The FT detectors, each 0.5 m 2 in area, are placed on a lattice of 7.5 m spacing with D detectors surrounding them. By detecting the electromagnetic components in an AS, each detector of the Tibet-III array records the arrival times and densities of the detected particles. Then, the arrival direction and energy of each primary cosmic-ray event can be reconstructed from these data. The air-shower direction can be estimated with an error of about 0.2 • above 100 TeV, and the primary energy resolution is estimated to be 12% at energies around 10 15 eV by our simulation [2].
The YAC-II consists of 124 detectors, covering a total area of 500 m 2 , as shown in figure 1. It is aiming to distinguish the primary cosmic-ray nuclei through recording the electromagnetic showers induced by high energy electrons and photons in the EAS cores. As shown in figure 2, each detector of the YAC-II has a lead plate of 3.5 cm thickness (corresponding to 6.3 radiation length), an iron plate of 0.9 cm thickness and a plastic scintillator of 1 cm thickness. In order to record the electromagnetic showers in the energy range from several GeV to several 10 TeV, a wide dynamic range from 10 MIPs to 10 6 MIPs of PMT is required. Therefore, each of YAC-II detectors has two photomultiplier tubes, a high-gain PMT (HG-PMT, HAMAMATSU R4125) and a low-gain PMT (LG-PMT, HAMAMATSU R5325) that are responsible for the range of 1-3 × 10 3 MIPs and 10 3 -10 6 MIPs, respectively. At present, the YAC-II array is constructed near the center of the Tibet-III, and it has been operated simultaneously with Tibet-III and MD. When a YAC event is triggered, its accompanying air shower is simultaneously recorded. The matching between YAC, AS and MD event is made by their arrival time stamps. The inner 100 counters of 80 cm × 50 cm each are deployed in a 10 × 10 matrix form with a separation of 1.9 m. Another outer 24 counters of 100 cm × 50 cm each surround the inner detector matrix to reject showers whose cores fall far away from the YAC-II array. The details of the hardware of YAC-II are described in previous papers [11,12].  The MD array now consists of 4 pools set up 2.5 m underground, each with 16 cells, covering a total area of ∼ 3,400 m 2 , as shown in figure 1. 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 [13], equipped with two downward-facing 20-inch-in-diameter PMTs (HAMAMATSU: R3600) on the ceiling, and the inside of each cell is -4 -painted with white epoxy resin for waterproof and efficient reflection of the water-Cherenkov light. The electromagnetic component is shielded by the soil overburden corresponding to ∼ 19 radiation lengths, while the energy threshold for muons to penetrate the soil is approximately 1 GeV.

Monte Carlo simulations
A full Monte Carlo simulation has been carried out on the development of EAS in the atmosphere and the response in the Tibet hybrid experiment (YAC-II + Tibet-III + MD) array. The simulation code CORSIKA (version 6.204) including QGSJET01c and SIBYLL2.1 hadronic interaction models and CORSIKA (version 7.7410) including SIBYLL2.3d and EPOS-LHC hadronic interaction models are used to generate air-shower events with primary energy above 50 TeV [14]. In this work, three primary composition models (helium poor model (He-poor) [15], helium rich model (He-rich) [16] and Gaisser's fit model (H4a model)) [17] are used due to the uncertainty of the primary chemical composition model. The fractions of the component of the three composition models in different energy regions are listed in table 1. Details of primary composition models are described in the previous paper [6]. Primaries are injected into the top of the atmosphere isotropically within the zenith angle from 0 • to 60 • . Secondary particles are traced to the altitude of 4300 m till 1 MeV. The generated secondary particles are fed into the detector simulation of the (YAC-II + Tibet-III + MD) array developed by using the GEANT4 code (version v4.10.00) [18]. As mentioned above, the generated 4 hadronic interaction models under 3 primary composition models are combined in order to check the uncertainties caused by the use of different primary cosmic-ray composition models. In the paper, we used 6 combinations of hadronic interaction models and primary composition models whose details are described in section 5. The statistics of generated primary cosmic-ray events are listed in table 2. And the number of selected core events satisfying the criteria described in the next section is also shown in table 2, respectively.

Analysis
In our hybrid experiment, the number of charged particles detected by each detector is defined as the PMT output (charge) divided by that of the single peak, which is determined by a probe calibration. According to the MC, the single peak value of the energy deposit for a single particle in each detector of YAC-II and Tibet-III is calculated as 1.98 MeV and 5.6 MeV, and the similar single peak value of photon yield of MD is 23.6 photons (for the details, please see [6,19]). Therefore, we can estimate the number of charged particles for each detector. The number of charged particles under the lead plate of a YAC detector ( , )(later referred to as "burst size") [6], detected particle density of a Tibet-III detector ( ) and the number of detected muon ( , ) can be calculated by the total energy deposit (or photons yield) of each primary cosmic ray divided by the corresponding value of single peak. Based on the basic quantities, we can reconstruct multiple parameters of the (YAC-II+Tibet-III+MD) array used to characterize an air-shower core event.

Reconstructed parameters by (YAC-II+Tibet-III+MD)
Based on the basic quantity , of YAC-II detector, we can get the following parameters from YAC-II array: • Nhit YAC : the number of 'fired' YAC-II detector units with ≥ a given threshold value; • top : the maximum burst size among the 'fired' YAC-II detectors; • : the total burst size of all 'fired' YAC-II detector units; • are the burst size in the th fired YAC-II detector unit and the lateral distance from the burst center, which is calculated as weighted center of gravity, to the center of the th fired detector, respectively.
Based on the basic quantity of Tibet-III detector, we can get the following parameters from Tibet-III array: • : the arrival direction of the air shower; • : the air-shower size, which is proportional to the energy of the primary cosmic ray; • core , core : reconstructed core location: ( core , core ) = ( Σ Σ , Σ Σ ), where denotes the particle density at the th detector, and the weight is an energy-dependent parameter; • AS : the distance from the reconstructed core location( core , core ) to the center of the YAC-II array; • AS : the lateral distance from the 'fired' Tibet-III detector units to the reconstructed core location.
Based on the basic quantity , of MD detector, we can get the muon numbers ( ) from MD array: • : the sum of muon numbers measured by all the 'fired' MD detector units.
In order to obtain proton events and reject events falling far from the YAC-II array, the final event selection criteria are: Nhit YAC ≥ 4, ≥ 3000, ≥ 10 5 , 0 ≤ 15, where 0 is the distance from the reconstructed core location to each side of the YAC-II array. The cut on 0 was chosen to reduce the influence of systematic effects from the edges of the YAC-II array.

Behavior of high-energy core events
After the above process, we studied sensitive parameters to separate different primary cosmic-ray nuclei obtained by (YAC-II+Tibet-III) array. The merit of YAC detector is to distinguish the primary cosmic-ray nuclei through recording the electromagnetic showers induced by high energy electrons and photons in the EAS forward region. Among primary particles, proton with the long interaction mean free path can penetrate deep in the atmosphere and produce air-shower cores near the observation level, resulting in a larger burst size and smaller lateral spread. Thus, the burst size of high-energy core events originated by proton is relatively large, while iron-induced event shows a smaller value. Figure 3(a) shows the behavior of the maximum burst size( top ) which is normalized by air-shower size ( ) based on SIBYLL2.1+He-poor model. It can be seen that the value of top produced by proton is larger than iron due to its inelastic cross section, indicating that the YAC-II array is capable to clarify chemical composition. Figure 3(b) is the correlation between AS and / which also indicates that the parameter of observed by YAC-II is quite diverse for different primary nuclei. It means that (YAC-II+Tibet-III) array can be used to separate proton events from other nuclei by detecting high-energy core in the EAS forward region. Therefore, we can check muon numbers simultaneously recorded by Tibet-MD with the selected proton samples in this work. In the following, we carried out the method to check the description of by the Tibet hybrid experiment (YAC-II + Tibet-III + MD).

Selection of proton-induced events by (YAC-II+Tibet-III) array
Firstly, we studied the separation of primary mass composition by (YAC-II+Tibet-III) array by using a feed-forward ANN method, whose applicability is shown to be quite effective as confirmed by our previous works [20]. In this ANN analysis, we use the following 8 quantities: (1) Nhit YAC , (2) top , (3) , (4) , (5) , (6) , (7) and (8) AS . We have also tested to include other parameters or the combination of used parameters as ANN-input, however, there was no significant improvement in the mass separation. To train the ANN in separating proton from other nuclei, the input patterns for proton and others are set to 0 and 1, respectively. Then, the ANN output pattern value (T) is a real number from 0 to 1. Figure 4 shows the ANN output distribution trained using the SIBYLL2.1+He-poor model. One can see that the proton-induced events are clearly separated from other nuclei. As seen in this figure, the events with ≤ 0.3 could be regarded as the proton-like events, and the average selection purity (p) and efficiency ( ) over the whole energy range ( ≥ 50 TeV) of proton-like events are ∼ 81% and 43%, where purity (p) represents the ratio of true proton events contained in the ANN trained core events with ≤ = 0.3, efficiency ( ) represents the ratio between true proton events with ≤ = 0.3 and all true proton events in the sample. We have also compared the ANN test results by (YAC-II+Tibet-III) array under the other 3 hadronic interaction models and found that the difference in the purity (p) is less than 5%, which indicates that the difference among the hadronic interaction models on the particle production in the forward region was not serious, as measured by LHC [8,9].

Overcoming the influence of primary cosmic-ray composition models
From above, we have studied the ability of mass separation by the ANN method based on SIBYLL2.1+He-poor model. Here, we fix SIBYLL2.1 interaction model as an example to discuss the effect of different cosmic-ray composition models on the check of muon numbers. The 3 combinations of SIBYLL2.1+He-poor, SIBYLL2.1+He-rich and SIBYLL2.1+H4a are fed to ANN to separate the proton events from others individually. After the training of ANN, we can obtain proton-like events with high accuracy under different primary cosmic-ray composition models, so as to strictly limit the systematic error caused by the uncertainty of the primary cosmic-ray composition models. At the same time, the accompany are fed out. Figure 5(a) shows the scatter plots between and on each proton-like event based on SIBYLL2.1+He-poor, SIBYLL2.1+He-rich and SIBYLL2.1+H4a model with zenith angles smaller than 40 • , respectively. In order to distinguish the difference more clearly, we fit the Gauss distribution of at different bins. Figure 5(b) shows the distribution between fitted by the Gaussian function and the air-shower size of proton-like events under different composition models. The errors bars indicate the fitted errors. Surprisingly, it is found that the -distribution of proton-like events is consistent within the error range of 5% after ANN training under different primary cosmic-ray composition models based on SIBYLL2.1 interaction model. In fact, we have also tested the effect of different primary cosmic-ray composition models based on other hadronic interaction models and obtained similar conclusions. That is, the effect of primary cosmic-ray composition models is very small (within 5%) when we fix a certain hadronic interaction model due to the excellent mass separation ability of (YAC-II+Tibet-III) and the YAC-II array also has an excellent ability to observe high-energy core events in the EAS forward region.
The reason for this result is that there is no serious difference among the current interaction models on the particle production in the EAS forward region at the knee energy region since all the models are well tuned using recent accelerator data including LHC [8,9]. Therefore, the YAC-II array plays the most important role in this work. Because YAC-II can well observe the EAS forward region, which not only overcomes the dependence of the ground-based experimental observation results on the hadronic interaction models, but also strictly limits the systematic error caused by the uncertainty of the primary cosmic-ray composition models at the knee region, so as to successfully check the number of muons in the hadronic interaction models by using the Tibet hybrid experiment (YAC-II + Tibet-III + MD). That is, we can check muon numbers in hadronic interaction models with proton-like events almost independently of the used primary composition models. In the following, we study the -distribution of proton-like events in 4 hadronic interaction models to check muon numbers under He-poor composition model as an example.

Check muon numbers in hadronic interaction models
In the same way, the 4 combinations of QGSJET01c+He-poor, SIBYLL2.1+He-poor, SIBYLL2.3d+Hepoor and EPOS-LHC+He-poor are fed to ANN to separate the proton events from others individually. After the ANN training, we can also obtain the distribution between fitted by the Gaussian function and the air-shower size of proton-like events with zenith angles smaller than 40 • under different hadronic interaction models, as shown in figure 6(a). It is seen from this figure that there exists some differences among the 4 hadronic interaction models. In order to investigate the difference more clearly, the -distributions obtained by 4 MC models are compared to each other by taking the ratio to that by the SIBYLL2.1+He-poor model in figure 6(b). The red solid line (Ratio = 1) denotes SIBYLL2.1+He-poor model. It is seen that the muon numbers of proton-like events generated by SIBYLL2.3d+He-poor model gives the highest ratio which indicates that the muon numbers contained in EAS described by SIBYLL2.3d is ∼ 20% larger than SIBYLL2.1. The muon numbers of proton-like events described by EPOS-LHC is ∼ 10% larger than SIBYLL2.1. Although the difference between QGSJET01c and SIBYLL2.1 is about 5%, we can investigate the distribution of of proton-like events at certain fixed energy region, as shown in figure 7(a). Some discrepancies can be seen in the shape of distribution between QGSJET01c and SIBYLL2.1 at around 1 PeV (the knee region). Figure 7(b) and figure 7(c) shows the similar distribution of of proton-like events when we compare SIBYLL2.3d and EPOS-LHC with SIBYLL2.1 at around 1 PeV (the knee region), respectively. This suggests that if we further impose Tibet-MD experimental data on the correlation of -and the distributions of at certain fixed energy region, the muon numbers described in different hadronic interaction models can be well systematics-checked with the selected proton-like events to some extent.

Summary
In this paper, we develop a method to check muon numbers by the Tibet hybrid experiment (YAC-II + Tibet-III + MD) to avoid the ambiguity of the cosmic-ray mass composition around the knee energy region. With the unique settle of YAC-II detectors, we can distinguish the proton events with high precision by detecting the high-energy core events in the EAS forward region in coincidence with an accompanied air shower by (YAC-II + Tibet-III). On the other hand, there is no serious difference among the current interaction models on the particle production in the forward region at the knee -11 -energy region since all the models are well tuned using recent accelerator data including LHC [8,9]. Therefore, we can separate primary proton events with high precision from other nuclei almost independently of the used hadronic interaction models [10]. At the same time, its accompanying is simultaneously recorded. With the selected proton samples, the description of among current hadronic interaction models can be well systematics-checked to avoid the ambiguity of the cosmic-ray mass composition. Our results show that, the muon numbers of proton-like events generated by SIBYLL2.3d+He-poor model gives the highest ratio which indicates that the muon numbers contained in EAS described by SIBYLL2.3d is about 20% larger than SIBYLL2.1. The muon numbers of proton-like events described by EPOS-LHC is about 10% larger than SIBYLL2.1. Although the difference between QGSJET01c and SIBYLL2.1 is about 5%, there exsist some discrepancies in the shape of distribution at around 1 PeV (the knee region). The key to the success of this work is that YAC-II array plays an extremely important role, because YAC-II array can well observe the EAS forward region, which not only overcomes the dependence of the ground-based experimental observation results on the hadronic interaction models, but also strictly limits the systematic error caused by the uncertainty of the primary cosmic-ray composition models at the knee region.