Muon Puzzle in Cosmic Ray Showers

• Muon Puzzle is a cosmic-ray phenomenon where observed muon counts exceed simulation predictions by 20–60%, exposing flaws in current hadronic models.
• Experimental programs like Auger, TA, and NEVOD–DECOR reveal the discrepancy through measures such as local muon density and the z-parameter, indicating energy-dependent composition effects.
• Proposed resolutions, including strangeness enhancement and thermalized core formation, adjust forward particle production models and are being tested via collider experiments.

The term muon puzzle most commonly denotes the persistent excess of muons measured in extensive air showers relative to the predictions of modern hadronic-interaction models. In this usage, air-shower experiments report ground-level muon counts up to 30–60% above simulations for primary energies $E_{\rm prim}\gtrsim10^{8}\,\mathrm{GeV}$ or around $10^{19}\,\mathrm{eV}$, while composition indicators such as $X_{\max}$ often favor a light or mixed primary composition rather than the extremely heavy composition that would otherwise be required to match the muon data (Ohashi et al., 28 Sep 2025, Li, 22 Jun 2026). The same phrase has also been used in other muon-related contexts, notably for the BMW-versus-dispersive hadronic-vacuum-polarization tension in the muon anomalous magnetic moment, where it appears as the “new muon $g-2$ puzzle” (Luzio et al., 2021). The dominant contemporary meaning, however, is the cosmic-ray extensive-air-shower anomaly.

1. Terminology and scope

In cosmic-ray physics, the muon puzzle is the statement that simulations based on state-of-the-art generators such as QGSJET, EPOS, and SIBYLL underpredict the number of muons produced in hadronic cascades in the atmosphere. The discrepancy is not a small correction: multiple summaries characterize it as a smooth, logarithmically rising deficit that becomes prominent above the cosmic-ray knee and reaches the 20–60% level depending on energy, observable, and model choice (Scaria et al., 2023, Albrecht et al., 2021).

The phrase has also acquired narrower specialized meanings. In precision muon phenomenology, the “new muon $g-2$ puzzle” denotes the mismatch between the BMW lattice-QCD determination of the leading hadronic vacuum polarization and the low-energy $e^+e^-\to\mathrm{hadrons}$ data customarily used in dispersive evaluations of $a_\mu$ (Luzio et al., 2021). In underground cosmic-muon studies, the title “Gran Sasso muon puzzle” was used for an unexpected 10–11 year modulation in TeV-scale muon flux, anticorrelated with the solar cycle, in addition to the standard annual modulation (Fernandez-Martinez et al., 2012).

A recurring source of confusion is the assumption that all of these usages refer to the same anomaly. They do not. The air-shower puzzle concerns hadronic cascade development and forward particle production; the $g-2$ puzzle concerns hadronic vacuum polarization in precision electroweak observables; the Gran Sasso usage concerns time-series structure in underground muon flux. The common element is the muon, not the underlying dynamics.

2. Experimental status in extensive air showers

The air-shower muon puzzle is supported by several independent observational programs. A recent working-group synthesis cited by later studies combines results from Auger, TA, KASCADE-Grande, IceCube, Yakutsk and other experiments and reconfirms a significant muon deficit in standard simulations (Li, 22 Jun 2026). In the Auger-centered formulation used in phenomenological studies, the observed muon count at ground exceeds the model prediction by up to 30–60%, with a representative Auger result quoted as

$$\langle \ln R_\mu\rangle(10^{19}\,\mathrm{eV})=0.601\pm0.016\ \mathrm{(stat)}\ ^{+0.167}_{-0.201}\ \mathrm{(sys)},$$

where $R_\mu\equiv N_\mu/N_\mu^{\rm model}$ and $10^{19}\,\mathrm{eV}$0 (Ohashi et al., 28 Sep 2025).

The NEVOD–DECOR program provides an especially detailed bundle-based characterization. Over May 2012–March 2021, with a live time of 58.3 khr for $10^{19}\,\mathrm{eV}$1 plus 6.3 khr for $10^{19}\,\mathrm{eV}$2, it recorded about 99.6 k events with muon multiplicity $10^{19}\,\mathrm{eV}$3 and $10^{19}\,\mathrm{eV}$4, together with 30.4 k events in the $10^{19}\,\mathrm{eV}$5–$10^{19}\,\mathrm{eV}$6 band (Bogdanov et al., 2022). Their local-muon-density analysis found that at moderate zenith angles, corresponding to $10^{19}\,\mathrm{eV}$7, the spectra lie close to proton simulations, whereas at $10^{19}\,\mathrm{eV}$8, probing $10^{19}\,\mathrm{eV}$9, the data rise above even the iron-initiated predictions. In WHISP $X_{\max}$0-parameter language, NEVOD–DECOR reports $X_{\max}$1 increasing from approximately 0 at $X_{\max}$2 to nearly 1 by $X_{\max}$3, implying that standard models would require an “extremely heavy” composition to reproduce the muon-bundle intensity (Bogdanov et al., 2022).

NEVOD–DECOR also measured the average muon energy in inclined bundles. For primary energies from 10 PeV to 1000 PeV, the measured $X_{\max}$4 rises from roughly 80–100 GeV to 140–180 GeV, while the model bands show a mild decrease or leveling off. In the highest-density bins the excess reaches $X_{\max}$5 for proton and $X_{\max}$6 for iron in QGSJET-II-04 (Bogdanov et al., 2022). This strengthens the anomaly from a pure multiplicity problem into a joint multiplicity-plus-spectrum problem.

Not all analyses interpret the discrepancy identically. A Yakutsk-based comparison argued that the Yakutsk muon densities are consistent with fluorescence-based composition inferences and that part of the apparent conflict may arise from energy-scale intercalibration. In that analysis, a 10% downward rescaling of Yakutsk $X_{\max}$7 brings both surface and muon densities into agreement with QGSJetII-04 proton predictions, while a $X_{\max}$8 upward shift of Auger’s $X_{\max}$9 would substantially reduce the Auger excess (Glushkov et al., 2023). This does not eliminate the broader puzzle, but it shows that the experimental status includes an unresolved calibration controversy rather than a universally identical anomaly across all observatories.

3. Observables and theoretical framework

The air-shower problem is usually formulated in terms of hadronic cascade transport or its Heitler–Matthews reduction. In a full cascade description, the number density $g-2$0 of particles of type $g-2$1 evolves with atmospheric depth $g-2$2 through coupled interaction and decay terms,

$g-2$3

with inclusive production kernels $g-2$4 and decay kernels $g-2$5 controlling population transfer between species (Albrecht et al., 2021).

The Heitler–Matthews approximation isolates the dominant control parameters. After

$g-2$6

generations, one obtains

$g-2$7

with $g-2$8, $g-2$9, and $g-2$0 as a typical hadronic-energy-retention fraction (Albrecht et al., 2021). This parameterization makes explicit why the puzzle is so sensitive to secondary composition: small modifications of the fraction of energy flowing into the hadronic branch are multiplied over several cascade generations.

Detailed simulation studies summarized in the air-shower literature identify the decisive control variable as the hadronic-to-electromagnetic energy partition, rather than the inelastic cross-section alone. A convenient measure is

$g-2$1

because reducing early $g-2$2 feed increases the energy retained in the muon-producing hadronic sector (Scaria et al., 2023).

For bundle experiments such as NEVOD–DECOR, the key observable is the local muon density

$g-2$3

with $g-2$4 the counted multiplicity and $g-2$5 the projected detector area normal to the bundle direction. The corresponding differential local density spectrum,

$g-2$6

maps the measured muon bundles onto the primary-energy scale (Bogdanov et al., 2022).

A widely used composition-discriminating observable is the $g-2$7-parameter,

$g-2$8

where $g-2$9 and $e^+e^-\to\mathrm{hadrons}$0 are the pure-proton and pure-iron expectations (Albrecht et al., 2021). Values $e^+e^-\to\mathrm{hadrons}$1, where $e^+e^-\to\mathrm{hadrons}$2 is inferred from $e^+e^-\to\mathrm{hadrons}$3, expose the muon excess independently of absolute normalization conventions.

A central empirical point is the energy scale at which the discrepancy becomes visible. The relation

$e^+e^-\to\mathrm{hadrons}$4

implies that $e^+e^-\to\mathrm{hadrons}$5 corresponds to $e^+e^-\to\mathrm{hadrons}$6, which is the scale identified in LHC-oriented reviews as the onset region where the muon deficit becomes apparent (Albrecht et al., 2021). This is why forward LHC measurements are treated as directly relevant rather than merely suggestive.

4. Hadronic-interaction resolutions

A large fraction of the current literature attempts to resolve the muon puzzle without invoking fundamentally new particles, by modifying hadronization, forward particle production, or the balance between electromagnetic and hadronic energy flow. Three broad mechanisms recur: strangeness enhancement, thermalized core formation, and global retuning of hadronic event generators.

Scenario	Main modification	Stated impact
Strangeball model	Pion$e^+e^-\to\mathrm{hadrons}$7kaon swapping with $e^+e^-\to\mathrm{hadrons}$8	Resolves $e^+e^-\to\mathrm{hadrons}$9 without shifting $a_\mu$0; implies $a_\mu$1 at collider energies (Manshanden et al., 2022)
Thermalized core / EPOS core–corona	Statistical hadronization and strangeness enhancement in a core with $a_\mu$2 increasing with $a_\mu$3	Reduces muon deficit from $a_\mu$4 to $a_\mu$5 at $a_\mu$6 (Scaria et al., 2023)
$a_\mu$7 swap $a_\mu$8 model	Forward pseudorapidity- and energy-dependent swapping probability	Constant $a_\mu$9–$g-2$0 reproduces Auger-level excess (Sciutto et al., 2023)
EPOS.LHC-R	Perfect isospin in fragmentation, extra neutral resonances, retuned $g-2$1, increased elasticity and forward multiplicity	$g-2$2 shifts by $g-2$3; $g-2$4 rises by $g-2$5; residual muon deficit $g-2$6 at $g-2$7 (Pierog et al., 9 Aug 2025)

In the strangeball formulation, the original fireball model was restricted to pure strangeness enhancement by setting $g-2$8, thereby eliminating the large inelasticity shift that would otherwise spoil $g-2$9. The effective hadronic energy fraction becomes

$$\langle \ln R_\mu\rangle(10^{19}\,\mathrm{eV})=0.601\pm0.016\ \mathrm{(stat)}\ ^{+0.167}_{-0.201}\ \mathrm{(sys)},$$0

and explicit parameter sets were found that reconcile $$\langle \ln R_\mu\rangle(10^{19}\,\mathrm{eV})=0.601\pm0.016\ \mathrm{(stat)}\ ^{+0.167}_{-0.201}\ \mathrm{(sys)},$$1 with $$\langle \ln R_\mu\rangle(10^{19}\,\mathrm{eV})=0.601\pm0.016\ \mathrm{(stat)}\ ^{+0.167}_{-0.201}\ \mathrm{(sys)},$$2. The resulting fits imply strangeball probabilities of about 0.35–0.40 at Tevatron energies and 0.40–0.45 at LHC energies, corresponding to a 5–9% increase in the energy fraction retained in the hadronic cascade (Manshanden et al., 2022).

In the thermalized-core approach, the collision is split into a dense core and a dilute corona. In EPOS LHC, the estimated core fraction grows from approximately 0.3–0.5 at few-TeV energies to 0.6–0.8 at effective UHECR energies around $$\langle \ln R_\mu\rangle(10^{19}\,\mathrm{eV})=0.601\pm0.016\ \mathrm{(stat)}\ ^{+0.167}_{-0.201}\ \mathrm{(sys)},$$3, while statistical hadronization of the core enhances strange and heavier hadron production (Scaria et al., 2023). In that framework, the hadronic energy fraction rises from about 0.45 at $$\langle \ln R_\mu\rangle(10^{19}\,\mathrm{eV})=0.601\pm0.016\ \mathrm{(stat)}\ ^{+0.167}_{-0.201}\ \mathrm{(sys)},$$4 to about 0.55 at $$\langle \ln R_\mu\rangle(10^{19}\,\mathrm{eV})=0.601\pm0.016\ \mathrm{(stat)}\ ^{+0.167}_{-0.201}\ \mathrm{(sys)},$$5, extrapolating to about 0.60–0.65 above $$\langle \ln R_\mu\rangle(10^{19}\,\mathrm{eV})=0.601\pm0.016\ \mathrm{(stat)}\ ^{+0.167}_{-0.201}\ \mathrm{(sys)},$$6, and the muon deficit at $$\langle \ln R_\mu\rangle(10^{19}\,\mathrm{eV})=0.601\pm0.016\ \mathrm{(stat)}\ ^{+0.167}_{-0.201}\ \mathrm{(sys)},$$7 is reduced from about 30% to about 10–15% (Scaria et al., 2023). The mechanism is explicitly tied to strangeness enhancement, with $$\langle \ln R_\mu\rangle(10^{19}\,\mathrm{eV})=0.601\pm0.016\ \mathrm{(stat)}\ ^{+0.167}_{-0.201}\ \mathrm{(sys)},$$8 rising from about 0.08 at low multiplicity to about 0.18 at high multiplicity in the model.

The $$\langle \ln R_\mu\rangle(10^{19}\,\mathrm{eV})=0.601\pm0.016\ \mathrm{(stat)}\ ^{+0.167}_{-0.201}\ \mathrm{(sys)},$$9 phenomenology packages the same intuition in a more agnostic way. It post-processes hadronic events by converting forward pions into kaons with a probability $R_\mu\equiv N_\mu/N_\mu^{\rm model}$0. A constant benchmark with $R_\mu\equiv N_\mu/N_\mu^{\rm model}$1–0.5 is sufficient to reproduce the Auger muon excess, whereas adding elaborate energy or pseudorapidity ramps reduces efficiency and tends to force $R_\mu\equiv N_\mu/N_\mu^{\rm model}$2 (Sciutto et al., 2023).

EPOS.LHC-R takes a different route. It retains a global QCD framework but changes the correlation between measured mid-rapidity data and forward particle production. The model restores exact isospin in the corona fragmentation, adds $R_\mu\equiv N_\mu/N_\mu^{\rm model}$3 and $R_\mu\equiv N_\mu/N_\mu^{\rm model}$4 resonances that preferentially feed $R_\mu\equiv N_\mu/N_\mu^{\rm model}$5, lowers the $R_\mu\equiv N_\mu/N_\mu^{\rm model}$6-air inelastic cross section by about 5%, raises the elasticity by about 20%, and increases forward charged-hadron multiplicity by about 30% at $R_\mu\equiv N_\mu/N_\mu^{\rm model}$7 (Pierog et al., 9 Aug 2025). The net effect is a reduction of the electromagnetic fraction,

$R_\mu\equiv N_\mu/N_\mu^{\rm model}$8

together with a deeper shower maximum and a residual disagreement with Auger of only $R_\mu\equiv N_\mu/N_\mu^{\rm model}$9 in muon number at $10^{19}\,\mathrm{eV}$00 (Pierog et al., 9 Aug 2025). This suggests that at least part of the anomaly may reflect incomplete forward-QCD modeling rather than a fundamentally missing degree of freedom.

5. Collider tests and nonstandard scenarios

The current phase of the subject is increasingly falsification-oriented. Rather than treating the muon puzzle as an atmospheric anomaly alone, several works ask whether proposed fixes can be tested directly at colliders.

For the strangeness-enhancement scenario, the key control variable is the forward kaon-to-pion ratio. A three-parameter model,

$10^{19}\,\mathrm{eV}$01

was confronted with Auger data using MCEQ (Ohashi et al., 28 Sep 2025). Matching the Auger central value requires $10^{19}\,\mathrm{eV}$02, while $10^{19}\,\mathrm{eV}$03 is sufficient at the $10^{19}\,\mathrm{eV}$04 level. The same study concludes that 10.8% precision on $10^{19}\,\mathrm{eV}$05 at LHCb and 8.4% at FASER would be sufficient to test the viable parameter space, and that a null result in both experiments would exclude nearly the entire Auger-compatible region except models with $10^{19}\,\mathrm{eV}$06 (Ohashi et al., 28 Sep 2025).

The Forward Physics Facility extends this logic. Because future neutrino detectors there will sample the extreme forward region, they can constrain the pion and kaon flux normalizations at the sub-percent level in idealized statistical projections, and even a benchmark $10^{19}\,\mathrm{eV}$07 swap with $10^{19}\,\mathrm{eV}$08 would be detectable or excludable with high significance in a FLArE-like setup (Sciutto et al., 2023). This establishes a direct collider-to-air-shower feedback loop: the same forward hadron chemistry that controls muon production in cascades should leave observable traces in forward neutrino yields.

More speculative work treats the muon puzzle as a possible signal of physics beyond standard hadronic modeling. One example is the D-foam scenario, in which subluminal Lorentz violation for photons suppresses Bethe–Heitler pair production in the electromagnetic subshowers initiated by $10^{19}\,\mathrm{eV}$09 (Li, 22 Jun 2026). In that picture the average electron count at ground is reduced,

$10^{19}\,\mathrm{eV}$10

so the reconstructed primary energy is biased low while the muon content remains essentially unaffected. Requiring a 30% apparent muon excess fixes $10^{19}\,\mathrm{eV}$11, corresponding to about 20% suppression of $10^{19}\,\mathrm{eV}$12 at that energy and an effective scale $10^{19}\,\mathrm{eV}$13 in the linearized model (Li, 22 Jun 2026). The same work emphasizes that a definitive assessment requires full shower simulations with modified cross sections and that the existence of the anomaly itself should continue to be tested observationally.

A common misconception is that all new-physics proposals are equally unconstrained. The literature instead shows the opposite trend: the more a scenario makes quantitative contact with shower observables, the more it becomes testable through LHCb, FASER, FPF, or dedicated shower-level Monte Carlo studies.

6. Related precision-muon anomalies and broader usage

Outside cosmic-ray physics, the best-known precision-phenomenology usage is the “new muon $10^{19}\,\mathrm{eV}$14 puzzle”. In that context, the Standard Model prediction

$10^{19}\,\mathrm{eV}$15

combined with the low-energy $10^{19}\,\mathrm{eV}$16 determination

$10^{19}\,\mathrm{eV}$17

yields a $10^{19}\,\mathrm{eV}$18 discrepancy with the experimental world average

$10^{19}\,\mathrm{eV}$19

whereas the BMW lattice value

$10^{19}\,\mathrm{eV}$20

reduces the discrepancy to $10^{19}\,\mathrm{eV}$21 (Luzio et al., 2021). The specific new-physics hypothesis examined there posits a sub-GeV vector boson $10^{19}\,\mathrm{eV}$22 that contaminates the measured $10^{19}\,\mathrm{eV}$23 cross section. The analysis concludes that this route is excluded by LEP II, BaBar, electron $10^{19}\,\mathrm{eV}$24, electroweak precision observables, and isospin-breaking constraints, so that the resolution must instead come from improved lattice results or independent HVP measurements such as MUonE (Luzio et al., 2021).

Muon spectroscopy and low-$10^{19}\,\mathrm{eV}$25 scattering are closely connected to this broader muon-sector anomaly program. The Mu-MASS experiment aims at a 1000-fold improvement in the $10^{19}\,\mathrm{eV}$26–$10^{19}\,\mathrm{eV}$27 transition frequency of muonium. Its projected precision corresponds to $10^{19}\,\mathrm{eV}$28 and $10^{19}\,\mathrm{eV}$29, with a Phase-2 total systematic of $10^{19}\,\mathrm{eV}$30 (Crivelli, 2018). A plausible implication is that improved muonium spectroscopy will remove parametric ambiguity in the interpretation of $10^{19}\,\mathrm{eV}$31 and in the extraction of the proton charge radius from muonic systems.

The MUSE experiment addresses the proton-radius problem through simultaneous $10^{19}\,\mathrm{eV}$32 and $10^{19}\,\mathrm{eV}$33 scattering over

$10^{19}\,\mathrm{eV}$34

with projected relative cross-section uncertainties at the 0.1–0.3% level and a radius precision of about 0.01 fm (Gilman et al., 2013). Because the charge asymmetry between $10^{19}\,\mathrm{eV}$35 and $10^{19}\,\mathrm{eV}$36 isolates two-photon exchange, MUSE provides a direct test of muon-specific versus purely hadronic explanations of the proton-radius discrepancy. This broader precision program is distinct from the air-shower muon puzzle, but it explains why “muon puzzle” can denote more than one unresolved problem in the literature.

Taken together, the contemporary record supports a differentiated conclusion. In cosmic-ray physics, the muon puzzle is a robust discrepancy between measured and simulated air-shower muons, sharpened by forward-QCD and composition tensions and only partially mitigated by present retunings (Pierog et al., 9 Aug 2025). In precision phenomenology, the phrase can denote a different HVP inconsistency in the $10^{19}\,\mathrm{eV}$37 program, where contamination of $10^{19}\,\mathrm{eV}$38 by new physics appears excluded (Luzio et al., 2021). In both domains, the field has moved from anomaly cataloguing to targeted cross-checks: forward hadron measurements for air showers, and spectroscopy plus low-$10^{19}\,\mathrm{eV}$39 scattering for the muon-sector precision anomalies.