Measurement of the mean number of muons with energies above 500 GeV in air showers detected with the IceCube Neutrino Observatory

TL;DR

The paper addresses the Muon Puzzle by measuring the mean multiplicity of TeV muons ($E_\mu>500$ GeV) in near-vertical air showers from primaries in $E\in[2.5,100]$ PeV using IceTop–IceCube coincidences. A neural-network framework jointly reconstructs primary energy and muon content from surface and in-ice signals, followed by Monte Carlo corrections to account for composition and model dependencies across Sibyll 2.1, QGSJet-II.04, and EPOS-LHC. The results show general agreement with Sibyll 2.1 predictions, with about a 5% model-dependent variation, while EPOS-LHC exhibits tensions when compared to low-energy muon data, underscoring the complexities of high-energy hadronic interactions in air showers. These measurements provide a critical testbed for hadronic-interaction models and help calibrate mass-composition inferences, with significant implications for future large-scale detectors and accelerator-physics extrapolations. The work demonstrates the utility of combining surface and deep detectors, machine learning-based reconstructions, and MC-driven corrections to probe the high-energy muon content of air showers and advance our understanding of cosmic-ray physics.

Abstract

We present a measurement of the mean number of muons with energies larger than 500 GeV in near-vertical extensive air showers initiated by cosmic rays with primary energies between 2.5 PeV and 100 PeV. The measurement is based on events detected in coincidence between the surface and in-ice detectors of the IceCube Neutrino Observatory. Air showers are recorded on the surface by IceTop, while a bundle of high-energy muons ("TeV muons") from the shower can subsequently produce a track-like event in the IceCube in-ice array. Results are obtained assuming the hadronic interaction models Sibyll 2.1, QGSJet-II.04, and EPOS-LHC. The measured number of TeV muons is found to be in agreement with predictions from air-shower simulations. The results have also been compared to a measurement of low-energy muons by IceTop, indicating an inconsistency between the predictions for low- and high-energy muons in simulations based on the EPOS-LHC model.

Figures (18)

• Figure 1: Schematic drawing of a cosmic-ray air shower observed in coincidence between IceTop, the surface component of the IceCube Neutrino Observatory, and the IceCube in-ice array.
• Figure 2: Correlation plots showing the relation between the true and neural-network reconstructed values for primary cosmic-ray energy (left) and high-energy muon multiplicity (right). The muon multiplicity has been divided by $E^\beta$ to reduce the effect of its underlying energy dependence. The histograms include Sibyll 2.1 simulations of all four mass groups (p, He, O, and Fe).
• Figure 3: Bias and resolution of the neural network reconstructions, defined as the mean and standard deviation of the difference between the logarithms of the reconstructed and true values, shown separately for different primary cosmic-ray masses. Primary energy reconstruction is shown on the left, muon multiplicity ($E_\mu > 500G\eV$) reconstruction on the right. The analysis includes events between 6.4 and 8.0 in $\log_{10} (E/\mathrm{GeV})$, corresponding to $0.7 \lesssim \log_{10} (N_\mu) \lesssim 2.8$ for proton showers and $1.5 \lesssim \log_{10} (N_\mu) \lesssim 3.1$ for iron showers.
• Figure 4: Top: comparison between the average reconstructed TeV muon number in bins of reconstructed primary energy and the true muon number in bins of true energy in Sibyll 2.1 simulation for four different primaries. Bottom: ratio of the reconstructed and true values from the top plot, fitted with quadratic functions defining correction factors based on Sibyll 2.1.
• Figure 5: Comparison between the initial reconstruction of the average muon number (in bins of reconstructed energy), the corrected results, and the true muon number (in bins of true energy) obtained from air showers simulated with Sibyll 2.1. The values have been scaled according to Eq. (\ref{['eq:z']}), so that the true value for proton is at zero and the true value for iron is at 1. The brackets show the uncertainty assigned to the final result to account for the typical offsets remaining after the correction. This comparison is shown for four pure composition scenarios.