---
title: Muon Puzzle in Cosmic Ray Showers
url: https://www.emergentmind.com/topics/muon-puzzle
type: topic
---

# Muon Puzzle in Cosmic Ray Showers

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 [2509.24021][2606.22855]. 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”** [2112.08312]. 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 [2304.00294][2105.06148].

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\) [2112.08312]. 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 [1204.5180].

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 [2606.22855]. 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 \(N_\mu^{\rm model}\approx1.455\times10^7\) [2509.24021].

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 \(\theta\ge55^\circ\) plus **6.3 khr** for \(40^\circ\le\theta<55^\circ\), it recorded about **99.6 k** events with muon multiplicity \(m\ge5\) and \(\theta\ge55^\circ\), together with **30.4 k** events in the \(40^\circ\)–\(55^\circ\) band [2208.05926]. Their local-muon-density analysis found that at moderate zenith angles, corresponding to \(E_0\sim10^{16}\,\mathrm{eV}\), the spectra lie close to proton simulations, whereas at \(\theta\gtrsim77^\circ\), probing \(E_0\sim10^{18}\,\mathrm{eV}\), the data rise above even the iron-initiated predictions. In WHISP \(z\)-parameter language, NEVOD–DECOR reports \(z\) increasing from approximately **0** at \(10^{16}\,\mathrm{eV}\) to nearly **1** by \(10^{18}\,\mathrm{eV}\), implying that standard models would require an **“extremely heavy”** composition to reproduce the muon-bundle intensity [2208.05926].

NEVOD–DECOR also measured the average muon energy in inclined bundles. For primary energies from **10 PeV to 1000 PeV**, the measured \(\langle E_\mu\rangle\) 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 **\(3.8\sigma\)** for proton and **\(2.6\sigma\)** for iron in **QGSJET-II-04** [2208.05926]. 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 \(E_0\) brings both surface and muon densities into agreement with **QGSJetII-04** proton predictions, while a **\(\sim25\%\)** upward shift of Auger’s \(E_0\) would substantially reduce the Auger excess [2304.13095]. 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 \(n_k(E,X)\) of particles of type \(k\) evolves with atmospheric depth \(X\) through coupled interaction and decay terms,
\[
\frac{\partial n_k(E,X)}{\partial X} = -\Bigl[\lambda_{\rm int,k}^{-1}(E)+\lambda_{\rm dec,k}^{-1}(E,X)\Bigr]\,n_k -\frac{\partial}{\partial E}\bigl[\mu_k(E)\,n_k\bigr] +\cdots,
\]
with inclusive production kernels \(c_{\ell\to k}\) and decay kernels \(d_{\ell\to k}\) controlling population transfer between species [2105.06148].

The Heitler–Matthews approximation isolates the dominant control parameters. After
\[
k=\frac{\ln(E_0/\xi_h)}{\ln m}
\]
generations, one obtains
\[
N_\mu(E_0,A)\approx A^{\,1-\beta}\Bigl(\frac{E_0}{\xi_h}\Bigr)^\beta,\qquad
\beta=\frac{\ln(\alpha m)}{\ln m}\lesssim1,
\]
with \(\xi_h\sim10\,\mathrm{GeV}\), \(m\sim50\), and \(\alpha\approx2/3\) as a typical hadronic-energy-retention fraction [2105.06148]. 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
\[
f_{\rm had}\equiv \frac{E_{\rm had}}{E_{\rm had}+E_{\rm EM}},
\qquad
R(\eta)\equiv \frac{\langle dE_{\rm EM}/d\eta\rangle}{\langle dE_{\rm had}/d\eta\rangle},
\qquad
f_{\rm had}=\frac{1}{1+R},
\]
because reducing early \(\pi^0\to\gamma\gamma\) feed increases the energy retained in the muon-producing hadronic sector [2304.00294].

For bundle experiments such as NEVOD–DECOR, the key observable is the **local muon density**
\[
D=\frac{m}{S_{\rm det}(\theta)},
\]
with \(m\) the counted multiplicity and \(S_{\rm det}(\theta)\) the projected detector area normal to the bundle direction. The corresponding differential local density spectrum,
\[
\frac{dF(D,\theta)}{dD}\sim\sum_{\rm events}\delta\!\left(D-\frac{m}{S_{\rm det}}\right),
\]
maps the measured muon bundles onto the primary-energy scale [2208.05926].

A widely used composition-discriminating observable is the \(z\)-parameter,
\[
z(E_0)=\frac{\ln N_\mu-\ln N_\mu^{(p)}}{\ln N_\mu^{(\mathrm{Fe})}-\ln N_\mu^{(p)}},
\]
where \(N_\mu^{(p)}\) and \(N_\mu^{(\mathrm{Fe})}\) are the pure-proton and pure-iron expectations [2105.06148]. Values \(z>z_{\rm mass}\), where \(z_{\rm mass}\) is inferred from \(X_{\max}\), 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
\[
\sqrt{s_{\rm NN}}=\sqrt{2m_N E_0}
\]
implies that \(E_0\simeq3\times10^7\,\mathrm{GeV}\) corresponds to \(\sqrt{s}\simeq8\,\mathrm{TeV}\), which is the scale identified in LHC-oriented reviews as the onset region where the muon deficit becomes apparent [2105.06148]. 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\(\leftrightarrow\)kaon swapping with \(f_{\rm thres}=1\) | Resolves \(\langle R_\mu\rangle\) without shifting \(\langle X_{\max}\rangle\); implies \(\Delta r\approx5\%-9\%\) at collider energies [2208.04266] |
| Thermalized core / EPOS core–corona | Statistical hadronization and strangeness enhancement in a core with \(f_{\rm core}\) increasing with \(\sqrt{s}\) | Reduces muon deficit from \(\sim30\%\) to \(\sim10\%-15\%\) at \(10^{19}\,\mathrm{eV}\) [2304.00294] |
| \(\pi\leftrightarrow K\) swap \(F_s\) model | Forward pseudorapidity- and energy-dependent swapping probability | Constant \(f_s\approx0.4\)–\(0.5\) reproduces Auger-level excess [2307.08634] |
| EPOS.LHC-R | Perfect isospin in fragmentation, extra neutral resonances, retuned \(\sigma_{\rm inel}\), increased elasticity and forward multiplicity | \(X_{\max}\) shifts by \(\Delta X_{\max}\approx+25\,\mathrm{g/cm^2}\); \(N_\mu\) rises by \(+7\%-10\%\); residual muon deficit \(\lesssim10\%\) at \(10^{19}\,\mathrm{eV}\) [2508.07105] |

In the **strangeball** formulation, the original fireball model was restricted to pure strangeness enhancement by setting \(f_{\rm thres}=1\), thereby eliminating the large inelasticity shift that would otherwise spoil \(\langle X_{\max}\rangle\). The effective hadronic energy fraction becomes
\[
r_{\rm eff}(E)=\bigl[1-p(E)\bigr]\,r_{\rm SM}+p(E)\,r_{\rm sb},
\]
and explicit parameter sets were found that reconcile \(\langle R_\mu\rangle\) with \(\langle X_{\max}\rangle\). 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 [2208.04266].

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 \(100\,\mathrm{TeV}\), while statistical hadronization of the core enhances strange and heavier hadron production [2304.00294]. In that framework, the hadronic energy fraction rises from about **0.45** at \(\sqrt{s}=7\,\mathrm{TeV}\) to about **0.55** at \(\sqrt{s}=50\,\mathrm{TeV}\), extrapolating to about **0.60–0.65** above \(100\,\mathrm{TeV}\), and the muon deficit at \(10^{19}\,\mathrm{eV}\) is reduced from about **30%** to about **10–15%** [2304.00294]. The mechanism is explicitly tied to strangeness enhancement, with \(K/\pi\) rising from about **0.08** at low multiplicity to about **0.18** at high multiplicity in the model.

The **\(F_s\)** phenomenology packages the same intuition in a more agnostic way. It post-processes hadronic events by converting forward pions into kaons with a probability \(F_s(E^{\rm (proj)},\eta)\). A constant benchmark with \(f_s\approx0.4\)–**0.5** is sufficient to reproduce the Auger muon excess, whereas adding elaborate energy or pseudorapidity ramps reduces efficiency and tends to force \(f_s^{(\max)}\to1\) [2307.08634].

**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 **\(\eta'\)** and **\(f_0\)** resonances that preferentially feed \(\rho^0\), lowers the \(p\)-air inelastic cross section by about **5%**, raises the elasticity by about **20%**, and increases forward charged-hadron multiplicity by about **30%** at \(\sqrt{s}=13\,\mathrm{TeV}\) [2508.07105]. The net effect is a reduction of the electromagnetic fraction,
\[
R_{\rm new}\simeq R_{\rm old}-0.05 \quad\Rightarrow\quad \Delta N_\mu\simeq+7\%,
\]
together with a deeper shower maximum and a residual disagreement with Auger of only \(\lesssim10\%\) in muon number at \(10^{19}\,\mathrm{eV}\) [2508.07105]. 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,
\[
\zeta(E_{\rm proj},x_{\rm lab})
=\Theta(E_{\rm proj}-E_{\rm start})\,
\Theta(x_{\rm lab}-x_{\rm thr})\,
\min\!\left[1,\frac{\kappa}{100}\log_{10}\!\left(\frac{E_{\rm proj}}{E_{\rm start}}\right)\right],
\]
was confronted with Auger data using **MCEQ** [2509.24021]. Matching the Auger central value requires \(\zeta(10^{10}\,\mathrm{GeV})\gtrsim0.5\), while \(\zeta\gtrsim0.3\) is sufficient at the \(1\sigma\) level. The same study concludes that **10.8%** precision on \(K/\pi\) 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 \(x_{\rm thr}>0.2\) [2509.24021].

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 \(\pi\leftrightarrow K\) swap with \(f_s=0.1\) would be detectable or excludable with high significance in a FLArE-like setup [2307.08634]. 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 \(\pi^0\to2\gamma\) [2606.22855]. In that picture the average electron count at ground is reduced,
\[
\langle \widetilde N_e\rangle = F(E_\gamma)\langle N_e\rangle,\qquad F(E_\gamma)<1,
\]
so the reconstructed primary energy is biased low while the muon content remains essentially unaffected. Requiring a **30%** apparent muon excess fixes \(F(E_\gamma\approx10^{17}\,\mathrm{eV})\approx0.8\), corresponding to about **20%** suppression of \(\sigma_{\gamma N}\) at that energy and an effective scale \(M_s/\xi\sim5\times10^{17}\,\mathrm{eV}\) in the linearized model [2606.22855]. 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 \(g-2\) puzzle”**. In that context, the Standard Model prediction
\[
a_\mu^{\rm SM}=116\,591\,810(43)\times10^{-11}
\]
combined with the low-energy \(e^+e^-\to\mathrm{hadrons}\) determination
\[
(a_\mu^{\rm HVP})_{e^+e^-}^{\rm TI}=6\,931(40)\times10^{-11}
\]
yields a **\(4.2\sigma\)** discrepancy with the experimental world average
\[
a_\mu^{\rm exp}=116\,592\,061(41)\times10^{-11},
\]
whereas the BMW lattice value
\[
(a_\mu^{\rm HVP})_{\rm BMW}=7\,075(55)\times10^{-11}
\]
reduces the discrepancy to **\(1.6\sigma\)** [2112.08312]. The specific new-physics hypothesis examined there posits a sub-GeV vector boson \(Z'\) that contaminates the measured \(e^+e^-\to\mathrm{hadrons}\) cross section. The analysis concludes that this route is excluded by **LEP II**, **BaBar**, electron \(g-2\), 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** [2112.08312].

Muon spectroscopy and low-\(Q^2\) scattering are closely connected to this broader muon-sector anomaly program. The **Mu-MASS** experiment aims at a **1000-fold improvement** in the \(1S\)–\(2S\) transition frequency of muonium. Its projected precision corresponds to \(\delta m_\mu/m_\mu\simeq1\times10^{-9}\) and \(\delta R_\infty/R_\infty\simeq1\times10^{-12}\), with a Phase-2 total systematic of \(\lesssim3\,\mathrm{kHz}\) [1811.00310]. A plausible implication is that improved muonium spectroscopy will remove parametric ambiguity in the interpretation of \(a_\mu\) and in the extraction of the proton charge radius from muonic systems.

The **MUSE** experiment addresses the proton-radius problem through simultaneous \(\mu^\pm p\) and \(e^\pm p\) scattering over
\[
Q^2\approx0.002\text{–}0.07\ \mathrm{GeV}^2,
\]
with projected relative cross-section uncertainties at the **0.1–0.3%** level and a radius precision of about **0.01 fm** [1303.2160]. Because the charge asymmetry between \(\ell^+p\) and \(\ell^-p\) 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 [2508.07105]. In precision phenomenology, the phrase can denote a different HVP inconsistency in the \(g-2\) program, where contamination of \(e^+e^-\to\mathrm{hadrons}\) by new physics appears excluded [2112.08312]. In both domains, the field has moved from anomaly cataloguing to targeted cross-checks: forward hadron measurements for air showers, and spectroscopy plus low-\(Q^2\) scattering for the muon-sector precision anomalies.

Source: https://www.emergentmind.com/topics/muon-puzzle