High-energy cosmic rays interacting in the Earth's atmosphere produce large cascades of secondary particles called extensive air showers, which are observed with large ground-based detector arrays. To determine the properties of the primary particle, the observations are compared to detailed simulations of the air-shower development in the atmosphere and the detector response. This indirect method leads to large uncertainties resulting from the incomplete description of the high-energy hadronic interactions in the atmosphere, as demonstrated by the various discrepancies between observations and predictions reported in the recent past [1]. This is in particular the case for the muonic component of air showers, which originates predominantly from the decay of charged pions and kaons produced in the hadronic cascade.
The IceCube Neutrino Observatory [2] can provide unique input to this problem by observing air-shower muons in different energy regimes. Its surface detector, IceTop [3], is sensitive to the electromagnetic component as well as the low-energy (O(GeV)) muon component of the shower. The in-ice detector can observe muons with energies above several 100 GeV from the same air showers. In this contribution, we will present results from two recent IceCube analyses of the lowand high-energy muon content in near-vertical air showers [4,5]. The IceCube Neutrino Observatory (IceCube) is a kilometerscale particle detector located in Antarctica [2]. It instruments a cubic-kilometer of ice with about 5000 Digital Optical Modules (DOMs), capturing the Cherenkov light produced by relativistic charged particles at depths of about 1.5 km to 2.5 km below the surface. The DOMs are deployed on 86 vertical strings forming a hexagonal grid with a horizontal spacing of about 125 m. In addition, IceCube includes a surface air-shower array, IceTop, instrumenting an area of about 1 km 2 with cylindrical ice-Cherenkov tanks containing two DOMs each [3]. The tanks are deployed in pairs, on approximately the same grid as the in-ice strings.
IceTop is located at an elevation of about 2.8 km above sea level, corresponding to an average atmospheric depth of about 690 g cm -2 , where it observes cosmic-ray air showers in the PeV to EeV primary energy range. Due to the high altitude, IceTop signals are dominated by the electromagnetic shower component, with the signal produced by low-energy muons becoming visible only at large lateral distances. Muons with energies above several 100 GeV in the air shower can penetrate all the way to the IceCube in-ice array, producing the so-called coincident events illustrated in Fig. 1.
Muons in air showers with IceCube
Air showers observed with IceTop are reconstructed with a typical likelihood method fitting the lateral charge and time distributions, providing the core position and arrival direction [3]. The reconstruction also provides the shower size 𝑆 125 , the signal strength at a reference distance of 125 m. The IceTop signals are calibrated in units of Vertical Equivalent Muon (VEM), the typical signal produced by a muon vertically traversing a tank. The reconstructed 𝑆 125 is strongly correlated with the primary cosmic-ray energy [6] and is used as an energy estimator in both muon analyses described below.
Both analyses are limited to near-vertical showers with a reconstructed zenith of cos 𝜃 > 0.95 or 𝜃 18 • . They include only air showers whose core is contained within the boundary of the IceTop array, and for the coincident events, whose shower axis goes through the in-ice array so that a high-energy muon bundle is observed. For more details about the events selection, see Refs. [4] and [7]. The density of muons at large lateral distance at the surface has been derived in an analysis using only the IceTop detector [4]. The density of these mainly low-energy (O(GeV)) muons was determined at distances of 600 m and 800 m from the shower axis. Predictions for these quantities obtained from COR-SIKA [8] simulations of proton and iron showers are shown in Fig. 2. The high-energy hadronic interaction models included are Sibyll 2.1 [9], QGSJet-II.04 [10], and EPOS-LHC [11], and the simulations were performed using an atmosphere approximately describing the yearly average South Pole atmosphere. The post-LHC models QGSJet-II.04 and EPOS-LHC predict about 20%-35% more muons than the older Sibyll 2.1.
While IceTop signals near the shower core are dominated by the electromagnetic shower component, at large distances from the core the contribution from muons becomes significant. This has been exploited in a statistical analysis in which a large number of events are combined in bins of 𝑆 125 . In each bin, a two-dimensional histogram of observed charges versus the lateral distance of the tanks is created, as shown in Fig. 3. At large distances, the muon component can be observed as a distinct population of signals around 1 VEM. To determine the average muon density in the sample as a function of distance, vertical slices in the histograms are fit with a multi-component signal model that includes an electromagnetic component, a muon component, as well as uncorrelated background signals. Based on these fits, the muon density values are determined at radial distances of 600 m and 800 m. To account for small differences between reconstructed and true muon densities observed in simulations, corrections are derived from the simulations and applied to the data.
The results derived using correction factors based on the models Sibyll 2.1, QGSJet-II.04, and EPOS-LHC are shown in Fig. 4. The muon densities at 600 m are derived for primary energies
Muons in air showers with IceCube between 2.5 PeV and 40 PeV, those at 800 m for energies between 9 PeV and 120 PeV. The results are compared to the predictions from simulated proton and iron showers using the corresponding hadronic interaction model. Systematic uncertainties related to the derivation of the correction factors, the unknown mass composition, the energy scale, and the electromagnetic signal model in the likelihood fit are taken into account. The difference between the results obtained with different hadronic models is about 15% and below; the lighter composition implied in the QGSJet-II.04 and EPOS-LHC plots results mainly from the increase in the predicted muon densities, as seen in Fig. 2.
IceCube has also determined the average number of muons with energies above 500 GeV in near-vertical showers [5]. This analysis uses coincident events, i.e. those detected with IceTop with a corresponding high-energy muon bundle going through the in-ice array. A large fraction of muons above this energy threshold can make it to the detector, and a neural network is used to relate the signal observed in the in-ice detector to the muon number at the surface. Fig. 5 shows the predictions from CORSIKA simulations for this observable, simulated using the same conditions
Muons in air showers with IceCube Stef Verpoest as those described in Section 3. QGSJet-II.04 predicts slightly more high-energy muons than Sibyll 2.1, while EPOS-LHC predicts slightly less. The analysis utilizes a reconstruction of the energy loss profile of the muon bundle in the in-ice detector. The deposited energy is unfolded in track segments of 20 m length [12], along a seed track which is given by the shower axis as reconstructed by IceTop. The reconstructed energy losses are then fed into a recurrent neural network layer. The output of the recurrent neural network is combined into a fully-connected network layer together with the shower size 𝑆 125 and the zenith angle 𝜃 from the IceTop reconstruction. This layer finally outputs an estimate of the primary cosmic-ray energy and of the number of muons with 𝐸 𝜇 > 500 GeV in the shower at the surface. The relation between true and reconstructed muon number is shown in Fig. 6.
The event-by-event estimates are used to obtain the average muon number 𝑁 𝜇 as a function of the cosmic-ray energy. A small Monte-Carlo based correction is applied to the obtained 𝑁 𝜇 to take into account systematic biases resulting from the biases the muon number reconstruction as well as from energy bin migration. Several checks have been performed that verified the robustness of the results under different approaches for the reconstructions, such as using separate neural networks using only IceTop information for the energy reconstruction and only in-ice information for the muon number reconstruction. Results are derived assuming three hadronic interaction models for the correction. The high-energy muon numbers derived from data are shown in Fig. 7, compared to predictions from simulated proton and iron showers using the corresponding hadronic interaction model. The systematic uncertainty is dominated by uncertainties in the simulated ice model and the DOM efficiency, and has smaller contributions related to the snow accumulation on IceTop, the IceTop charge calibration, the simulated atmosphere, and the Monte-Carlo based correction factors. The results for all three models are in between the proton and iron predictions; the difference between the results derived using different models is only about 8% and below.
Muons in air showers with IceCube
Figure 7: High-energy (𝐸 𝜇 > 500 GeV) muon number measured in coincident events assuming different hadronic interaction models. Shown for comparison are the predictions for proton and iron from the corresponding model, as shown earlier in Fig. 5. The total systematic uncertainty is shown by the shaded area, statistical uncertainties are too small to be visible.
A common way of comparing muon measurements to predictions from simulations is through the definition of the so-called z-value [13],
where 𝑁 𝜇 is the result obtained from data, and 𝑁 𝜇 p and 𝑁 𝜇 Fe are the predictions from simulated proton and iron showers respectively (and equivalent for the muon density 𝜌 𝜇 ). Both the muon density and high-energy muon multiplicity results are shown in this representation in Fig. 8, where they are also compared to predictions based on different cosmic-ray flux models. The muon density results for Sibyll 2.1 agree well with the expectations from these models, while QGSJet-II.04 and EPOS-LHC imply a lighter composition. For the high-energy muons, all results agree within uncertainties with the predictions, with the EPOS-LHC result indicating only a slightly heavier composition than that inferred from the Sibyll 2.1 and QGSJet-II.04 results.
If the air-shower simulations give an accurate description of reality, the z-values for the muon density and the high-energy muon number results should be consistent with each other in their overlapping energy range of 2.5 PeV to 100 PeV. While good agreement is found for the results based on Sibyll 2.1, a tension is observed for the results based on QGSJet-II.04 and EPOS-LHC, resulting from the increased production of low-energy muons in these models (Fig. 2). This finding is in line with an earlier preliminary IceCube result [14], which additionally has indications for an inconsistency between the muon results and the slope of the lateral charge distribution in IceTop for Sibyll 2.1.
In this contribution, we have presented results from recent analyses demonstrating the unique opportunities that the IceCube Neutrino Observatory provides for measurements of the muon content of extensive air showers. An analysis of the density of mainly low-energy (O(GeV)) muons at lateral distances of 600 m and 800 m using IceTop only was presented, as well as an analysis of
Muons in air showers with IceCube 1). Top row: z-values for the measured muon densities as presented in Fig. 4. Bottom row: z-values for the measured high-energy muon numbers as presented in Fig. 7. Shown for comparison are the expectations from the cosmic-ray flux models H3a [15], GST [16], and GSF [17].
the number of high-energy (O(TeV)) muons in coincident events observed with both the IceTop and IceCube in-ice detectors. The measurement of the high-energy muons agrees within uncertainties with predictions for all tested hadronic interaction considering various cosmic-ray flux models. On the other hand, the predictions of the low-energy muon content in air showers differ more between models, leading to differences in the interpretation of the measurements, with the results based on QGSJet-II.04 and EPOS-LHC yielding a lighter-than-expected mass composition. This also results in a tension between the low-and high-energy muon results for these hadronic interaction models, indicating that they do not give a consistent description of the air-shower development.
Future improvements in the analysis are expected to provide stronger tests of the description of muon production in air showers. The high-energy muon analysis will benefit from smaller systematic uncertainties related to the modeling of the Antarctic ice. The low-energy muon analysis will benefit from the development of an event-by-event estimator of the low-energy muon content [18]. Both analyses could also be extended towards higher energies and more inclined showers. We also plan to develop analyses of the fluctuations in the muon numbers, as well as a combined analysis studying the correlation between the low-and high-energy muon content in air showers.
The existing and future muon analyses are also expected to benefit from the ongoing and future upgrades of the observatory. The installation of additional surface instrumentation, namely scintillation detectors and radio antennas, will bring improved separation between the electromagnetic and muonic shower components, as well as an independent measure of the shower energy and the depth of shower maximum [19]. The IceCube Upgrade [20] may benefit muon bundle studies with the deployment of a denser core of strings and will improve our understanding of the ice with the help of various new calibration devices. Furthermore, the plans for IceCube-Gen2 [21], with an eightfold increase in volume of the in-ice detector and a corresponding surface array [22], will bring increased statistics at high energies and a larger opening angle for coincident events.