Muon Lateral Distribution Function (MLDF)

• Muon Lateral Distribution Function (MLDF) is the radial profile of muon density in extensive air showers, essential for energy reconstruction and composition estimation.
• Various parameterizations, such as Greisen-type forms and scaling laws, are applied to describe and fit the muon distributions with reference distances to reduce fluctuations.
• Comparative studies of MLDF reveal shape discrepancies and normalization challenges, highlighting their impact on hadronic model tests and the muon puzzle.

The muon lateral distribution function (MLDF) is the radial dependence of the muon component of an extensive air shower at ground, usually expressed as a muon density $\rho_\mu(r)$, an expected detector count $\mu(r)$, or a detector signal $S_\mu(r)$ as a function of distance from the shower axis. It is a central observable for reconstructing primary energy, inferring mass composition, and testing hadronic interaction models, because both the normalization and the shape of the muon distribution depend on shower development, detector threshold, and geometry. Across the literature, the MLDF appears in several mathematically distinct but operationally related forms: empirical mean lateral distributions, Greisen- or NKG-like parameterizations anchored at a reference distance, one-parameter scaling laws, and full likelihood models that treat electromagnetic and muonic signals jointly on an event-by-event basis (Raikin et al., 2018, Weyrauch et al., 2023, Covilakam et al., 2024).

1. Formal definitions and parameterization families

A widely used scaling representation writes the muon density as

$$\rho_\mu(r;E,X)=\frac{N_\mu(E,X)}{R_{0,\mu}^2(E,X)}\,F_\mu\!\left(\frac{r}{R_{0,\mu}(E,X)}\right),$$

where $N_\mu(E,X)$ is the total number of muons above threshold, $R_{0,\mu}(E,X)$ is a radial scale factor, and $F_\mu$ is an invariant shape function. In the same framework, a canonical baseline for the invariant shape is

$$F(x)=C\,x^{-\alpha}(1+x)^{-(\beta-\alpha)}\left(1+\left(\frac{x}{10}\right)^{\gamma}\right)^{-\delta},$$

with $C=0.28$, $\alpha=1.2$, $\mu(r)$0, $\mu(r)$1, and $\mu(r)$2; the fitted $\mu(r)$3 is reported to depend only weakly on the explicit choice of $\mu(r)$4 (Raikin et al., 2018).

Other analyses use fixed-shape reference-distance forms rather than a universal scaling function. In the SUGAR analysis, the empirical muon LDF is

$\mu(r)$5

with $\mu(r)$6, $\mu(r)$7, $\mu(r)$8, and

$\mu(r)$9

SUGAR used this form to determine $S_\mu(r)$0 and then normalized densities by $S_\mu(r)$1 to study shape independently of absolute normalization (Kalmykov et al., 2022).

Yakutsk analyses illustrate both empirical and analytic strategies. For $S_\mu(r)$2 eV, the Yakutsk array constructed empirical mean MLDFs from binned detector responses and emphasized reference-distance quantities such as $S_\mu(r)$3 and $S_\mu(r)$4 rather than introducing a closed-form global $S_\mu(r)$5 fit (Glushkov et al., 2023). For showers around $S_\mu(r)$6 eV, Yakutsk also used an explicit reference-distance parameterization,

$S_\mu(r)$7

so that the amplitude is directly the muon density at $S_\mu(r)$8 m (Glushkov et al., 2023).

Reference-distance Greisen-type forms also appear in IceTop and Auger-related reconstructions. IceTop parameterized large-distance muon densities as

$S_\mu(r)$9

with $$\rho_\mu(r;E,X)=\frac{N_\mu(E,X)}{R_{0,\mu}^2(E,X)}\,F_\mu\!\left(\frac{r}{R_{0,\mu}(E,X)}\right),$$0 m, whereas segmented-counter analyses for AMIGA-like detectors used KASCADE-Grande-like forms with fixed $$\rho_\mu(r;E,X)=\frac{N_\mu(E,X)}{R_{0,\mu}^2(E,X)}\,F_\mu\!\left(\frac{r}{R_{0,\mu}(E,X)}\right),$$1 and $$\rho_\mu(r;E,X)=\frac{N_\mu(E,X)}{R_{0,\mu}^2(E,X)}\,F_\mu\!\left(\frac{r}{R_{0,\mu}(E,X)}\right),$$2 and free normalization and slope parameters (Gonzalez, 2015, Ravignani et al., 2014). This variety suggests that “MLDF” denotes a family of closely related observables rather than a single canonical formula.

2. Measurement frameworks and reconstruction strategies

Experimentally, the MLDF is detector-defined. At Yakutsk, two particle components were analyzed simultaneously: $$\rho_\mu(r;E,X)=\frac{N_\mu(E,X)}{R_{0,\mu}^2(E,X)}\,F_\mu\!\left(\frac{r}{R_{0,\mu}(E,X)}\right),$$3 from underground scintillation muon detectors with threshold approximately $$\rho_\mu(r;E,X)=\frac{N_\mu(E,X)}{R_{0,\mu}^2(E,X)}\,F_\mu\!\left(\frac{r}{R_{0,\mu}(E,X)}\right),$$4 GeV, calibrated with the cosmic muon background, and $$\rho_\mu(r;E,X)=\frac{N_\mu(E,X)}{R_{0,\mu}^2(E,X)}\,F_\mu\!\left(\frac{r}{R_{0,\mu}(E,X)}\right),$$5 from surface scintillation detectors. Events were selected with axes within a $$\rho_\mu(r;E,X)=\frac{N_\mu(E,X)}{R_{0,\mu}^2(E,X)}\,F_\mu\!\left(\frac{r}{R_{0,\mu}(E,X)}\right),$$6 km radius around the array center, axis reconstruction accuracy better than $$\rho_\mu(r;E,X)=\frac{N_\mu(E,X)}{R_{0,\mu}^2(E,X)}\,F_\mu\!\left(\frac{r}{R_{0,\mu}(E,X)}\right),$$7 m, zenith angles $$\rho_\mu(r;E,X)=\frac{N_\mu(E,X)}{R_{0,\mu}^2(E,X)}\,F_\mu\!\left(\frac{r}{R_{0,\mu}(E,X)}\right),$$8–$$\rho_\mu(r;E,X)=\frac{N_\mu(E,X)}{R_{0,\mu}^2(E,X)}\,F_\mu\!\left(\frac{r}{R_{0,\mu}(E,X)}\right),$$9, energy bins of width $N_\mu(E,X)$0, and radial binning $N_\mu(E,X)$1; CORSIKA simulations with three high-energy hadronic models, FLUKA, EGS4, and thinning at $N_\mu(E,X)$2–$N_\mu(E,X)$3 were used to build averaged detector-response MLDFs (Glushkov et al., 2023).

In IceTop’s large-distance muon analysis, the MLDF was reconstructed from tank charge histograms at fixed $N_\mu(E,X)$4, $N_\mu(E,X)$5, and $N_\mu(E,X)$6. The per-muon charge response was obtained from Geant4, the number of muons per tank was modeled as a Poisson mixture,

$N_\mu(E,X)$7

and the mean density followed from

$N_\mu(E,X)$8

This procedure was applied in bins of $N_\mu(E,X)$9, $R_{0,\mu}(E,X)$0, and $R_{0,\mu}(E,X)$1, with emphasis on $R_{0,\mu}(E,X)$2 m where the electromagnetic component is strongly suppressed (Gonzalez, 2015).

More recent IceTop reconstructions use an explicitly two-component model,

$R_{0,\mu}(E,X)$3

with a Double Logarithmic Parabola for the electromagnetic component and a Greisen-based muon term normalized at $R_{0,\mu}(E,X)$4 m. The likelihood includes HLC and SLC tanks, silent and saturated detectors, Geant4-based muon signal PDFs, and snow attenuation applied to the electromagnetic component only (Weyrauch et al., 2023). A later reconstruction variant places the muon normalization at $R_{0,\mu}(E,X)$5 m and calibrates it to the true low-energy muon number for $R_{0,\mu}(E,X)$6 MeV and $R_{0,\mu}(E,X)$7 km (Weyrauch, 16 Sep 2025).

Underground segmented counters lead to a different statistical formulation. For AMIGA-like or Auger UMD analyses, the expected mean number of muons in a station is

$R_{0,\mu}(E,X)$8

with Poisson multiplicity

$R_{0,\mu}(E,X)$9

If a station has $F_\mu$0 segments and $F_\mu$1 fired segments in a time bin, the occupancy likelihood is

$F_\mu$2

Binary, time-resolved profile, integrator, and joint binary+ADC likelihoods have all been developed for this detector class, with the most recent combined approach using both acquisition modes simultaneously at each station (Ravignani et al., 2014, Ravignani et al., 2015, Covilakam et al., 2024, Supanitsky et al., 3 Oct 2025).

3. Reference distances and composition-sensitive observables

A recurring theme in MLDF work is the use of reference distances chosen to reduce fluctuations, stabilize fits, or suppress electromagnetic contamination. Yakutsk centers its high-energy composition analysis on $F_\mu$3 m, extracting $F_\mu$4 and $F_\mu$5, while $F_\mu$6 m is used for direct comparison with Auger’s underground muon measurement at $F_\mu$7 zenith (Glushkov et al., 2023). IceTop likewise uses a large reference radius for the muon component, with $F_\mu$8 m in the large-distance MLDF analysis and $F_\mu$9 m in the two-component reconstruction because this reduces electromagnetic background while avoiding excessive fluctuations (Gonzalez, 2015, Weyrauch et al., 2023). In the newer IceTop muon-number reconstruction, the normalization is instead defined at $$F(x)=C\,x^{-\alpha}(1+x)^{-(\beta-\alpha)}\left(1+\left(\frac{x}{10}\right)^{\gamma}\right)^{-\delta},$$0 m (Weyrauch, 16 Sep 2025). For the $$F(x)=C\,x^{-\alpha}(1+x)^{-(\beta-\alpha)}\left(1+\left(\frac{x}{10}\right)^{\gamma}\right)^{-\delta},$$1 m-spaced Auger underground array, $$F(x)=C\,x^{-\alpha}(1+x)^{-(\beta-\alpha)}\left(1+\left(\frac{x}{10}\right)^{\gamma}\right)^{-\delta},$$2 m is used because the relative fluctuations of the inferred muon density are minimized there (Covilakam et al., 2024).

Once a reference distance is fixed, MLDF studies commonly convert shape information into composition-sensitive scalar observables. Yakutsk defines the muon fraction at $$F(x)=C\,x^{-\alpha}(1+x)^{-(\beta-\alpha)}\left(1+\left(\frac{x}{10}\right)^{\gamma}\right)^{-\delta},$$3 m as

$$F(x)=C\,x^{-\alpha}(1+x)^{-(\beta-\alpha)}\left(1+\left(\frac{x}{10}\right)^{\gamma}\right)^{-\delta},$$4

and then constructs the inter-experiment scaling variable

$$F(x)=C\,x^{-\alpha}(1+x)^{-(\beta-\alpha)}\left(1+\left(\frac{x}{10}\right)^{\gamma}\right)^{-\delta},$$5

Because $$F(x)=C\,x^{-\alpha}(1+x)^{-(\beta-\alpha)}\left(1+\left(\frac{x}{10}\right)^{\gamma}\right)^{-\delta},$$6 is formed from the same shower sample, it is designed to be insensitive to the adopted absolute energy scale within a bin (Glushkov et al., 2023).

SUGAR uses an analogous normalization logic, but with eventwise densities multiplied by $$F(x)=C\,x^{-\alpha}(1+x)^{-(\beta-\alpha)}\left(1+\left(\frac{x}{10}\right)^{\gamma}\right)^{-\delta},$$7 so that the subsequent comparison isolates the radial shape rather than the total muon number (Kalmykov et al., 2022). IceTop’s large-distance analysis fits $$F(x)=C\,x^{-\alpha}(1+x)^{-(\beta-\alpha)}\left(1+\left(\frac{x}{10}\right)^{\gamma}\right)^{-\delta},$$8 and $$F(x)=C\,x^{-\alpha}(1+x)^{-(\beta-\alpha)}\left(1+\left(\frac{x}{10}\right)^{\gamma}\right)^{-\delta},$$9 in each $C=0.28$0 bin and reports that, for near-vertical showers, the reference density follows

$C=0.28$1

in agreement with the Akeno result (Gonzalez, 2015). In the event-by-event IceTop reconstruction, $C=0.28$2 or $C=0.28$3 plays the role of a low-energy muon proxy, while $C=0.28$4 remains the energy proxy (Weyrauch et al., 2023, Weyrauch, 16 Sep 2025).

These reference observables are not interchangeable across experiments. Their numerical values depend on detector threshold, shielding, zenith-angle handling, and whether the fitted quantity is a density, an expected count, or a tank signal. A plausible implication is that cross-experiment MLDF comparisons are most stable when they are expressed through explicitly defined reference-distance observables rather than through nominally similar but detector-specific normalizations.

4. Universality, scaling, and asymmetry

One strand of MLDF research argues that the muon profile obeys a one-parameter scaling law over a broad radial interval. For vertical showers at sea level with $C=0.28$5 MeV, the scaling description for muons is reported to be accurate over

$C=0.28$6

with relative uncertainties mostly within $C=0.28$7 across $C=0.28$8–$C=0.28$9 eV, proton and iron primaries, and both EPOS LHC and QGSJet-II-04. Within this picture, $\alpha=1.2$0 and $\alpha=1.2$1 are strongly anticorrelated event by event, and the rate of change $\alpha=1.2$2 is practically the same for EPOS LHC and QGSJet-II-04. By contrast, a single-parameter scaling of the total charged-particle LDF fails in the radial region where electrons and muons contribute comparably, with uncertainties up to $\alpha=1.2$3 (Raikin et al., 2018).

That result directly addresses a common simplification in surface-array work: the charged-particle LDF is not, in general, an adequate substitute for separate electromagnetic and muonic lateral distributions. This suggests that arrays with shielding, segmentation, or hybrid measurements have a structural advantage for composition analyses, because they can constrain the MLDF without forcing the mixed $\alpha=1.2$4 region into a single radial template.

For non-vertical showers, even a symmetric muon LDF in the shower plane does not remain symmetric in the ground plane. A dedicated analytical treatment finds that lateral and polar muon densities are asymmetric on the ground, that concentric shower-plane circles map to shifted ground-plane ellipses, and that the shift can be characterized by a gap length

$\alpha=1.2$5

The corresponding asymmetric ground-plane muon LDF replaces the radial coordinate by an effective ordinate $\alpha=1.2$6 and yields fits with mean absolute percentage error up to $\alpha=1.2$7 for muon lateral density distributions and up to $\alpha=1.2$8 for muon polar distributions at $\alpha=1.2$9 m with $\mu(r)$00 m (Basak et al., 2023). In this formulation, symmetric LDFs such as NKG-like expressions remain useful as shower-plane baselines, but they are inadequate for direct reconstruction in the ground plane when zenith-angle asymmetry is appreciable.

5. Large-radius structure and the hard muon component

At sufficiently large lateral separations, the measured muon distribution can cease to be governed solely by soft shower physics. IceCube measured laterally separated muons out to approximately $\mu(r)$01 m, a factor of $\mu(r)$02 larger than previous measurements by MACRO, and treated the observable as a separation distribution rather than a per-shower density. Under azimuthal symmetry one may relate it to a surface density through $\mu(r)$03, but the underlying measurement integrates over primary energies, zenith angles, and selection criteria rather than reconstructing a normalized MLDF for each shower (Collaboration et al., 2012).

The corrected separation spectrum exhibits a transition from an exponential soft component to a power-law tail. IceCube fitted it with

$\mu(r)$04

obtaining a transition scale around $\mu(r)$05 m, tail exponent $\mu(r)$06, and $\mu(r)$07, whereas a pure exponential fit gave $\mu(r)$08 (Collaboration et al., 2012). Through the small-angle relation between lateral separation and parent transverse momentum, this tail is interpreted as the manifestation of a hard component associated with perturbative-QCD production at $\mu(r)$09 GeV/c.

This large-radius hard component should be distinguished from the low-energy muon MLDFs measured by surface arrays such as IceTop, Yakutsk, or SUGAR. IceTop’s large-distance muon densities between roughly $\mu(r)$10 and $\mu(r)$11 m remain well described by a Greisen-like form anchored at $\mu(r)$12 m and exhibit the empirical scaling $\mu(r)$13 for near-vertical showers (Gonzalez, 2015). The juxtaposition suggests that “large-radius muons” can mean two different regimes: the classical GeV-scale surface MLDF and the much rarer, TeV-scale, laterally separated component that directly probes high-$\mu(r)$14 hadron production.

6. Model comparisons, shape discrepancies, and the muon puzzle

Several MLDF studies now emphasize that disagreements with simulations are not limited to the total muon count; the radial shape itself can differ. In the SUGAR reanalysis over $\mu(r)$15–$\mu(r)$16 eV, the normalized muon density is higher than QGSJET-II-04 and EPOS-LHC predictions at small $\mu(r)$17 and lower at larger $\mu(r)$18, so that the observed MLDF falls faster with increasing core distance than in simulation. The quoted goodness-of-fit values are $\mu(r)$19 and $\mu(r)$20 for EPOS-LHC proton and iron, and $\mu(r)$21 and $\mu(r)$22 for QGSJET-II-04 proton and iron; the probability that the observed trend arises by chance is $\mu(r)$23 (Kalmykov et al., 2022). This indicates that MLDF shape mismodeling can contaminate any inference that relies on muon densities measured over a restricted radial range.

Yakutsk reaches a different conclusion at higher energies. For $\mu(r)$24 eV, the measured surface and muon responses at $\mu(r)$25 m, when plotted as density divided by nominal $\mu(r)$26, are both lower than proton-model predictions, with the muon densities described as “significantly lower.” The authors argue that reducing the proportionality constant in the Yakutsk energy calibration by approximately $\mu(r)$27 for QGSJET, $\mu(r)$28 for QGSJET-II, or $\mu(r)$29 for EPOS-LHC would bring both SD and MD responses into agreement with simulations within experimental uncertainties at $\mu(r)$30 eV. In the same study, the derived $\mu(r)$31 values remain near zero within uncertainties, for example $\mu(r)$32 at $\mu(r)$33 under QGSJET, implying a light or light-mixed composition rather than an iron-dominated one (Glushkov et al., 2023).

Yakutsk’s dedicated $\mu(r)$34 eV analysis adds further nuance. Near-vertical bins are consistent with proton predictions within uncertainties, but more inclined bins show a muon deficit relative to proton simulations, corresponding to roughly $\mu(r)$35–$\mu(r)$36 lower measured densities at $\mu(r)$37 m and negative values of the proton-to-iron interpolation variable $\mu(r)$38. In a two-component interpretation, the mean weights are reported as $\mu(r)$39 and $\mu(r)$40 for a proton+$\mu(r)$41 mixture (Glushkov et al., 2023). This does not reproduce the standard “muon excess” narrative; instead it points to a detector- and geometry-dependent pattern of agreement and deficit.

The cross-experiment tension becomes sharp in the Yakutsk–Auger comparison. Yakutsk converts its $\mu(r)$42 values to $\mu(r)$43 using the Auger attenuation function

$\mu(r)$44

with

$\mu(r)$45

and finds agreement with proton-model simulations at $\mu(r)$46 eV, plus indications of a mixed composition at lower energies. The same work argues that if Auger’s energy scale were renormalized upward by a factor $\mu(r)$47, the Auger muon densities normalized by $\mu(r)$48 would drop by $\mu(r)$49, alleviating the discrepancy between direct muon measurements and optical composition indicators (Glushkov et al., 2023). The resulting “muon puzzle” is therefore not only a question of hadronic modeling; it is also a question of how MLDF observables are normalized, which reference distance is used, and how detector response and energy scale are harmonized across experiments.