Hadronic interaction model dependence in cosmic Gamma-ray flux estimation using an extensive air shower array with a muon detector

Observation techniques of high-energy gamma rays using air showers have remarkably progressed via the Tibet ASγ, HAWC, and LHAASO experiments. These observations have significantly contributed to gamma-ray astronomy in the northern sky’s sub-PeV region. Moreover, in the southern sky, the ALPACA experiment is underway at 4,740 m altitude on the Chacaltaya plateau in Bolivia. This experiment estimates the gamma-ray flux from the difference between the number of on-source and off-source events by real data, utilizing the gamma-ray detection efficiency calculated through Monte Carlo simulations, which in turn depends on the hadronic interaction models. Even though the number of cosmic-ray background events can be experimentally estimated, this model dependence affects the estimation of gamma-ray detection efficiency. However, previous reports have assumed that the model dependence is negligible and have not included it in the error of gamma-ray flux estimation. Using ALPAQUITA, the prototype experiment of ALPACA, we quantitatively evaluated the model dependence on hadronic interaction models for the first time. We evaluate the model dependence on hadronic interactions as less than 3.6 % in the typical gamma-ray flux estimation performed by ALPAQUITA; this is negligible compared with other uncertainties such as energy scale uncertainty in the energy range from 6 to 300 TeV, which is dominated by the Monte Carlo statistics. This upper limit of 3.6 % model dependence is expected to apply to ALPACA.

gamma rays. Hadronic particles, such as protons, nuclei, and pions, are largely produced by collisions with the atmospheric nuclei when the incident particle is a hadronic cosmic ray. Collisions are repeated in the atmosphere, causing a hadronic cascade shower that amplifies the number of particles. Knowledge of hadronic particle interactions in the atmosphere is imperative in the Monte Carlo simulation process. The Monte Carlo simulation is based on various phenomenological models. Currently, the differences between these models are not negligible, and the number of muons contained in a hadronic air shower differs depending on the phenomenological model. Thus, differences in hadronic interaction models can cause systematic uncertainties in the detection efficiency of the gamma-ray induced air showers using the number of muons as a selection criterion. In order to avoid this dependence, Tibet ASγ experimentally estimates the background hadronic cosmic rays using air showers from directions away from the gamma-ray point source to measure flux [4].
Conversely, when the primary particle is a gamma ray or an electron, the air shower component is a mixture of gamma rays, electrons, and positrons due to the repetition of electron-positron pair creation and their bremsstrahlung, causing an electromagnetic cascade shower. However, even in this case, a small number of muons are produced through photonuclear interactions and muon-pair creation, among others. Pions are produced in gamma-ray-induced air showers via photonuclear reactions, which then interact in the atmosphere, undergoing hadronic interactions. The charged pions decay into muons. Additionally, muons are also produced in muon-pair production during electromagnetic interactions. The crosssection for photonuclear reactions is roughly two orders of magnitude larger than that of muon pair production in the energy range of GeV to TeV. This means that the muon component produced by photonuclear reactions is significantly more dominant than that from the muon pair production process in gamma-ray-induced air showers. Therefore, the muon component in the gamma-ray-induced air showers could depend on hadronic interaction models. However, the hadronic interaction model dependence of the gamma-ray-induced air showers is believed to be small and may be ignored, then the model dependence has yet to be quantified.
Herein, we quantitatively evaluate the model dependence of gamma-ray flux estimation for the first time. Section 2 clarifies the characteristics of air shower muons generated by a Monte Carlo simulation with some interaction models and their model differences using a vertical gamma-ray-induced air shower at the ALPACA altitude. Section 3 evaluates the systematic differences in the detection efficiency of gammaray-induced air showers caused by the γ -CR separation in the energy range from a few TeV to several hundred TeV with a small air shower array (ALPAQUITA [18]), which is the prototype of ALPACA [16,17].

Characteristics of air shower muons at high altitudes
ALPACA and Tibet ASγ measure the muon component above ∼1 GeV using an underground muon detector for γ -CR separation [4,18]. Therefore, we investigated the characteristics of the muons above 1 GeV in gamma-ray-induced air showers depending on hadronic interaction models.

Air shower simulation and simulation setting
We simulated air showers induced by vertical incident gamma rays with energies between 10 TeV and 100 TeV, using the CORSIKA 7.6000 code [26]. It employs the EGS4 code [27] to simulate electromagnetic interactions. The calculation of the photonuclear reactions was added to the EGS4 code [26]. The muon-pair production is also incorporated in the EGS4 code using the analogy of electron-positron pair production [26,28]. There is no definitive model for calculating the behavior of hadronic interactions. Thus, even in gamma-ray-induced air showers, the muon features are model-dependent.
A summary of the simulation setting is presented in Table 1. The secondary gamma rays and electrons are tracked down to 1 MeV while the secondary hadrons and muons are tracked to 1 GeV, the minimum kinetic energy above which they can reach the water surface of MD [18]. In addition, the particle identification, position, energy, timing, and directional vector of each secondary particle were recorded at an altitude of 4,740 m.

Model dependence of the energy spectrum
The calculated differential energy spectrum of muons >1 GeV in gamma-rayinduced air showers for each model combination is shown in Fig. 1a for 10 TeV and Fig. 1b for 100 TeV, respectively. For E γ = 10 TeV, the energy spectrum of muons peaks at about a few GeV. Approximately 85 % of the total number of muons are below 10 GeV. The energy spectrum shows a decrease above 100 GeV with a powerlaw index of −2.62, containing 99 % of the total number of muons below 100 GeV. The trend is similar for E γ = 100 TeV.
To evaluate the differences in the spectra caused by hadronic interaction models, we used the relative change R(F i ) in the differential spectrum to the QGSJET-II + FLUKA model using the following equation: where F i is an arbitrary physical quantity, i = 1, 2, 3, and 4 represent QGSJET-II + FLUKA, QGSJET-II + UrQMD, EPOS LHC + FLUKA and SIBYLL + FLUKA, respectively. For 10 TeV, the value of R for muons below a few GeV is 3.6 × 10 −2 or less. The model differences tend to be larger as the muon energy increases, and the EPOS LHC + FLUKA model has the most significant difference of 1.3×10 −1 at 100 GeV. We can test the two low-energy models by comparing the results of QGSJET-II + FLUKA and QGSJET-II + UrQMD for E γ = 10 TeV. The difference between the two sets of results is less than ∼ 2.6 × 10 −2 below 100 GeV.
In addition, the trend of R for E γ = 100 TeV is almost identical to that of E γ = 10 TeV for QGSJET-II + FLUKA, QGSJET-II + UrQMD, EPOS LHC + FLUKA and SIBYLL + FLUKA, and the differences below 10 GeV are typically of a few percentages. The integral energy spectra show the effects of these low-energy muon differences more clearly. The results of the energy spectrum of muons above 1 GeV in gamma-ray-induced air showers of 10 and 100 TeV are shown in Fig. 2. The majority of the total number of muons (over 85 %) are below 10 GeV as shown by the bottom panels in Fig. 2. The difference between the models in several GeV regions is small  ( Fig. 2). The difference in the integrated spectrum above 1 GeV is less than 2.2 % for the 10 TeV gamma-ray shower and less than 3.1 % for 100 TeV.

Model dependence of the lateral spread of muons
The ALPAQUITA underground muon detector has an area of 900 m 2 [18], which measures a portion of the muons in a widely spread air shower. Tibet ASγ and ALPACA have significantly larger areas than ALPAQUITA. For γ -CR separation, however, all three experiments use the muon density within a few hundred meters of the air shower core. Figure 3 shows the total number of muons above 1 GeV within a certain radius from the shower core, Fig. 3a for E γ = 10 TeV and Fig. 3b for E γ = 100 TeV. The number of muons in a 100 m radius circle centered on the air shower core is ∼ 9.6 × 10 −1 for 10 TeV gamma rays and 15 muons for 100 TeV gamma rays, depending on gamma-ray energy. Independent of energy, the distribution shows the same trend in all models.
The Rs relative to the QGSJET-II + FLUKA model, calculated using the same form as (1) is shown in Fig. 3. The number of muons produced by the QGSJET-II + UrQMD is smaller than that of QGSJET-II + FLUKA. However, the difference is small, about 2.0 × 10 −2 over the whole distance range. As the distance from the core approaches zero, the numbers of muons produced by the SIBYLL + FLUKA and the EPOS LHC + FLUKA becomes smaller than that by QGSJET-II + FLUKA. Furthermore, the maximum Rs is 8.9 × 10 −2 for both energies. However, when comparing the total numbers of muons within 100 m (∼Molière unit), the Rs at 100 m are small, ±7.7×10 −3 for 10 TeV gamma rays and ±1.8×10 −2 for 100 TeV gamma rays. The differences in the number of muons due to the hadronic interaction model options are only a few percent if we locate a large muon detector array within a circular area of 300 m in radius

Model dependence of the total number of muons for an air shower
ALPAQUITA uses only the total number of measured muons, but not the lateral distribution for the γ -CR separation. We investigated the total number of muons over 1 GeV per air shower to clarify the model dependence.
The upper figures in Fig. 4a and b show the total number of muons over 1 GeV per air shower for 10 and 100 TeV gamma rays, respectively. The leftmost bin in Fig. 4a for E γ = 10 TeV corresponds to less than 10 muons, and the number of events in the bin constitutes the large majority (∼90 %). For E γ = 100 TeV, the leftmost bin corresponds to less than 100 muons and contains ∼90 %.
As the typical gamma-ray efficiency in ALPAQUITA is assumed to be 50 % to 90 %, the hadronic model dependence of the leftmost bin in Fig. 4 is significant in this work. The lower figures in Fig. 4 show the Rs relative to the QGSJET-II +  FLUKA model. For E γ = 10 TeV, the Rs in the leftmost bin (<10 muons) are less than 4.6 × 10 −3 , and for E γ = 100 TeV, the Rs in the leftmost bin (< 100 muons) are less than 3.1 × 10 −3 . The differences are also small. To estimate the gamma-ray flux, ALPAQUITA counts the number of excess (on-source minus off-source) events by real data. To convert the excess counts to gamma-ray flux, we need the gamma-ray detection efficiency calculated by a Monte Carlo simulation, where we expect that the gamma-ray detection efficiency will have some model dependence on hadronic interaction models. The next section will discuss the model dependence with a Monte Carlo simulation, assuming the ALPAQUITA configuration.

Cosmic gamma-ray flux estimation with ALPAQUITA
A conversion factor from the measured lateral distribution of air shower particles into incident particle energy, as well as the detection efficiency obtained by a Monte Carlo simulation of air showers with proper detector response, are required to estimate gamma-ray flux from an experimental data. For example, the differential flux ( where E γ is the energy of a gamma ray, N on is the number of events in an onsource window, N off is the number of events in an off-source window, T obs is the observation period, and S Eff is the effective area of the air shower array. R surv i (E γ ) is the survival ratio of gamma rays after the γ -CR separation using the number of measured muons. R surv i (E γ ) is calculated as Here, N rec sim (E γ ) is the number of reconstructed gamma-ray events, and N rec,γ -like sim (E γ ) denotes the number of events deemed gamma-like among N rec sim (E γ ). The index i = 1, 2, 3, and 4 represents QGSJET-II + FLUKA, QGSJET-II + UrQMD, EPOS LHC + FLUKA and SIBYLL + FLUKA, respectively. From (2) and (3), three factors S Eff , R surv i (E γ ), and E γ have potential hadronic interaction model dependence. As a gamma-ray-induced air shower is almost entirely composed of electromagnetic particles, the hadronic interaction model dependence on S Eff and E γ is regarded as negligible compared to that in R surv i (E γ ). Therefore, we will discuss only the hadronic interaction model dependence in R surv i (E γ ).

Air shower simulation and simulation setting
Assuming a gamma-ray source with a power-law energy spectrum with a spectral index of -2 in the direction of RX J1713.7-394, we generated 10 8 gamma-ray air showers in the energy range from 300 GeV to 10 PeV with CORSIKA 7.6400. Furthermore, we generated cosmic rays assuming the same source orbit and considered their isotropic characteristics by correcting the number of events using the weighting method described in a separate study [18]. As a chemical composition model of cosmic rays, we adopted the Shibata model [19], which was obtained by assuming a rigidity-dependent acceleration limit at the CR source. This model reproduces the results of direct observation experiments up to a few TeV and gradually becomes dominant by heavier nuclei from the region of tens of TeV. We produce 10 9 cosmicray events ranging from 300 GeV to 10 PeV. Table 2 summarizes the setup conditions. Regarding the hadronic interaction models in air shower generation, the four models explained in Section 2 are employed in this section. Each incident particle is randomly injected within a 300 m radius from the center of the air shower array. The particle information is traced to 1 MeV for secondary electrons and gamma rays, 50 MeV for muons, and 1 GeV for hadrons and is recorded at an altitude of 4,740 m.

Detector simulation of ALPAQUITA
Using a Monte Carlo simulation of ALPAQUITA [18] implemented in GEANT 4.10.02, we investigated the effect of the four hadronic interaction models described in Section 2 on gamma-ray efficiency. ALPAQUITA is a prototype detector for a new experiment ALPACA [16] [17] and is currently under construction on the Chacaltaya plateau (4,740 m a.s.l., 16 • .23 S, 68 • .08 W) in Bolivia. At 15 m intervals, 97 plastic scintillation detectors with an area of 1 m 2 are deployed, which is constituting an effective array area of 18,450 m 2 . A water Cherenkov-type muon detector (MD) is installed 2 m below the air shower array. The MD comprises 16 water tanks (total area of 900 m 2 ), and each tank has an air layer of 0.9 m, a water depth of 1.5 m, an area of about 56 m 2 . The concrete ceiling thickness is 20 cm, the thickness of concrete walls is 30 cm. The Cherenkov light photons are diffusely reflected by the walls and floor, which have an 80 % reflectance, and are detected by a downward-facing 20-inch-diameter PMT installed on the water tank's ceiling. The details of ALPAQUITA are described in a separate study [18],

General performance of the air shower array
The secondary particles generated with the CORSIKA are input to the GEANT4 simulation described in Section 3.2, and the behavior of each detector is simulated. The air shower event analysis method used in this study is the same as that described in [18]. We briefly summarize the general performance of ALPAQUITA below. The air shower array measures electrons, gamma rays, and charged particles in an air shower, their total energy loss is converted to the particle number density ρ (m −2 ). The primary energy is calculated as the sum of particle number densities from all detectors, ρ. Based on the relative hit timing of each detector, we reconstructed the arrival direction of the air shower using an air shower front surface, approximated with a cone [4]. The selection conditions described in [18] are required to analyze the air shower array data.

Separation method
The MD detects Cherenkov light photons reflected by the tank walls and floor by one PMT mounted downward on the ceiling [18]. The number of photoelectrons peaks at ∼24 when a single muon passes through the water tank. We define this value as one particle, calculate the number of muons N μ , and use it to select a gamma-ray-induced air shower. Figure 5a shows a ρ vs. N μ scatter plot using the QGSJET-II + FLUKA model, with cosmic rays distributed in the upper region and the gamma rays distributed in the lower region. Figure 5b, c, and d show the same trend as the QGSJET-II + FLUKA model, with the QGSJET-II + UrQMD model, EPOS LHC + FLUKA model, and SIBYLL + FLUKA model, respectively. Figure 6 shows the number of muon distributions for 56.2 < ρ < 100, corresponding to a representative energy of ∼28.8 TeV. As shown in Fig. 6, the N μ of a gamma-ray-induced air shower is significantly smaller than that of a cosmic-ray-induced air shower with the same ρ. Events with N μ < 0.1 are artificially piled up at N μ = 0.01, where majority gamma rays are contained. For the QGSJET-II + FLUKA model, the number determined by the method to discriminate between gamma-ray and cosmic-ray events, described in [18]. The three survival lines in each figure are the same as determined by the QGSJET-II + FLUKA model of gamma-ray events in this bin was 21.8, corresponding to 62.4 % of the total number, whereas the number for the other three models ranged from 21.3 to 21.9. We define gamma-like events below the threshold and cosmic-ray-like events above the threshold using a threshold value of N μ ( N μ,cut ). The portion of total gamma-ray events that remained below the threshold (survival ratio) is 90 % when N μ,cut is set to 2.57 in Fig. 6a. For each ρ bin, we determine different thresholds to satisfy survival ratios of N μ,cut s satisfying the survival ratio at 50 % and 90 %.

Optimal separation
When N μ,cut is high, we can maintain a high detection efficiency for gamma rays, but rejecting cosmic-ray events is challenging. Conversely, when N μ,cut is low, the rejection power of cosmic-ray events is high, but the survival ratio of the gamma-ray events is small. Therefore, there is an optimal N μ,cut value. We used the quality factor (Q factor = N γ / N B + N γ ) [18] to determine the optimal threshold value, where N γ is the number of gamma-ray events and N B is the number of cosmic-ray events. The optimal value of N μ,cut for the QGSJET-II + FLUKA model in Fig. 6a is 0.63 at 28.8 TeV. For each ρ bin, we determine the optimal N μ,cut and obtain the relation between ρ and N μ,cut . We fitted the relationship to the following equation: In the case of QGSJET-II + FLUKA, the optimal line is the optimal survival line shown in Fig. 5, where b becomes 1.6 with a fixed ρ 0 = 31.2 [18]. Utilizing the optimal survival line shown in Fig. 6, we obtained a survival ratio of ∼0.7 for all models at the representative energy of ∼28.8 TeV.
We employ a survival line similar to the optimal survival line in the realistic case. Therefore, we estimate the final hadronic interaction model dependence, assuming the optimal survival line. The final hadronic interaction model dependence includes conservatively differences among the models in the high-energy region [the average fitting result R(R surv i=4 ) = 5.5 × 10 −3 ± 1.2 × 10 −2 MCstat and bin-by-bin maximum difference (R(R surv i=4 )) max = 2.2 × 10 −2 ± 1.1 × 10 −2 MCstat ] and differences among the models in the low-energy region [the average fitting result R(R surv i=2 ) = 2.0 × 10 −3 ± 9.5 × 10 −3 MCstat and bin-by-bin maximum difference (R(R surv i=2 )) max = 1.8 × 10 −2 ± 1.0 × 10 −2 MCstat ]. We then estimate the final hadronic interaction model dependence, adding the differences above in quadrature, We obtained the final hadronic interaction model dependence of less than 3.6 × 10 −2 , which is dominated by Monte Carlo statistics.

Conclusion
Using an extensive air shower array and an underground MD to measure high-energy gamma ray, we studied the performance of a hybrid experiment. The experiment separates gamma-ray-induced muon-poor air showers from cosmic-ray-induced muon-rich air showers using an underground muon detector. To investigate the performance of this type of experiment, such as ALPAQUITA, a Monte Carlo simulation of air showers based on hadronic interaction models is required, which might have some model dependence. Air shower experiments such as ALPAQUITA estimate the gamma-ray flux from the difference between on-source and off-source events by real data, using the gamma-ray detection efficiency calculated by a Monte Carlo simulation, which depends on the hadronic interaction models, whereas the off-source, background cosmic-ray events, can be estimated experimentally. In particular, the models affect the number of muons in the gamma-ray-induced air showers.
To evaluate the differences in the characteristics of air shower muons at 4,740 m above sea level, we simulate the gamma-ray-induced air showers with the four models-QGSJET-II + FLUKA, QGSJET-II + UrQMD, EPOS LHC + FLUKA, and SIBYLL + FLUKA. Thus, we evaluate the hadronic interaction model dependence in the gamma-ray flux measurements between 6 and 300 TeV.
First, before including detector simulation, we studied the hadronic interaction model dependence of gamma-ray-induced air showers at the CORSIKA level. In the 10 and 100 TeV gamma-ray-induced air showers over 85 % of muons have energies below 10 GeV, and model differences in the integral muon spectrum in several GeV regions are small, less than a few percent. For the lateral spread of muons, the differences were also small. They are ± 0.77 % for 10 TeV gamma-ray air showers and ± 1.8 % for 100 TeV gamma-ray air showers when we compare the total number of muons above 1 GeV within 100 m (∼Molière unit) in radius from the air shower core in the lateral distribution. Approximately 90% of 10 TeV gamma-ray-induced air showers contain less than 10 muons, whereas about 90 % of 100 TeV gammaray-induced air showers contain less than 100 muons. For both 10 TeV air showers with less than 10 muons and 100 TeV air showers with less than 100 muons, model differences in the total number of gamma-ray-induced air showers are less than 1 %.
Using the ALPAQUITA simulated detector, we evaluated the impact of these small differences in the characteristics of air shower muons on the gamma-ray flux estimation. For optimal survival, the survival ratio of gamma-ray events varied from 0.7 to 0.9 for ρ in the energy range from 6 to 300 TeV. Thus, the survival ratio's model dependence is less than 3.6 % in the energy range from 6 to 300 TeV. The survival ratio is included in the flux calculation in (2) in the form of (R surv i ) −1 ; thus, the model dependence in the flux estimation is of the same magnitude, which is dominated by Monte Carlo statistics. The contribution of the hadronic interaction model dependence to the gamma-ray flux estimation is negligible compared to other systematic errors, such as the energy scale uncertainty (typically ∼10 %) corresponding to the 20 %-30 % gamma-ray flux uncertainty. Furthermore, as Monte Carlo statistics account for most of the 3.6 % uncertainty due to hadronic interaction model dependence, it is projected to be significantly less than 3.6 %.