Muon Production Depth: Concepts & Applications

• Muon production depth is the longitudinal distribution (measured in g/cm²) indicating where muons are generated in air showers and underground settings.
• It is reconstructed through timing and geometric analyses that convert detector signals and atmospheric profiles into a detailed production profile.
• Understanding muon production depth aids in evaluating cosmic ray composition, hadronic interaction models, and influences underground muon and neutrino detection.

Searching arXiv for recent and foundational papers on muon production depth to ground the article. Muon production depth denotes the longitudinal distribution of locations at which muons are produced, but the precise meaning depends on context. In extensive air showers, it is the distribution along the shower axis of the atmospheric depth $X$ at which muons are produced, usually expressed in $\mathrm{g\,cm^{-2}}$, and summarized by observables such as the depth of maximum muon production, $X_\mu^{\max}$ (Collica, 2016). In underground and neutrino-telescope applications, related constructions describe either the slant-depth profile of energetic muon production in the atmosphere (Gaisser et al., 2021, Verpoest et al., 2021), the stopping-depth distribution of cosmic-ray muons in rock or water near the Earth’s surface (Guo, 2018), or the depth distribution of muons generated by $\nu_\mu$ charged-current interactions along a chord through the Earth (Takahashi et al., 2010). The shared idea is a mapping from production or stopping location to an experimentally relevant longitudinal depth coordinate, but the physics, observables, and reconstruction strategies differ substantially.

1. Terminology, coordinates, and distinct usages

The standard air-shower definition is the muon production depth distribution $P(X_\mu)\equiv dN_\mu/dX$, where $X$ is slant depth along the shower axis and the reconstructed distribution is operationally the apparent distribution of production depths for those muons that survive to ground (Collaboration, 2014, Collica, 2016). The corresponding peak position,

$X_\mu^{\max}=\arg\max_X P(X_\mu),$

is the principal summary observable in Pierre Auger analyses (Collaboration, 2014).

A separate usage appears in atmospheric-bundle studies for underground detectors, where slant depth $X$ is again the independent variable, but the object of interest is a depth-resolved production spectrum of muons above threshold,

$\frac{dN_\mu}{dX}(>E_\mu,E_0,A,\theta,X),$

parameterized in terms of primary energy $E_0$, mass $\mathrm{g\,cm^{-2}}$0, zenith angle $\mathrm{g\,cm^{-2}}$1, threshold $\mathrm{g\,cm^{-2}}$2, and atmospheric state (Gaisser et al., 2021, Verpoest et al., 2021). In this setting, $\mathrm{g\,cm^{-2}}$3 encodes where parent mesons decay into muons along the shower trajectory rather than the apparent production distribution of ground-reaching muons reconstructed event by event.

A further distinction is required for low-energy neutrino backgrounds from stopped cosmic-ray muons. There, the relevant longitudinal variable is not atmospheric production depth but the stopping depth $\mathrm{g\,cm^{-2}}$4 beneath the Earth’s surface in rock or seawater, obtained from sea-level muon intensities and continuous slowing down approximation ranges (Guo, 2018). The paper explicitly notes that, in cosmic-ray studies, “muon production depth” often refers to where in the atmosphere muons are produced, whereas its focus is the stopping-depth distribution in the Earth’s surface layers (Guo, 2018).

For upward-going neutrino-induced muons, the depth variable is the location along the neutrino chord where a $\mathrm{g\,cm^{-2}}$5 charged-current interaction produces a muon. In that case column depth is

$\mathrm{g\,cm^{-2}}$6

with $\mathrm{g\,cm^{-2}}$7 the Earth density along the path (Takahashi et al., 2010). The same phrase therefore spans at least four related but non-identical observables.

Context	Depth variable	Primary observable
EAS composition studies	Atmospheric slant depth $\mathrm{g\,cm^{-2}}$8	$\mathrm{g\,cm^{-2}}$9, $X_\mu^{\max}$0
Underground bundle studies	Atmospheric slant depth $X_\mu^{\max}$1	$X_\mu^{\max}$2
Stopped surface muons	Stopping depth $X_\mu^{\max}$3 in rock/water	$X_\mu^{\max}$4
Neutrino-induced muons in Earth	Column depth or distance along chord	Interaction-depth distribution

This multiplicity of meanings explains a common misconception: the same phrase may refer either to an atmospheric production profile or to a terrestrial stopping or interaction-depth distribution. The underlying variable is always longitudinal depth, but the physical process being localized is different (Guo, 2018, Takahashi et al., 2010).

2. Reconstruction from air-shower timing

The Pierre Auger Observatory reconstructs atmospheric muon production depths from the arrival-time structure recorded by water-Cherenkov surface detectors (Collica, 2016). The method assumes that muons travel on straight lines from their production points to the ground, that ultra-relativistic muons have $X_\mu^{\max}$5 but not exactly $X_\mu^{\max}$6, and that additional effects such as pion-decay geometry, geomagnetic bending, and elastic scattering are either corrected or reduced by selection cuts (Collica, 2016).

For a station at lateral distance $X_\mu^{\max}$7 from the shower axis, with measured delay $X_\mu^{\max}$8 relative to a plane-front reference, the total delay is decomposed as

$X_\mu^{\max}$9

where the geometric delay is

$\nu_\mu$0

and the kinematic delay is

$\nu_\mu$1

Defining $\nu_\mu$2, the geometric inversion gives

$\nu_\mu$3

Auger’s operational form includes the inclined-shower correction $\nu_\mu$4 and a pion-decay offset $\nu_\mu$5, yielding

$\nu_\mu$6

with $\nu_\mu$7 after subtraction of the average kinematic delay (Collica, 2016, Collaboration, 2014).

The conversion from production height to slant depth uses the atmospheric density profile:

$\nu_\mu$8

For uncertainty propagation Auger uses an exponential density model, $\nu_\mu$9, which gives the timing-induced depth uncertainty

$P(X_\mu)\equiv dN_\mu/dX$0

This expression makes explicit why near-core stations are problematic: the uncertainty scales strongly with decreasing $P(X_\mu)\equiv dN_\mu/dX$1 (Collica, 2016, Collaboration, 2014).

The practical reconstruction chain is built from selected surface-detector traces. In the later Auger analysis, events are taken in the zenith range $P(X_\mu)\equiv dN_\mu/dX$2–$P(X_\mu)\equiv dN_\mu/dX$3, energies $P(X_\mu)\equiv dN_\mu/dX$4–$P(X_\mu)\equiv dN_\mu/dX$5, and station-core distances $P(X_\mu)\equiv dN_\mu/dX$6; the lower bound ensures the kinematic delay is $P(X_\mu)\equiv dN_\mu/dX$7 of the geometric delay (Collica, 2016). In the earlier Auger study the working range is $P(X_\mu)\equiv dN_\mu/dX$8, $P(X_\mu)\equiv dN_\mu/dX$9, with a fixed lower cut $X$0 and a waveform threshold $X$1 of the peak to retain bins with muon fractions $X$2 (Collaboration, 2014). The apparent MPD profile is then fit, in the 2016 analysis, with the Universal Shower Profile, with $X$3 fixed as a function of zenith angle and $X$4 as the key fit parameter (Collica, 2016). The earlier analysis instead used a Gaisser–Hillas fit with $X$5 fixed because limited muon sampling made a four-parameter fit unstable (Collaboration, 2014).

Detector-response corrections are integral to the method. The Surface Detector samples at $X$6, with GPS timing accuracy of about $X$7 (Collica, 2016). A global time offset $X$8 is subtracted in the 2016 Auger analysis to compensate average smearing from light propagation and electronics (Collica, 2016), whereas the 2014 reconstruction used a station-level timing offset $X$9 for the 15% threshold (Collaboration, 2014).

3. Parameterized production profiles in the atmosphere

For underground and under-ice detectors, the problem is not event-by-event inversion of surface timing but forward parameterization of the atmospheric muon production profile above an energy threshold (Gaisser et al., 2021, Verpoest et al., 2021). Both 2021 papers generalize the Elbert formula by introducing a depth-differential profile based on the derivative of a Gaisser–Hillas function for the hadronic cascade combined with explicit meson decay–reinteraction competition.

The inclusive Elbert formula is

$X_\mu^{\max}=\arg\max_X P(X_\mu),$0

with $X_\mu^{\max}=\arg\max_X P(X_\mu),$1 for the TeV-scale fit and $X_\mu^{\max}=\arg\max_X P(X_\mu),$2 for the $X_\mu^{\max}=\arg\max_X P(X_\mu),$3 threshold fit (Gaisser et al., 2021). The depth-resolved generalization is written schematically as a derivative of a Gaisser–Hillas profile multiplied by a meson decay factor and a threshold factor (Gaisser et al., 2021). The closely related explicit parameterization of Verpoest and Gaisser writes the mean muon production profile as a derivative of Gaisser–Hillas for charged mesons, multiplied by pion and kaon decay kernels and the threshold suppression factor $X_\mu^{\max}=\arg\max_X P(X_\mu),$4 (Verpoest et al., 2021).

The key decay–reinteraction competition is governed by

$X_\mu^{\max}=\arg\max_X P(X_\mu),$5

with

$X_\mu^{\max}=\arg\max_X P(X_\mu),$6

and

$X_\mu^{\max}=\arg\max_X P(X_\mu),$7

The generalized formulation introduces the combined pion/kaon contribution

$X_\mu^{\max}=\arg\max_X P(X_\mu),$8

with $X_\mu^{\max}=\arg\max_X P(X_\mu),$9, $X$0, $X$1, and $X$2 (Gaisser et al., 2021). In the Verpoest–Gaisser parameterization, the same relative weights 0.92 and 0.08 appear directly in the profile formula (Verpoest et al., 2021).

The longitudinal-shape parameters are fitted primarily as functions of

$X$3

For the TeV-threshold parameterization, the fitted forms are

$X$4

with analogous forms for $X$5 and $X$6, using two regimes separated by $X$7 (Gaisser et al., 2021). The quoted TeV-scale parameters are: for $X$8, $X$9 and $\frac{dN_\mu}{dX}(>E_\mu,E_0,A,\theta,X),$0; for $\frac{dN_\mu}{dX}(>E_\mu,E_0,A,\theta,X),$1 in $\frac{dN_\mu}{dX}(>E_\mu,E_0,A,\theta,X),$2, $\frac{dN_\mu}{dX}(>E_\mu,E_0,A,\theta,X),$3 and $\frac{dN_\mu}{dX}(>E_\mu,E_0,A,\theta,X),$4; for $\frac{dN_\mu}{dX}(>E_\mu,E_0,A,\theta,X),$5, $\frac{dN_\mu}{dX}(>E_\mu,E_0,A,\theta,X),$6 and $\frac{dN_\mu}{dX}(>E_\mu,E_0,A,\theta,X),$7; for $\frac{dN_\mu}{dX}(>E_\mu,E_0,A,\theta,X),$8, $\frac{dN_\mu}{dX}(>E_\mu,E_0,A,\theta,X),$9 and $E_0$0 (Gaisser et al., 2021). The low-threshold ($E_0$1) parameter set is likewise given explicitly in that paper (Gaisser et al., 2021).

A central implication is that atmospheric temperature enters the MPD not merely as an external nuisance but directly through $E_0$2, $E_0$3, and the depth–altitude mapping

$E_0$4

so warmer atmospheres favor meson decay at higher altitude, increase single-muon rates, and widen multi-muon bundles (Gaisser et al., 2021, Verpoest et al., 2021).

4. Composition sensitivity and hadronic-model discrimination

In ultra-high-energy cosmic-ray studies, $E_0$5 is sensitive both to primary mass and to hadronic interaction properties (Collica, 2016, Collaboration, 2014). Simulations at EeV energies give iron showers as developing earlier, with shallower $E_0$6 than proton showers by about $E_0$7 (Collica, 2016). This composition lever arm is large enough that the Auger systematic uncertainty of $E_0$8 across the full angular range is about $E_0$9 of the simulated proton–iron separation (Collica, 2016).

The model dependence is comparably important. The post-LHC models EPOS-LHC and QGSJetII-04 predict similar elongation rates but differ in absolute $\mathrm{g\,cm^{-2}}$00 by amounts comparable to the proton–iron separation (Collica, 2016). The earlier Auger study reports a muonic elongation rate

$\mathrm{g\,cm^{-2}}$01

which disfavors a pure proton $\mathrm{g\,cm^{-2}}$02 or pure iron $\mathrm{g\,cm^{-2}}$03 trend at about $\mathrm{g\,cm^{-2}}$04 and $\mathrm{g\,cm^{-2}}$05, respectively, under the quoted hadronic assumptions (Collaboration, 2014). The same study finds that QGSJetII-04 proton and iron predictions bracket the Auger $\mathrm{g\,cm^{-2}}$06 data, whereas EPOS-LHC iron predictions lie above the data (Collaboration, 2014).

The performance of the reconstruction is sufficient for these comparisons. In the 2016 Auger analysis, the resolution on $\mathrm{g\,cm^{-2}}$07 is about $\mathrm{g\,cm^{-2}}$08 at $\mathrm{g\,cm^{-2}}$09 improving to about $\mathrm{g\,cm^{-2}}$10 at $\mathrm{g\,cm^{-2}}$11 for $\mathrm{g\,cm^{-2}}$12–$\mathrm{g\,cm^{-2}}$13, and about $\mathrm{g\,cm^{-2}}$14 improving to about $\mathrm{g\,cm^{-2}}$15 for $\mathrm{g\,cm^{-2}}$16–$\mathrm{g\,cm^{-2}}$17 (Collica, 2016). The main contributors are discrete muon sampling at ground, detector time resolution, core/angle reconstruction, and a kinematic-delay parameterization whose contribution is $\mathrm{g\,cm^{-2}}$18 (Collica, 2016). In the 2014 study, the RMS of $\mathrm{g\,cm^{-2}}$19 improves from about $\mathrm{g\,cm^{-2}}$20 for proton and $\mathrm{g\,cm^{-2}}$21 for iron at $\mathrm{g\,cm^{-2}}$22 to about $\mathrm{g\,cm^{-2}}$23 at the highest energies (Collaboration, 2014).

A persistent interpretation issue is the relation between MPD constraints and the broader muon deficit problem. The Auger papers state that MPD is complementary to fluorescence-detector $\mathrm{g\,cm^{-2}}$24 because the electromagnetic profile is driven by early $\mathrm{g\,cm^{-2}}$25 production, whereas $\mathrm{g\,cm^{-2}}$26 tracks the charged hadron cascade down to pion critical energies (Collaboration, 2014, Collica, 2016). This suggests that mismatches between measured and simulated MPDs probe multiplicities, inelasticities, baryon production, and charge ratios in hadronic models rather than composition alone.

5. Temperature, primary mass, and bundle observables

The atmospheric MPD formalism developed for underground detectors connects longitudinal production depth directly to seasonal modulation, bundle multiplicity, and transverse size (Gaisser et al., 2021, Verpoest et al., 2021). The physical mechanism is the temperature dependence of the decay–reinteraction balance of charged pions and kaons. Since

$\mathrm{g\,cm^{-2}}$27

higher temperature corresponds to larger scale height and larger meson critical energies, increasing the probability that mesons decay before interacting (Gaisser et al., 2021). The seasonal-modulation relation is written as

$\mathrm{g\,cm^{-2}}$28

with

$\mathrm{g\,cm^{-2}}$29

where $\mathrm{g\,cm^{-2}}$30 is the production spectrum folded with the primary spectrum (Gaisser et al., 2021).

The generalized profile reproduces several detector-specific trends. For IceCube-like events with $\mathrm{g\,cm^{-2}}$31 and primary energies about $\mathrm{g\,cm^{-2}}$32–$\mathrm{g\,cm^{-2}}$33, the muon multiplicity varies by up to about $\mathrm{g\,cm^{-2}}$34 around the yearly mean, while the bundle radius varies seasonally by about $\mathrm{g\,cm^{-2}}$35, with a summer maximum due to higher production altitude (Gaisser et al., 2021). In the South Pole vertical-shower use case of Verpoest and Gaisser, for $\mathrm{g\,cm^{-2}}$36 showers and $\mathrm{g\,cm^{-2}}$37, the predicted seasonal variation of multiplicity is about $\mathrm{g\,cm^{-2}}$38 and the estimated transverse size varies by roughly $\mathrm{g\,cm^{-2}}$39 around the average (Verpoest et al., 2021).

The link between MPD and lateral spread is expressed geometrically through

$\mathrm{g\,cm^{-2}}$40

with a transverse momentum distribution

$\mathrm{g\,cm^{-2}}$41

and $\mathrm{g\,cm^{-2}}$42 (Gaisser et al., 2021). Averaging over the production profile yields a mean bundle radius,

$\mathrm{g\,cm^{-2}}$43

This formalism is used to show that heavy primaries develop higher in the atmosphere and therefore produce both larger muon multiplicities and wider bundles (Gaisser et al., 2021).

The composition dependence is quantitatively significant. For IceCube-like thresholds, the calculation indicates that iron-induced bundles are about $\mathrm{g\,cm^{-2}}$44 larger than proton-induced bundles and have about $\mathrm{g\,cm^{-2}}$45 the muon multiplicity (Gaisser et al., 2021). In compact underground detectors, the same altitude effect explains why multiple-muon rates can anti-correlate with temperature even while inclusive single-muon rates correlate positively: warmer atmospheres shift production higher, widening bundles and reducing the chance that closely separated muons jointly trigger a compact detector (Gaisser et al., 2021). This interpretation is explicitly advanced for MINOS and NOvA, although the paper also states that full quantitative agreement requires more detailed simulations of acceptance, propagation, and realistic bundle geometry (Gaisser et al., 2021).

6. Stopping depth in rock and water, and interaction depth through the Earth

A distinct terrestrial variant of the problem concerns low-energy neutrinos from cosmic-ray muons that stop near the Earth’s surface (Guo, 2018). The starting point is the sea-level muon intensity

$\mathrm{g\,cm^{-2}}$46

with

$\mathrm{g\,cm^{-2}}$47

where

$\mathrm{g\,cm^{-2}}$48

For $\mathrm{g\,cm^{-2}}$49, a flat approximation is used,

$\mathrm{g\,cm^{-2}}$50

These reproduce the standard total sea-level flux

$\mathrm{g\,cm^{-2}}$51

with muons above $\mathrm{g\,cm^{-2}}$52 contributing $\mathrm{g\,cm^{-2}}$53 (Guo, 2018).

Using tabulated CSDA ranges for standard rock and water, the stopping-depth mapping is

$\mathrm{g\,cm^{-2}}$54

and the stopping-depth rate density is constructed as

$\mathrm{g\,cm^{-2}}$55

Most stopped muons have shallow depths $\mathrm{g\,cm^{-2}}$56 in both standard rock and water (Guo, 2018). For stopped $\mathrm{g\,cm^{-2}}$57, the decay and nuclear capture probabilities depend on material composition: in upper continental crust, $\mathrm{g\,cm^{-2}}$58 and $\mathrm{g\,cm^{-2}}$59, whereas in water, treated as an oxygen target, $\mathrm{g\,cm^{-2}}$60 and $\mathrm{g\,cm^{-2}}$61 (Guo, 2018). This difference feeds directly into the low-energy neutrino yields, particularly the $\mathrm{g\,cm^{-2}}$62 component.

The detector-depth-dependent neutrino flux from these stopped muons is

$\mathrm{g\,cm^{-2}}$63

with

$\mathrm{g\,cm^{-2}}$64

and the approximation

$\mathrm{g\,cm^{-2}}$65

for $\mathrm{g\,cm^{-2}}$66 in meters (Guo, 2018). For $\mathrm{g\,cm^{-2}}$67 and $\mathrm{g\,cm^{-2}}$68, the resulting $\mathrm{g\,cm^{-2}}$69, $\mathrm{g\,cm^{-2}}$70, $\mathrm{g\,cm^{-2}}$71, and $\mathrm{g\,cm^{-2}}$72 fluxes are averagely $\mathrm{g\,cm^{-2}}$73, $\mathrm{g\,cm^{-2}}$74, $\mathrm{g\,cm^{-2}}$75, and $\mathrm{g\,cm^{-2}}$76 of the corresponding atmospheric neutrino fluxes, and these results increase by a factor of 1.4 if $\mathrm{g\,cm^{-2}}$77 (Guo, 2018). Most neutrinos come from within $\mathrm{g\,cm^{-2}}$78 and from near-horizontal directions (Guo, 2018).

At much higher energies, neutrino-induced muon production through the Earth defines yet another depth distribution (Takahashi et al., 2010). The chord length for an upward-going neutrino is

$\mathrm{g\,cm^{-2}}$79

and the local interaction probability depends on the Earth density profile through

$\mathrm{g\,cm^{-2}}$80

The charged-current muon-production distribution is influenced by layered Earth densities and by neutral-current regeneration of the parent neutrino (Takahashi et al., 2010). For long chords, the interaction-depth distributions exhibit abrupt changes near about $\mathrm{g\,cm^{-2}}$81 and about $\mathrm{g\,cm^{-2}}$82 from the detector, corresponding to entering and exiting the core (Takahashi et al., 2010). The paper emphasizes that, at $\mathrm{g\,cm^{-2}}$83, “all neutrino events which are produced within 100 km from the detector are produced via neutral currents,” meaning that near-detector muon production is dominated by sequences with one or more prior neutral-current interactions before the final charged-current event (Takahashi et al., 2010). This usage of production depth therefore tracks the longitudinal structure of neutrino interactions in the Earth rather than atmospheric shower development.

7. Limitations, ambiguities, and interpretive boundaries

Several limitations recur across the literature. In Auger reconstructions, the per-muon kinematic delay is not measured and must be parameterized from simulations; the method therefore relies on cuts that keep the kinematic component subdominant compared with the geometric delay (Collica, 2016, Collaboration, 2014). Electromagnetic contamination must be suppressed either by timing-structure algorithms or threshold cuts, and detector-response corrections such as the global time shift are essential to remove systematic bias (Collica, 2016, Collaboration, 2014). The quoted Auger systematics incorporate atmospheric variability, reconstruction biases, and modeling dependence, but the absolute interpretation of $\mathrm{g\,cm^{-2}}$84 remains model dependent (Collica, 2016).

In the atmospheric parameterizations for underground detectors, the dependence on $\mathrm{g\,cm^{-2}}$85 is stated to be approximate, with residual $\mathrm{g\,cm^{-2}}$86 dependence still visible (Gaisser et al., 2021, Verpoest et al., 2021). The papers recommend optimizing the fitted parameter sets for the target energy range and detector conditions (Gaisser et al., 2021). The geometric treatment of lateral spread captures the altitude effect but omits or only approximately treats geomagnetic bending and multiple Coulomb scattering in some applications (Gaisser et al., 2021, Verpoest et al., 2021).

For stopped cosmic-ray muons in the Earth, the calculation uses tabulated CSDA ranges and therefore neglects stochastic fluctuations; the paper states that for shallow ranges this is adequate (Guo, 2018). The sea-level intensity model is based on the Reyna parameterization with a flat low-momentum piece, and uncertainties from alternative parameterizations are not explored there (Guo, 2018). The rock calculation assumes an upper continental crust composition, while water is treated as an oxygen target because $\mathrm{g\,cm^{-2}}$87 quickly migrates to $\mathrm{g\,cm^{-2}}$88 (Guo, 2018).

For neutrino-induced muons through the Earth, the key uncertainties arise from deep inelastic scattering cross sections at very small Bjorken-$\mathrm{g\,cm^{-2}}$89, the Earth density model, and muon energy-loss parameters (Takahashi et al., 2010). A plausible implication is that any interpretation of the interaction-depth distribution near a detector is inseparable from uncertainties in both neutrino cross sections and neutral-current regeneration.

Taken together, these caveats delimit the scope of the term. “Muon production depth” is not a single universal observable but a family of depth-localized descriptions adapted to different transport and detection problems. In atmospheric cosmic-ray physics it is primarily a probe of hadronic cascade development and composition through observables such as $\mathrm{g\,cm^{-2}}$90 and $\mathrm{g\,cm^{-2}}$91 (Collica, 2016, Collaboration, 2014). In underground muon and neutrino studies it becomes a forward-model ingredient for rates, bundle geometry, seasonal modulation, and low-energy or ultra-high-energy backgrounds (Gaisser et al., 2021, Verpoest et al., 2021, Guo, 2018, Takahashi et al., 2010). This suggests that the unity of the concept lies in its depth-coordinate formalism rather than in any single experimental implementation.