# No muon excess in extensive air showers at 100–500 PeV primary energy: EAS–MSU results

###### Abstract

Some discrepancies have been reported between observed and simulated muon content of extensive air showers: the number of observed muons exceeded the expectations in HiRes-MIA, Yakutsk and Pierre Auger Observatory data. Here, we analyze the data of the Moscow State University Extensive Air Shower (EAS–MSU) array on GeV muons in showers caused by eV primary particles and demonstrate that they agree with simulations (QGSJET-II-04 hadronic interaction model) once the primary composition inferred from the surface-detector data is assumed.

###### keywords:

extensive air showers, muon data, hadronic interaction models^{†}

^{†}journal: Astroparticle Physics

## 1 Introduction

Ultra-high-energy cosmic rays provide a unique laboratory to study hadronic interactions at the center-of-mass energies and in kinematical regimes not accessible at colliders. Modelling of the development of an extensive air shower (EAS), a cascade process in the terrestrial atmosphere initiated by an energetic cosmic particle, requires an extrapolation of verified interaction models. Not surprisingly, this often results in discrepancies between measured and simulated EAS properties, or between physical properties of the primary particle reconstructed by different methods. A well-known result, possibly related to the lack of understanding of the EAS development New-PAO (), is the systematic difference between the primary energies reconstructed by the fluorescence detectors and by surface arrays for very same events, as seen by the Pierre Auger Observatory (PAO) Engel () and the Telescope Array (TA) experiment TA-E (). It may or may not be related to the apparent excess of muons (GeV) in EAS reported at eV by the PAO New-PAO (); Engel (); PAOmu2 () and Yakutsk Yak-mu () experiments. A similar excess had been observed earlier by the HiRes/MIA experiment at eV HiRes-MIA (). The purpose of the present study is to compare observed and simulated densities of GeV muons in air showers induced by eV primaries, based on the EAS-MSU data.

A subtle point of all comparisons of this kind is that the muon content of a EAS depends strongly on the type of a primary particle. As a result, the average muon content in the MC set depends not only on the hadronic interaction model used, but also on the primary composition assumed at the simulation. Therefore, for a meaningful comparison, one needs an independent estimator of the primary composition in the very same data set for which the muon data are analysed. An estimator of this kind is often missing. In this work, we take advantage of the knowledge of the primary composition obtained from the surface-detector data only, as discussed below.

The rest of the paper is organized as follows. In Sec. 2, a brief description of the installation and of the data set is given, together with references to previous more detailed publications about the EAS-MSU array. In Sec. 3, we discuss the analysis performed in this work. Section 4 presents our results, while Sec. 5 summarizes our conclusions.

## 2 Data

### Installation

The EAS-MSU array EAS-MSU (), located in Moscow (190 m a.s.l., corresponding to the atmospheric depth g/cm), operated in 1955–1990 in various configurations. A detailed description of the array in the ultimate configuration, whose data are discussed here, may be found in Ref. JINST (). The total area of the array, 0.5 km, was covered with 76 charged-particle detector stations, each consisting of multiple Geiger–Mueller counters. A unique feature of the installation was the presence of large-area muon detectors. The main muon detector, located in the array center, had a total area of 36.4 m and consisted of 1104 similar Geiger-Mueller counters, each with an area of 0.033 m, located at a depth of 40 meters of water equivalent underground. It is the data of this detector what we study here.

The geometry of the array, the trigger system and the reconstruction procedure are described in detail in Ref. JINST (). The surface detector stations allow to reconstruct the lateral distribution function (LDF) for charged particles, parametrized as

(1) |

where is the particle density at the detector plane (we study vertical showers, , and thus neglect the asymmetry caused by attenuation effects), is the distance to the shower axis, m is the Moliere radius, is the modified age parameter of the EAS, is the normalization coefficient, which is calculated numerically, and is the correction to determined empirically and presented e.g. in Ref. JINST (), ; the normalization determines the total number of charged particles.

### Basic selection cuts

The following cuts were imposed for the high-energy surface-detector data sample studied in this work:

1. Convergence of the reconstruction.

2. The LDF age parameter belongs to the interval .

3. The reconstructed zenith angle satisfies .

4. The reconstructed shower axis is within 240 m from the array center, where the muon detector is located.

5. The reconstructed shower size is , which corresponds to eV, assuming protons as primaries .

For this study, we use the data recorded in the period 1984 — 1990 (1372 days). After application of the cuts, the surface-detector data sample contains 922 events. This sample is representative in the sense that the event set agrees with the Monte-Carlo simulations based on the EAS-MSU spectrum, as discussed in Ref. JINST (). The efficiency of the installation is greater than 95% for primary energies eV for these cuts.

### Muon detectors

As it has been mentioned above, the muon detectors are shielded by a thick
layer of soil which effectively absorbs all the shower particles but
muons: the fraction of other particles in the detector is far below the
level. The threshold energy of vertical muons registered
by the underground detectors JINST ()
is^{1}^{1}1As it has been discussed in
Ref. JINST (), the detector is shielded by the soil of
40 m.w.e. The energy-dependent muon absorption is treated in the
continuous-slowing-down approximation CSDA (). For zenith
angles studied here (), the threshold
energy scales obviously with . GeV.
The main muon detector consisted of 32 independent sections. The muon
density is estimated as , where
is the number of counters in the muon detector, is the number of
fired counters and is the area of the counter^{2}^{2}2This
formula follows from the binomial distribution. Indeed, the probability
that one counter does not fire is estimated as . The probability to have counters fired is
where .
It is maximized for , which results in .. In our analysis, we further
restrict the data sample to the days when not less than 28 of 32 sections
were operational and no failure of the muon detector was recorded in
24 hours. This removes 168 days of data, so the final
sample for muon studies contains 809 events recorded
during the live time of 14060 hours, so that the overall exposure is km s sr.

## 3 Analysis

The muon density in an air shower decreases with the distance from the axis, , so to compare this quantity between data and MC, one often uses the LDF to recalculate the density to a particular fixed value of . In this study, we use , the muon surface density recalculated to the typical core-detector distance of m with the help of the EAS-MSU muon LDF determined in previous works EAS-MSU-muon-LDF (),

(2) |

where is fixed as in Eq. (1) and . This muon LDF agrees reasonably well with the data as well as with MC simulations. Figure 1

compares the LDF used in our study with simulations performed for proton and iron primaries. Allowing for variations of the parameter in Eq. (2), we obtain for primary protons, for primary iron and for the best-fit mixture (57% Fe, 43% p) describing the surface-detector data, see below. The same fitting procedure for the data gives . Data and simulations agree well with the value fixed in the analysis.

The limited number of events in the data set and the limited range used make the experimental study of the possible dependence in the muon LDF hardly possible. In Fig. 2,

we compare Monte-Carlo simulations (QGSJET-II-04) for vertical () and the most inclined () events. Separate fits for these ranges result in (vertical) and (inclined), which supports the use of -independent muon LDF in our study.

Comparison of the distributions in for the real and
simulated data constitutes the main part of this work. To produce a
reliable simulation, we make use of the full Monte-Carlo (MC) model of the
installation. This model, developed and described in detail in
Ref. JINST (), accounts for the air-shower production
and development in the atmosphere, its detection and reconstruction.
The air-shower simulations have been performed with the CORSIKA 7.400.1
code CORSIKA (), using QGSJET-II-04 QGSJET-II-04 () as the
high-energy hadronic interaction model, FLUKA 2011.2c FLUKA () as the
low-energy hadronic interaction model and EGS-4 EGS-4 () as the
electromagnetic-interaction model. No thinning was employed; the Central
European atmosphere for October 14, 1993 (CORSIKA model number 7) was
used^{3}^{3}3The mean atmospheric pressure in Moscow
for 1984–1991, recalculated to the sea level, is hPa yak (). Correction to the installation altitude of
190 m gives hPa. The CORSIKA model 7 gives
1020.1 hPa for the sea level and 997.7 hPa for 190 m.. The detector
response to individual particles was modelled as described in
Ref. JINST (). We used 1370 independent CORSIKA showers (852 protons and 518 iron
nuclei), from which we derived, by random moving of the shower axis, our
MC sample (see Ref. JINST () for more details). After applying all
cuts, the MC sample contains 4468 proton-induced and 1093 iron-induced
showers.

The artificial events are recorded in the same data format as the real ones and are processed with the same reconstruction software. This means that all selection effects related to the trigger, threshold and efficiency are accounteed for in the simulation in the same way as in the real data. The resulting distributions of the reconstructed surface-detector parameters agree well between data and MC JINST (). The estimated uncertainty (68% containment radius around the mean) in the reconstruction of the axis position is 5.7 m; for the arrival direction it is 1.1 and for it is 16.5%.

The muon density in air showers is a composition-dependent observable: heavier primary nuclei produce more muons. Many surface-detector observables are degenerate with respect to the primary composition, and a good description of these data might, in principle, be achieved for various assumptions about the composition. That would lead to quite different predictions for the muon content, however. To perform a detailed comparison of muon data with simulations, one therefore requires the knowledge of the primary composition fully independent on muon detector readings. This is not always easy to achieve, and even in modern experiments, surface-detector composition studies are complicated and not very precise, see e.g. Refs. PAO-SD-comp (); TA-SD-comp (). Fortunately, the EAS-MSU setup provided for a required observable. Its surface array was very dense in its central part, and the slope of the LDF, parametrized by the age parameter in Eq. (1), is determined with a precision sufficient to determine the average primary composition with a reasonable accuracy JINST (). Indeed, an event in the sample we consider has, on average, 15 detectors used for the LDF reconstruction, which guarantees the required accuracy. We note that, in principle,the LDF slope may also, like the muon number, demonstrate disagreement between data and simulations with present-day interaction models. This fact, however, is not very relevant for the present study, since most of the events lay on the distance not exceeding few hundred meters from the center of the array. The MC and measured age parameters are in agreement in this radial distance range.

In Ref. JINST (), working in the frameworks of the two-component mixture of primary particles (protons and iron), we determined, for , the best-fit composition which describes the surface-detector data, including the distribution. Within the statistical uncertainties, the composition does not change with energy in the narrow range of energies we study. We use this mixture (43% protons and 57% iron, following the common power-law spectrum ) in MC simulations and obtain the expected distribution in which we compare to that observed in the real data. Then, in order to quantify potential muon excess in data over simulations, we introduce the coefficient by which the muon number is scaled in simulated showers. By definition, corresponds to the muon number predicted in our simulations with the QGSJET-II-04 hadronic interaction model and with the determined above primary composition. By means of the binned chi-squared test, we compare distributions derived for various with the data and determine, consequently, the allowed range of . The scaling of the muon number was implemented only for muon density measured by underground detectors, but not for the surface-detector observables. We discuss the justification for this assumption in Appendix.

## 4 Results

Figure 3

compares the distributions in obtained from the simulations described above () and from the data. One may see that the distributions are in a good agreement.

To further quantify this agreement, one may proceed in two ways. Firstly, note that the muon content depends strongly on the assumed primary composition, see Fig. 4(a).

We find that the proton/iron mixture which describes the distribution is iron, see Fig. 4(b) for the best binned-likelihood fit, which agrees with 57% iron determined from the surface-detector data.

Secondly, we study how the scaling of the muon number, that is the variation of , affects the agreement between data and simulations for the distribution. This comparison is illustrated in Fig. 5,

where the normalized chi-squared is presented as a function of . The standard statistical analysis results in , so that no muon excess in data is observed and is excluded by the data at the 92% CL.

All statistical uncertainties (related to air-shower fluctuations and fluctuations in the detector response, as well as uncertainties in the reconstruction procedure) together with potential reconstruction biases are taken into account by the procedure we use in which data and MC are processed in the same way. Intrinsic fluctuations of the muon density are taken into account because we use no-thinning simulations. The uncertainty in the determination of the muon density is taken into account in simulations, since the probability of the counter to fire being hit by a muon is close to . However, there are potential sources of systematic uncertainties which we are turning to now.

### Muon LDF

The systematic uncertainty, introduced by the variations of to the Fe fraction determined from , does not exceed .

### Composition

The main result of the paper is based on the primary composition derived from surface-detector data for the same data set. We tested its stability when relaxing this assumption. Assuming the p/Fe mixture reproducing the results of KASCADE-Grande KASCADE-comp () (Fig. 7 of Ref. KASCADE-comp (), stars; 59% iron), we obtain ; for Tunka-133 Tunka-comp () (51% iron) it is . The confidence levels of exclusion of are 95% and 67%, respectively.

One may note that, suggested by results of other experiments, the composition at eV changes due to the “heavy knee” in the spectrum. However, working within a rather narrow energy interval, we do not have sufficient statistics to trace this change in the data. We also do not go beyond the two-component proton-iron mixture. Previous studies have demonstrated JPG () that the EAS-MSU experiment does not have enough statistics to distinguish the contribution of medium-weight nuclei. On the other hand working with the p-Fe mixture and constant composition is supported by the data/MC agreement.

### Interaction models

The main conclusion of the present study, the absence of an excess of muons in data over MC, was obtained for one particular hadronic-interaction model, QGSJET-II-04. To estimate the effect of the change of the interaction model on the result, we performed a simplified study with three other models, namely EPOS-LHC EPOS-LHC (), SIBYLL 2.1 SIBYLL-2.1 () and QGSJET-01 QGSJET-01 (). The aim of the study was to estimate the model-related systematic uncertainty of our conclusions. The simulations were performed with the thinning parameter and maximal weight restrictions of 40 for hadrons and 4000 for electromagnetic particles, as suggested by the method of “optimized thinning”, Ref. Kobal (). For each model, we simulated 400 showers with eV and 400 showers with eV for primary protons and the same amount of showers for primary iron. The age parameter was determined by averaging of particle densities in concentrical rings around the shower axis, the muon content of the shower was studied using the total number of muons, which is uniquely translated to muon density by the muon LDF.

It has

Model | mean Fe fraction | |

from , % | ||

QGSJET-II-04 | 57 | 0.920.06 |

EPOS-LHC | 42 | 0.960.06 |

SIBYLL 2.1 | 76 | 1.000.07 |

QGSJET-01 | 58 | 0.950.06 |

been shown by an explicit comparison of thinned and non-thinned showers that the procedure does not change the central values of LDF-related parameters, though it does introduce additional fluctuations livni (). We re-checked that the central values remain the same as in the full simulation by performing the same simplified simulations with our principal model, QGSJET-II-04. Assuming the measured spectrum, we estimated the iron fraction from , repeating the procedure for this case. We determined the coefficient in the same way as in the main part of the study. The results are presented in Table 1.

We see that though particular best-fit composition differs from model to model (by in terms of the iron fraction in the iron-proton mixture), the systematic uncertainty in , related to the interaction model, does not exceed . The muon excess is absent in all cases, in agreement with the observation that the fit never requires primaries heavier than iron.

## 5 Conclusions

We have analyzed the densities of muons ( GeV) registered by underground detectors of the EAS-MSU experiment. Starting from the Monte-Carlo simulation based on the primary composition inferred from the surface-detector data alone and on the QGSJET-II-04 hadronic interaction model, we obtain a good agreement between the simulations and the data. Assuming that the number of muons in air showers scales with a coefficient with respect to the simulation, we constrain , so that no muon excess () is observed and agrees with the data at the confidence level. Similar conclusions are obtained for primary composition assumptions favoured by the results of other experiments.

We note that the nice agreement between predicted and observed muon densities reported here does not necessarily mean that QGSJET-II-04 gives a correct description of the muon production in any case. The agreement observed here relates to eV, GeV and inner parts of the shower, . Previous results, collected in Table 2,

Experiment | altitude, | , | , eV | , | muon excess | ||

m a.s.l. | g/cm | GeV | (data over MC) | ||||

HiRes-MIA HiRes-MIA () | 1500 | 860 | N/A | yes | |||

PAO Engel (); PAOmu2 () | 1450 | 880 | 70 | yes | |||

Yakutsk Yak-mu () | 100 | 1020 | 45 | yes | |||

IceTop IceTop () | 2835 | 680 | 13 mean | no | |||

EAS-MSU (this work) | 190 | 990 | 30 | no |

have been obtained in various different regimes, and some at different altitudes. The muon excess reported by PAO Refs. Engel (); PAOmu2 () and Yakutsk Yak-mu () was observed at primary energies eV and muon energies GeV. HiRes-MIA HiRes-MIA () observed the excess for GeV at eV eV. Contrary, recent preliminary IceTop results IceTop () for GeV muons and eV eV suggest that no excess is seen. One should not forget also the important difference between our work and all these studies: here we investigate the inner parts of EAS, , while the results of other experiments refer to the outer parts, . Note that at even lower and higher TeV, the muon excess may be probed with the help of atmospheric muons DedJETPL (), and it has been reported in 1504.05853 () that QGSJET-II-04 seems to overestimate the number of muons in this regime. Preliminary results of the KASCADE-Grande experiment (110 m a.s.l., 1022 g/cm) at eV suggest KASCADE-muon-att () that the atmospheric attenuation of the muon number ( GeV, ) is underestimated by all high-energy hadronic interaction models studied there, including QGSJET-II-04. Clearly, further experimental and theoretical studies are required to understand the origin of the reported discrepancies and to arrive at a succesful model of the air-shower development.

## Appendix. Muon number scaling and the surface detectors

The scaling of the muon number in our Monte-Carlo simulations was implemented for underground detectors, but not for the surface-detector observables. This assumption was tested by the analysis of surface-detector observables in the Monte-Carlo showers with downscaled number of muons. The downscaling itself is a random removing of a fraction of muons from the CORSIKA showers. Comparing the showers with and , we found the mean change in the reconstructed of 2.5%. The root mean squares of the variations of the principal observables (primed quantities correspond to ) are: , , m, . We should note that the downscaling of muon number affects the SD observables stronger than its upscaling, so the presented values of variations can be considered as upper limits on what we can get with the upscaling of the muon number. At the same time, these values are smaller than the experimental uncertainties, see Ref. JINST (), therefore we can neglect the impact of the muon-number scaling on surface-detector observables.

## Acknowledgements

ST acknowledges interesting and stimulating discussions of the “muon excess” with Ralf Engel, Sergey Ostapchenko and Gordon Thomson. Monte-Carlo simulations have been performed at the computer cluster of the Theoretical Physics Department, Institute for Nuclear Research of the Russian Academy of Sciences. The experimental work of the EAS-MSU group was funded by the Government of the Russian Federation (agreement 14.B25.31.0010) and the Russian Foundation for Basic Research (project 14-02-00372). Development of the analysis methods and application of them to the EAS-MSU data were funded by the Russian Science Foundation (grant 14-12-01340).

## References

- (1) A. Aab et al. [Pierre Auger Collaboration], Phys. Rev. Lett. 117 (2016) 192001 [arXiv:1610.08509 [hep-ex]].
- (2) R. Engel [Pierre Auger Collaboration], arXiv:0706.1921 [astro-ph].
- (3) T. Abu-Zayyad et al. [Telescope Array Collaboration], Astrophys. J. 768 (2013) L1 [arXiv:1205.5067 [astro-ph.HE]].
- (4) A. Aab et al. [Pierre Auger Collaboration], Phys. Rev. D 91 (2015) 032003.Erratum: [Phys. Rev. D 91 (2015) 059901] [arXiv:1408.1421 [astro-ph.HE]].
- (5) A. V. Glushkov et al. JETP Lett. 87 (2008) 190 [arXiv:0710.5508 [astro-ph]].
- (6) T. Abu-Zayyad et al. [HiRes and MIA Collaborations], Phys. Rev. Lett. 84 (2000) 4276 [astro-ph/9911144].
- (7) S. N. Vernov et al., Bull. Russ. Acad. Sci. Phys. 44 (1980) 80 [Izv. Ross. Akad. Nauk Ser. Fiz. 44 (1980) 537].
- (8) Y. A. Fomin et al., JINST 11 (2016) T08005 [arXiv:1607.00309 [astro-ph.HE]].
- (9) D. E. Groom, N. V. Mokhov and S. I. Striganov, Atom. Data Nucl. Data Tabl. 78 (2001) 183.
- (10) G.B. Khristiansen, G.V. Kulikov, Yu.A. Fomin. Cosmic Rays of Superhigh Energies, Karl Thiemig Verlag, Munchen, 1980, 246 p.
- (11) D. Heck et al., CORSIKA: A Monte Carlo code to simulate extensive air showers, Tech. Rep. 6019, FZKA (1998)
- (12) S. Ostapchenko, Monte Carlo treatment of hadronic interactions in enhanced Pomeron scheme: I. QGSJET-II model, Phys. Rev. D 83 (2011) 014018 [arXiv:1010.1869 [hep-ph]].
- (13) G. Battistoni et al., The FLUKA code: Description and benchmarking, AIP Conf. Proc. 896 (2007) 31.
- (14) W. R. Nelson, H. Hirayama, D. W. O. Rogers, The EGS4 code system, Tech. Rep. 0265, SLAC (1985).
- (15) L. Kartsev, A. Sergeyev, V. Petrov, A. Gavrishev, Weather statistics for Russian cities, http://www.atlas-yakutia.ru (in Russian).
- (16) P. Abreu et al. [Pierre Auger Collaboration], arXiv:1107.4804 [astro-ph.HE].
- (17) G. I. Rubtsov et al. [Telescope Array Collaboration], J. Phys. Conf. Ser. 608 (2015) 012067.
- (18) P. Luczak et al., EPJ Web Conf. 99 (2015) 13001.
- (19) S. F. Berezhnev et al., Proc. 32nd ICRC, Beijing 1 (2011) 209.
- (20) Yu. A. Fomin et al., J. Phys. G 22 (1996) 1839.
- (21) T. Pierog et al., Phys. Rev. C 92 (2015) 034906 [arXiv:1306.0121 [hep-ph]].
- (22) E. J. Ahn et al., Phys. Rev. D 80 (2009) 094003 [arXiv:0906.4113 [hep-ph]].
- (23) N. N. Kalmykov and S. S. Ostapchenko, Phys. Atom. Nucl. 56 (1993) 346 [Yad. Fiz. 56 (1993) 105].
- (24) M. Kobal [Pierre Auger Collaboration], Astropart. Phys. 15 (2001) 259.
- (25) D. S. Gorbunov, G. I. Rubtsov and S. V. Troitsky, Phys. Rev. D 76 (2007) 043004 [astro-ph/0703546 [ASTRO-PH]].
- (26) J. G. Gonzalez [IceCube Collaboration], J. Phys. Conf. Ser. 718 (2016) 052017.
- (27) L. G. Dedenko, T. M. Roganova and G. F. Fedorova, JETP Lett. 100 (2014) 223 [Pisma Zh. Eksp. Teor. Fiz. 100 (2014) 247].
- (28) L. G. Dedenko, A. V. Lukyashin, G. F. Fedorova and T. M. Roganova, EPJ Web Conf. 99 (2015) 10003 [arXiv:1504.05853 [astro-ph.HE]].
- (29) J. C. Arteaga-Velazquez et al., EPJ Web Conf. 99 (2015) 12002.