Local Muon Density Spectra

• Local Muon Density Spectra are the distributions of muon counts per unit area in air showers, measured under fixed geometric and threshold conditions.
• They convert detector signals—such as hit multiplicities and charge distributions—into muon densities using simulation-based calibration and correction factors.
• Analyses of LMDS reveal key spectral features like the second knee and offer composition diagnostics by comparing observations with predictions from hadronic interaction models.

Searching arXiv for the cited LMDS and related muon-density papers to ground the article in current arXiv records. Local Muon Density Spectra (LMDS) are distributions of air-shower events as functions of the locally measured muon density at the observation level or, equivalently, on a plane perpendicular to the shower axis. In extensive air shower work, the observable is always conditioned by geometry and threshold: zenith angle, lateral distance from the shower core when relevant, and the detector’s effective muon-energy threshold. In that form, LMDS provide a detector-level bridge between measured muon content, the primary cosmic-ray spectrum, mass composition, and hadronic-interaction modeling. The concept is realized in different but compatible ways in the NEVOD complex, where muon-bundle multiplicities are converted to local densities, and in IceTop, where GeV muon densities are extracted from tank charge distributions at fixed reference radii (Kokoulin et al., 2017, Soldin, 2021, Cazon, 2020).

1. Definition and mathematical form

The local muon density $\rho_\mu$ is the number of muons per unit area in a small region around the observation point for a given shower and arrival direction. In near-vertical IceTop analyses, the corresponding quantity is the number of GeV-scale muons per unit area at lateral distance $r$ from the shower axis, measured on a plane perpendicular to the shower axis. If $\lambda(r)$ is the expected muon count in a tank at distance $r$ and $A_{\mathrm{eff}}$ is the tank’s projected effective area, then

$$\lambda(r)=\rho_\mu(r)\,A_{\mathrm{eff}}, \qquad \rho_\mu(r)=\frac{\lambda(r)}{A_{\mathrm{eff}}}.$$

For muon bundles in NEVOD, an approximate multiplicity-to-density conversion is

$\rho_\mu \approx \frac{m}{A_{\text{eff}}(\theta)},$

with $m$ the reconstructed bundle multiplicity and $A_{\text{eff}}(\theta)$ the detector effective area for the zenith angle $\theta$ (Soldin, 2021, Kokoulin et al., 2017).

The differential LMDS in a zenith bin can be written as

$r$0

while a more general notation for the local muon density spectrum is $r$1, the differential event rate as a function of the local density under specified geometry and threshold. At fixed IceTop reference distance $r$2, one may also define an integral spectrum $r$3 and a differential spectrum $r$4 (Kokoulin et al., 2017, Cazon, 2020, Soldin, 2021).

The relation between LMDS and more familiar muon-size observables follows from the muon lateral distribution function (LDF). If the LDF factorizes as

$r$5

then at fixed $r$6, or within a fixed annulus where $r$7 varies weakly over the interval considered, one has $r$8, and the LMDS slope coincides with the muon-size spectrum slope. In KASCADE language, the truncated muon number is

$r$9

so spectra in $\lambda(r)$0 are integral transforms of LMDS over that radial band (Bijay et al., 2015).

2. Instrumental realizations and observation conditions

LMDS have been reconstructed with markedly different detector layouts. The NEVOD experimental complex uses the coordinate detector DECOR for inclined muon bundles and the Calibration Telescope System (CTS) for near-vertical bundles. IceTop uses water Cherenkov tanks on the surface of IceCube and derives local muon density from tank charge distributions in near-vertical air showers (Kokoulin et al., 2017, Soldin, 2021).

Instrument	Muon observable	Key conditions
DECOR	Bundle multiplicity converted to $\lambda(r)$1	Effective zenith angles $\lambda(r)$2, $\lambda(r)$3, $\lambda(r)$4; $\lambda(r)$5; about $\lambda(r)$6 h live time
CTS	Bottom-plane hit multiplicity converted to $\lambda(r)$7	Effective zenith angle $\lambda(r)$8; 3–40 hit counters; about $\lambda(r)$9 h live time
IceTop	$r$0 at fixed $r$1 from SLC charge fits	$r$2; $r$3 m and $r$4 m; about $r$5 live days

DECOR consists of eight supermodules placed in building galleries around the Cherenkov water detector. Each supermodule contains eight vertical planes of streamer-tube chambers, each plane having area $r$6, and the reconstructed muon-track angular accuracy is better than $r$7. The analysis uses three inclined zenith bins, $r$8–$r$9, $A_{\mathrm{eff}}$0–$A_{\mathrm{eff}}$1, and $A_{\mathrm{eff}}$2–$A_{\mathrm{eff}}$3, with effective zenith angles of $A_{\mathrm{eff}}$4, $A_{\mathrm{eff}}$5, and $A_{\mathrm{eff}}$6 respectively. The sample comprises about $A_{\mathrm{eff}}$7 events with muon bundle multiplicity $A_{\mathrm{eff}}$8, taken over two periods, 2002–2007 and 2012–2016, for a total live time of about $A_{\mathrm{eff}}$9 h (Kokoulin et al., 2017).

CTS consists of two planes of scintillation counters arranged in a chess order over $$\lambda(r)=\rho_\mu(r)\,A_{\mathrm{eff}}, \qquad \rho_\mu(r)=\frac{\lambda(r)}{A_{\mathrm{eff}}}.$$0, with $$\lambda(r)=\rho_\mu(r)\,A_{\mathrm{eff}}, \qquad \rho_\mu(r)=\frac{\lambda(r)}{A_{\mathrm{eff}}}.$$1 counters per plane and counter dimensions $$\lambda(r)=\rho_\mu(r)\,A_{\mathrm{eff}}, \qquad \rho_\mu(r)=\frac{\lambda(r)}{A_{\mathrm{eff}}}.$$2. For LMDS, only the bottom plane is used, below an $$\lambda(r)=\rho_\mu(r)\,A_{\mathrm{eff}}, \qquad \rho_\mu(r)=\frac{\lambda(r)}{A_{\mathrm{eff}}}.$$3 m water layer that suppresses the electromagnetic component. Events are recorded when at least two bottom-plane counters fire, and LMDS are reconstructed from events with 3–40 hit counters. The effective zenith angle is taken as $$\lambda(r)=\rho_\mu(r)\,A_{\mathrm{eff}}, \qquad \rho_\mu(r)=\frac{\lambda(r)}{A_{\mathrm{eff}}}.$$4, assuming an angular distribution of muon bundles proportional to $$\lambda(r)=\rho_\mu(r)\,A_{\mathrm{eff}}, \qquad \rho_\mu(r)=\frac{\lambda(r)}{A_{\mathrm{eff}}}.$$5 (Kokoulin et al., 2017).

IceTop is located at altitude $$\lambda(r)=\rho_\mu(r)\,A_{\mathrm{eff}}, \qquad \rho_\mu(r)=\frac{\lambda(r)}{A_{\mathrm{eff}}}.$$6 km and atmospheric depth $$\lambda(r)=\rho_\mu(r)\,A_{\mathrm{eff}}, \qquad \rho_\mu(r)=\frac{\lambda(r)}{A_{\mathrm{eff}}}.$$7, and the measurement is quoted at an atmospheric depth of about $$\lambda(r)=\rho_\mu(r)\,A_{\mathrm{eff}}, \qquad \rho_\mu(r)=\frac{\lambda(r)}{A_{\mathrm{eff}}}.$$8. The array comprises 81 stations on a triangular grid with about $$\lambda(r)=\rho_\mu(r)\,A_{\mathrm{eff}}, \qquad \rho_\mu(r)=\frac{\lambda(r)}{A_{\mathrm{eff}}}.$$9 m spacing; each station has two water Cherenkov tanks about $\rho_\mu \approx \frac{m}{A_{\text{eff}}(\theta)},$0 m apart. Signals are calibrated in units of Vertical Equivalent Muon (VEM). The analysis uses data from May 31, 2010 to May 2, 2013, corresponding to about $\rho_\mu \approx \frac{m}{A_{\text{eff}}(\theta)},$1 live days and more than $\rho_\mu \approx \frac{m}{A_{\text{eff}}(\theta)},$2 million selected events. The near-vertical requirement is $\rho_\mu \approx \frac{m}{A_{\text{eff}}(\theta)},$3, and the reconstructed-energy threshold is $\rho_\mu \approx \frac{m}{A_{\text{eff}}(\theta)},$4 PeV (Soldin, 2021).

3. Reconstruction of local muon density

In DECOR, muon bundles are identified through multiple reconstructed tracks consistent in direction and crossing DECOR planes. The conversion to local density is based on the effective area projected onto the plane perpendicular to the bundle direction, with acceptance and efficiency corrections from the established DECOR LMDS methodology. In this configuration, background from non-muon components is negligible for inclined bundles (Kokoulin et al., 2017).

In CTS, the bottom-plane hit multiplicity is converted to $\rho_\mu \approx \frac{m}{A_{\text{eff}}(\theta)},$5 using the geometric counter coverage and plane area. A detector-specific correction is required because electromagnetic cascades and penetrating hadrons generated in building materials and water cause a systematic overestimation of $\rho_\mu \approx \frac{m}{A_{\text{eff}}(\theta)},$6 by about $\rho_\mu \approx \frac{m}{A_{\text{eff}}(\theta)},$7 in raw data. That factor was evaluated and corrected with Geant4 simulations of the CTS response (Kokoulin et al., 2017).

IceTop reconstructs the shower geometry, core, and size parameter $\rho_\mu \approx \frac{m}{A_{\text{eff}}(\theta)},$8 with Hard Local Coincidence (HLC) signals, while muons are extracted mainly from Soft Local Coincidence (SLC) signals at large lateral distance, where the charge-versus-distance distribution exhibits the “Muon Thumb” near about $\rho_\mu \approx \frac{m}{A_{\text{eff}}(\theta)},$9 VEM. In each $m$0 bin, the SLC charge distribution is fitted with a log-likelihood multi-component model containing a muon response model, an electromagnetic signal model, and an accidental-background model. The muon response includes charge probability density functions for up to $m$1 simultaneous muons; the EM component is represented by empirical models EM1 and EM2; and accidental background is modeled with a Poisson process based on off-time windows. The fit returns the mean number of muons per tank, $m$2, and dividing by the tank projected area gives $m$3 (Soldin, 2021).

The underlying stochastic model uses the Poisson law

$m$4

with $m$5, so that

$m$6

IceTop denotes the raw reconstructed density by $m$7. A small Monte Carlo bias is corrected with an energy-dependent factor

$m$8

and the final density is

$m$9

The applied correction is the average of the proton and iron corrections, and half of the proton–iron difference is assigned as a systematic uncertainty. The reconstructed $A_{\text{eff}}(\theta)$0 curves are then interpolated to the reference distances $A_{\text{eff}}(\theta)$1 m and $A_{\text{eff}}(\theta)$2 m (Soldin, 2021).

4. Spectral parameterizations and mapping to primary energy

The LMDS in a zenith bin are commonly represented by a power law,

$A_{\text{eff}}(\theta)$3

and, where a knee-like steepening is observed, by a broken power law

$A_{\text{eff}}(\theta)$4

In the NEVOD analysis, the two fitted density intervals correspond to primary-energy ranges $A_{\text{eff}}(\theta)$5–$A_{\text{eff}}(\theta)$6 eV and $A_{\text{eff}}(\theta)$7 eV, with the knee position $A_{\text{eff}}(\theta)$8 corresponding to a knee energy $A_{\text{eff}}(\theta)$9 eV. The mapping from $\theta$0 and $\theta$1 to primary energy is obtained from CORSIKA simulations; the contextual relation

$\theta$2

is cited, but explicit values of $\theta$3 and $\theta$4 are not quoted (Kokoulin et al., 2017).

A broader scaling framework is given in the combined analysis of eight air-shower experiments. The total muon number obeys

$\theta$5

with $\theta$6. For fixed $\theta$7, $\theta$8, and threshold, the local density scales similarly:

$\theta$9

If the cosmic-ray differential flux is $r$00, then with $r$01 one obtains the approximate LMDS scaling

$r$02

This expresses how the event spectrum in primary energy is transferred into a spectrum in local muon density (Cazon, 2020).

The same mapping appears in size-spectrum analyses. If

$r$03

then

$r$04

Under a factorized muon LDF, LMDS at fixed $r$05 follow

$r$06

In this sense, LMDS inherit the same spectral information as muon-size spectra, while retaining explicit control of detector radius or fiducial band (Bijay et al., 2015).

For IceTop, the conversion from measured $r$07 to a local muon density spectrum at fixed $r$08 is written as

$r$09

and

$r$10

Above $r$11 PeV, the IceTop efficiency is close to $r$12 for all masses, so $r$13 in the considered range (Soldin, 2021).

5. Measured spectral features and the second knee

The clearest direct LMDS feature reported in the supplied literature is a steepening above primary energies of about $r$14 eV in NEVOD data. Both DECOR and CTS show an increase in the LMDS slope above that energy, interpreted as a “second knee” in the local muon density spectra (Kokoulin et al., 2017).

Data set	Below-knee slope	Above-knee slope
DECOR, combined inclined bins	$r$15	$r$16
CTS, near-vertical	$r$17	$r$18

For DECOR, the combined independent inclined bins give

$r$19

which corresponds to about $r$20 statistical significance. The individual DECOR fits are also reported: at $r$21, $r$22 in the $r$23–$r$24 eV interval; at $r$25, $r$26 below the knee and $r$27 above; at $r$28, $r$29 above the knee. For CTS, the corresponding slope change is

$r$30

or about $r$31. In both setups the knee energy is again about $r$32 eV, although the corresponding $r$33 values are indicated in the spectra rather than tabulated (Kokoulin et al., 2017).

IceTop probes a lower energy domain with fixed-radius local densities rather than explicit broken-power-law LMDS fits. It reports $r$34 at $r$35 m for $r$36–$r$37 PeV and at $r$38 m for $r$39–$r$40 PeV. The measured densities increase monotonically with energy at both reference radii. The internal consistency check is that densities are lower at $r$41 m than at $r$42 m for the same energy, while both radii show the same monotonic energy trend. The published results are graphical, with statistical error bars and bracketed systematic uncertainties; the paper does not provide numerical tables or explicit power-law fit coefficients for $r$43 (Soldin, 2021).

6. Composition sensitivity and hadronic-interaction tests

Muon-density observables are explicitly composition sensitive. In IceTop, this is encoded in the parameter

$r$44

where $r$45 is the measured density and $r$46 and $r$47 are the proton and iron predictions of a given hadronic model. Values near $r$48 are proton-like and values near $r$49 are iron-like. The measured $r$50 values are compared to expectations from the cosmic-ray flux models GSF, GST, and H3a (Soldin, 2021).

At IceTop energies, Sibyll 2.1 predictions agree with the measured $r$51 within uncertainties for physically reasonable cosmic-ray flux models, at least up to about $r$52 PeV. By contrast, EPOS-LHC and QGSJet-II.04 predict higher muon densities than Sibyll 2.1 and than the data over much of the $r$53–$r$54 PeV interval. At low energies below about $r$55 PeV, the post-LHC models imply an unrealistically light composition if one attempts to force agreement through composition alone. In the IceTop analysis, Sibyll 2.1 yields $r$56 values consistent with GSF, GST, and H3a within uncertainties, whereas EPOS-LHC and QGSJet-II.04 yield too-light compositions (Soldin, 2021).

The combined analysis introduces a universal reference scale,

$r$57

for detector-level muon density estimates. If simulations are perfect, this scale maps linearly from $r$58 for pure proton to $r$59 for pure iron, with

$r$60

To isolate model deviations from composition trends, the analysis fits

$r$61

where $r$62. The baseline slopes are $r$63 for EPOS-LHC and $r$64 for QGSJet-II.04, and robustness studies give $r$65 to $r$66 for EPOS-LHC and $r$67 to $r$68 for QGSJet-II.04. In all cases the slope differs from model expectations at more than $r$69. Above about $r$70 PeV, most experimental data show a muon excess relative to simulations, so the measured LMDS are shifted to higher local densities than the corresponding Monte Carlo predictions (Cazon, 2020).

This composition sensitivity also clarifies why LMDS and charged-particle spectra can lead to different inferences near a knee. In the Monte Carlo study of simultaneous charged and muon spectra, the mapping

$r$71

implies that the muon-spectrum or LMDS break tracks the primary-spectrum break in a manner that is less sensitive to abrupt composition changes than the charged-particle spectrum. In the simulated imposed-knee examples, the muon-spectrum break remains at $r$72–$r$73, whereas the charged-particle break can vary much more strongly, reaching $r$74 in a $r$75 transition (Bijay et al., 2015).

7. Uncertainties, cross-calibration, and interpretive limits

The systematic budget depends strongly on the detector implementation. In IceTop, four dominant sources are identified. The energy scale and resolution translate to about $r$76 uncertainty in $r$77; the EM signal model choice, assessed by comparing EM1 and EM2, induces up to about $r$78 uncertainty; the Monte Carlo correction factor depends on composition and on hadronic model choice; and detector-related effects such as snow attenuation and VEM calibration are folded into the energy and EM-model systematics. Statistical uncertainties are shown as error bars and systematic uncertainties as brackets around the points (Soldin, 2021).

In NEVOD, DECOR quotes a track-reconstruction accuracy better than $r$79, and acceptance and efficiency follow established procedures, but systematic uncertainties on the fitted $r$80 values are not itemized separately in the paper. For CTS, the principal explicit correction is the Geant4-derived factor of about $r$81 accounting for accompanying particles under building materials and water. The residual systematic effect of that correction on the LMDS slopes is not numerically quoted (Kokoulin et al., 2017).

Cross-experiment synthesis introduces an additional calibration layer. Because $r$82 with $r$83, even a $r$84 energy-scale offset induces about an $r$85 offset in muon density. The combined analysis therefore cross-calibrates energy scales using the isotropic cosmic-ray flux as reference. The residual uncertainty of the cross-calibrated energy scale is at least $r$86, implying a collective uncertainty of about $r$87 in the $r$88 scale. The same study emphasizes that measurements span a wide range of $r$89 values, lateral distances, and effective production-energy thresholds, with $r$90 spanning about $r$91–$r$92 GeV; extreme zenith angles such as $r$93 add complications from atmospheric-density gradients, $r$94 critical energies, and effective core distance for muon bundles (Cazon, 2020).

Several formal limitations remain explicit in the source literature. IceTop does not tabulate $r$95 values or analytic fit coefficients for $r$96, so exact numerical LMDS at $r$97 m and $r$98 m require digitizing the plotted points or using a supplemental data release if available. In the NEVOD second-knee analysis, the knee positions $r$99 are indicated graphically rather than tabulated, and the specific hadronic-interaction models used in the CORSIKA energy mapping are not stated. These omissions do not alter the reported qualitative results, but they constrain exact reproduction of LMDS parameterizations from the text alone (Soldin, 2021, Kokoulin et al., 2017).

Taken together, the supplied studies establish LMDS as a compact but information-rich muon observable. At detector level, LMDS encode the local muon content of showers across fixed geometry and threshold conditions; through simulation or fitted scaling relations they map onto primary energy; and through proton–iron bracketing they become composition diagnostics. The published measurements show both a second-knee steepening around $\lambda(r)$00 eV in NEVOD LMDS and, in the broader PeV-to-EeV comparison, an energy-growing muon excess above about $\lambda(r)$01 PeV relative to recent hadronic-interaction models (Kokoulin et al., 2017, Cazon, 2020).