High-energy cosmic rays produce Extensive Air Showers (EAS) when interacting with the Earth’s atmosphere, generating, among other particles, muons that can reach underground detectors such as ALICE. Muon multiplicity provides information on the energy and composition of primary cosmic rays. However, discrepancies between experimental measurements and the predictions of hadronic interaction models, known as the Muon Puzzle, remain a challenge for the description of these air showers. During LHC Run 1, ALICE found that the rate of high-muon-multiplicity events was compatible with QGSJET-II-04 predictions under the assumption of a primary composition dominated by heavy elements. This result motivated the present work, whose objective was to analyze a larger data sample collected during Run 2 (2015–2018) and to evaluate the ability of different hadronic interaction models to describe the observed muon multiplicity.

In this work, we study the muon multiplicity distribution measured by ALICE during LHC Run 2 to evaluate whether the QGSJET-II-04, EPOS-LHC, and SIBYLL 2.3d models can describe the observations for different primary cosmic-ray compositions. The data are compared with Monte Carlo simulations performed using CORSIKA, considering proton and iron primaries as representatives of light and heavy compositions, respectively. The results show that QGSJET-II-04 reasonably reproduces the muon multiplicity distribution for $4 < N_{\mu} < 50$ when a heavy primary composition is considered, whereas SIBYLL 2.3d and EPOS-LHC exhibit muon deficits of approximately 20% and 30%, respectively, over different multiplicity intervals. For events with $N_{\mu} > 100$, the rates predicted by QGSJET-II-04 and SIBYLL 2.3d assuming an iron composition are compatible with the data, while EPOS-LHC underestimates the observed rate. Overall, the results are consistent with a predominantly heavy composition for high-muon-multiplicity events, associated with primary energies $E_{\mathrm{prim}}$ close to $10^{17}\,\mathrm{eV}$, and provide constraints on current hadronic interaction models and on possible more exotic mechanisms of muon production.