A combined analysis of muon density measurements from eight leading air-shower experiments, compared on a common simulation-based reference scale (the z-scale) with cross-calibrated energy scales, established that post-LHC hadronic interaction models underpredict the measured muon content of extensive air showers, with a deficit that starts above about 10 PeV (101610^{16} eV) and grows with shower energy at 8σ8\sigma significance. The open problem is to identify which specific properties of high-energy hadronic interactions — features of forward multiparticle production such as forward multiplicities, kaon/strangeness or baryon-antibaryon yields, leading-particle spectra (elasticity), rho-meson production, or transverse-momentum distributions — are responsible for this growing deficit, and to correct the interaction models accordingly, while remaining consistent with all air-shower observables (muon content and its energy trend, muon-number fluctuations, shower-maximum depth) and with accelerator constraints, or else to demonstrate quantitatively that a non-hadronic origin (energy-scale shifts, composition reference, experimental biases) accounts for the discrepancy.

Background

Cosmic rays with energies from PeV (101510^{15} eV, detected indirectly) to tens of EeV (101910^{19} eV) are studied through the extensive air showers (EAS) they trigger in the atmosphere: hadronic cascades whose electromagnetic component (electrons, positrons, photons from neutral-pion decay) dominates the energy deposit, and whose muonic component arises from the decay of charged pions and kaons produced in hadronic collisions. Interpreting air-shower measurements requires event generators ('hadronic interaction models' such as EPOS-LHC, QGSJet-II.04 and SIBYLL) that extrapolate accelerator data on multiparticle production to energies and, crucially, to forward emission angles that no collider covers.

In the Heitler-Matthews toy model of shower development, the muon number at ground scales as a power law in primary energy EE and rises with the mass number AA of the primary as Nμ∝A1−βEβN_\mu \propto A^{1-\beta} E^{\beta} with β≈0.9\beta \approx 0.9, so the muon content is a mass-composition estimator — but only as trustworthy as the hadronic model. The mean logarithmic mass ⟨ln⁡A⟩\langle\ln A\rangle can be inferred independently from Xmax⁡X_{\max}, the atmospheric depth of the shower maximum measured optically by fluorescence and Cherenkov detectors. In 2016 the Pierre Auger Observatory reported, with a method insensitive to the absolute energy calibration, that the hadronic signal of showers at 6-16 EeV exceeds the predictions of EPOS-LHC and QGSJet-II.04 by factors 1.33+/-0.16 and 1.61+/-0.21 respectively, with a corresponding muon excess; at 10^16-10^17 eV, however, KASCADE-Grande and EAS-MSU reported no muon discrepancy while SUGAR did, leaving a confused picture.

To resolve this, the Working group on Hadronic Interactions and Shower Physics (WHISP) combined muon density measurements from eight experiments (EAS-MSU, IceCube, KASCADE-Grande, NEVOD-DECOR, the Pierre Auger Observatory and its AMIGA muon-detector infill, SUGAR, Telescope Array, Yakutsk) covering PeV to tens of EeV (Dembinski et al., 2019, arXiv:1902.08124). Comparisons use the z-scale,

z=ln⁡(Nμdet/Nμ,pdet)ln⁡(Nμ,Fedet/Nμ,pdet)z = \frac{\ln\left(N_\mu^{\rm det}/N_{\mu,p}^{\rm det}\right)}{\ln\left(N_{\mu,\rm Fe}^{\rm det}/N_{\mu,p}^{\rm det}\right)}

where NμdetN_\mu^{\rm det} is the detected muon density and Nμ,pdetN_{\mu,p}^{\rm det}, Nμ,FedetN_{\mu,\rm Fe}^{\rm det} are simulated muon densities for proton and iron primaries after full detector simulation: z is 0 for pure proton showers and 1 for pure iron if data agree with simulations, and many detector-specific biases cancel. Experiment energy scales, independently calibrated only to 10-20%, were cross-calibrated using the isotropic all-particle cosmic-ray flux; the expected value for the actual mean composition is zmass=⟨ln⁡A⟩/ln⁡56z_{\rm mass} = \langle\ln A\rangle/\ln 56, and the deviation Δz=z−zmass\Delta z = z - z_{\rm mass} measures the model deficit. After cross-calibration the picture is remarkably consistent: muon measurements agree with the post-LHC models EPOS-LHC, QGSJet-II.04 and SIBYLL-2.3 up to about 101610^{16} eV; above that, z rises above expectation with slope 0.22-0.35 per decade in energy for EPOS-LHC and QGSJet-II.04, significant at 8 sigma. The analysis identified known limitations: the energy-scale reference carries an overall uncertainty that allows all z-values to shift collectively by about +/-0.25; EAS-MSU and KASCADE-Grande could not be cross-calibrated; and zmassz_{\rm mass} inherits model dependence through the conversion of Xmax⁡X_{\max} data to ⟨ln⁡A⟩\langle\ln A\rangle.

The successor review (Albrecht et al., 2022) traced the onset of the deficit to a nucleon-nucleon center-of-mass energy of about 8 TeV — inside LHC reach — and argued that the most likely explanation is a modification of hadron composition in soft multiparticle production that reduces the fraction of hadronic energy diverted to the electromagnetic cascade through neutral-pion decay, plausibly connected to the multiplicity-dependent strangeness enhancement observed by ALICE at mid-rapidity, but requiring confirmation in the forward region, at large pseudorapidity (η>2\eta > 2, the collider angular coordinate; large η\eta means emission close to the beam direction) where air-shower-relevant hadrons are produced. Because muon production compounds over many cascade generations, even small per-interaction deviations accumulate; the measured fluctuations of the muon number (Pierre Auger Collaboration, 2021) agree with simulations, indicating many small accumulating deviations rather than one gross error in the first interaction.

Problem Statement

Identify the physical origin of the energy-dependent muon deficit in extensive-air-shower simulations: determine which specific properties of high-energy hadronic interactions — features of forward multiparticle production such as forward particle multiplicities, kaon/strangeness and baryon-antibaryon yields, ρ0\rho^0-meson production, leading-particle spectra (elasticity), or transverse-momentum distributions, and at which energies and phase-space regions they act — are responsible for the growing shortfall of simulated muons relative to the cross-calibrated measurements above about 101610^{16} eV, and quantitatively account for the observed trend of the deviation Δz=z−zmass\Delta z = z - z_{\mathrm{mass}} (slope 0.22-0.35 per decade in energy, significant at 8σ8\sigma for EPOS-LHC and QGSJet-II.04).

A resolution of the problem takes one of the following forms, subject to the stated conditions:

  1. A demonstrated mechanism: an identified modification of hadron production, implemented in a hadronic interaction model, that simultaneously reproduces (i) the measured muon densities and their energy trend on the z-scale for the cross-calibrated experiments over the PeV-to-tens-of-EeV range (the six experiments whose energy-scale-corrected z-values enter the fitted Δz\Delta z trend; EAS-MSU and KASCADE-Grande, which could not be cross-calibrated, are included only for comparison), (ii) the measured fluctuations of the muon number, (iii) the electromagnetic observables (mean and fluctuations of Xmax), and (iv) existing accelerator constraints on particle production (central and forward LHC data and fixed-target data), with the mechanism's distinctive forward-region predictions stated so they can be tested against future collider measurements; or

  2. A quantitative demonstration that a non-hadronic origin — a global energy-scale shift within the roughly 10% cross-calibration uncertainty, the composition reference used for z_mass, or correlated experimental biases — accounts for the full observed Delta z trend across all cross-calibrated experiments simultaneously, which would refute the hadronic-origin interpretation; or

  3. A forward hadron-production measurement (e.g., forward kaon-to-pion, baryon, or ρ0\rho^0 yields at LHC energies or above) that, propagated through air-shower simulations, closes the deficit in a way consistent with conditions (i)-(iv) of form 1.

The scientific objective uniting these routes is to make the muon content of air showers a trustworthy cosmic-ray mass-composition observable again by bringing the hadronic interaction models into agreement with the measured muon data without breaking any other air-shower or accelerator observable. Partial phenomenological retunes that fit the muon data but leave the responsible physics unidentified do not by themselves resolve the problem; the identification must be explicit and testable.

References

  1. H. P. Dembinski, J. C. Arteaga-Velazquez, L. Cazon, R. Conceicao, J. Gonzalez, Y. Itow, D. Ivanov, N. N. Kalmykov, I. Karpikov, S. Mueller, T. Pierog, F. Riehn, M. Roth, T. Sako, D. Soldin, R. Takeishi, G. Thompson, S. Troitsky, I. Yashin, E. Zadeba, Y. Zhezher, for the EAS-MSU, IceCube, KASCADE-Grande, NEVOD-DECOR, Pierre Auger, SUGAR, Telescope Array, and Yakutsk EAS Array collaborations, 'Report on Tests and Measurements of Hadronic Interaction Properties with Air Showers' (2019), EPJ Web Conf. 210, 020004, DOI 10.1051/epjconf/201921002004, arXiv:1902.08124, https://arxiv.org/abs/1902.08124
  2. L. Cazon, 'Working Group Report on the Combined Analysis of Muon Density Measurements from Eight Air Shower Experiments' (2020), PoS(ICRC2019)214, arXiv:2001.07508, https://arxiv.org/abs/2001.07508
  3. D. Soldin, 'Update on the Combined Analysis of Muon Measurements from Nine Air Shower Experiments' (2021), PoS(ICRC2021)349, arXiv:2108.08341, https://arxiv.org/abs/2108.08341
  4. Pierre Auger Collaboration (A. Aab et al.), 'Testing Hadronic Interactions at Ultrahigh Energies with Air Showers Measured by the Pierre Auger Observatory' (2016), Phys. Rev. Lett. 117, 192001, DOI 10.1103/PhysRevLett.117.192001, arXiv:1610.08509, https://arxiv.org/abs/1610.08509
  5. J. Albrecht, L. Cazon, H. Dembinski, A. Fedynitch, K.-H. Kampert, T. Pierog, W. Rhode, D. Soldin, B. Spaan, R. Ulrich, M. Unger, 'The Muon Puzzle in cosmic-ray induced air showers and its connection to the Large Hadron Collider' (2022), Astrophys. Space Sci. 367, 27, DOI 10.1007/s10509-022-04054-5 (preprint 2021), arXiv:2105.06148, https://arxiv.org/abs/2105.06148
  6. Pierre Auger Collaboration (A. Aab et al.), 'Measurement of the fluctuations in the number of muons in extensive air showers with the Pierre Auger Observatory' (2021), Phys. Rev. Lett. 126, 152002, DOI 10.1103/PhysRevLett.126.152002, arXiv:2102.07797, https://arxiv.org/abs/2102.07797
  7. L. A. Anchordoqui, C. Garcia Canal, F. Kling, S. J. Sciutto, J. F. Soriano, 'An explanation of the muon puzzle of ultrahigh-energy cosmic rays and the role of the Forward Physics Facility for model improvement' (2022), J. High Energy Astrophys., DOI 10.1016/j.jheap.2022.03.004, arXiv:2202.03095, https://arxiv.org/abs/2202.03095
  8. S. J. Sciutto, L. A. Anchordoqui, C. Garcia Canal, F. Kling, J. F. Soriano, 'Gauging the cosmic ray muon puzzle with the Forward Physics Facility' (2023), arXiv:2307.08634, https://arxiv.org/abs/2307.08634
  9. F. Riehn, A. Fedynitch, R. Engel, 'Sibyll★' (phenomenologically modified versions of Sibyll 2.3d for muon-puzzle studies) (2024), Astropart. Phys. 160, 102964, DOI 10.1016/j.astropartphys.2024.102964, arXiv:2404.02636, https://arxiv.org/abs/2404.02636
  10. T. Pierog, K. Werner, 'EPOS LHC-R: a global approach to solve the muon puzzle' (2025), arXiv:2508.07105, https://arxiv.org/abs/2508.07105
  11. K. Ohashi, A. Fedynitch, H. Menjo, 'Constraining the strangeness enhancement scenario of the UHECR muon puzzle with LHC experiments' (2025), arXiv:2509.24021, https://arxiv.org/abs/2509.24021
  12. Pierre Auger Collaboration (A. Abdul Halim et al.), 'Measuring the muon content of inclined air showers using AERA and the water-Cherenkov detectors of the Pierre Auger Observatory' (2025), Phys. Rev. D 112, 123042, DOI 10.1103/2q9f-pbrp, arXiv:2507.02558, https://arxiv.org/abs/2507.02558
  13. Pierre Auger Collaboration (A. Abdul Halim et al.), 'Proton-air interaction properties at √s ≈ 100 TeV from shower-depth measurements with the Pierre Auger Observatory and their connection to the Muon Puzzle' (2026), arXiv:2607.22323, https://arxiv.org/abs/2607.22323
  14. E. dos Santos, 'Probing hadronic interactions at the 100 TeV scale with the Pierre Auger Observatory' (2025), arXiv:2512.14592, https://arxiv.org/abs/2512.14592
  15. IceCube Collaboration (R. Abbasi et al.), 'Measurement of the mean number of muons with energies above 500 GeV in air showers detected with the IceCube Neutrino Observatory' (2025), Phys. Rev. D 112, 082004, DOI 10.1103/lrjy-3hht, arXiv:2506.19241, https://arxiv.org/abs/2506.19241
  16. S. Verpoost, 'Multiplicity of TeV muons in extensive air showers detected with IceTop and IceCube' (2023), PoS(ICRC2023)207, arXiv:2307.14689, https://arxiv.org/abs/2307.14689
  17. M. Reininghaus, R. Ulrich, T. Pierog, 'Air shower genealogy for muon production' (2021), PoS(ICRC2021)463, DOI 10.22323/1.395.0463, arXiv:2108.03266, https://arxiv.org/abs/2108.03266
  18. A. M. Botti, I. A. Goos, M. Perlin, T. Pierog, 'Study of the muon component in the core-corona model using CONEX 3D' (2024), arXiv:2411.06918, https://arxiv.org/abs/2411.06918
  19. J. Matthews, 'A Heitler model of extensive air showers' (2005), Astropart. Phys. 22, 387-397, DOI 10.1016/j.astropartphys.2004.09.003, https://doi.org/10.1016/j.astropartphys.2004.09.003