Honda and HKKM: Modeling the Atmospheric Neutrino Flux

Atmospheric neutrinos opened the era of neutrino oscillation. The 1998 Super-Kamiokande observation of muon-neutrino disappearance as a function of zenith angle established that neutrinos oscillate, and a generation of experiments — Super-K’s continued refinement of the result, IceCube’s high-energy precision measurement, and now the upcoming KM3NeT ORCA and Hyper-Kamiokande — continues to mine the atmospheric flux for oscillation and other physics. None of these analyses would be possible without an accurate model of the neutrino flux entering the detector before any oscillation has happened. The flux model is what you compare data against.

The standard parameterisation in use across the field is the HKKM model, named after Honda, Kajita, Kasahara and Midorikawa, who have produced and updated the calculation since the mid-1990s. The model takes the measured cosmic-ray spectrum, the known atmospheric profile, and a Monte Carlo treatment of hadronic interactions, and outputs the all-flavor neutrino flux at Earth as a function of energy, zenith angle, and geomagnetic location. Its tables are the input every atmospheric-neutrino analysis starts from.

This post is about how the atmospheric neutrino flux is constructed, what physical inputs the model needs, and where its remaining uncertainties live.

From cosmic rays to neutrinos

A 100 GeV proton arriving at the top of the atmosphere has a mean free path of about 80 g/cm² before its first inelastic interaction. The Earth’s atmosphere presents 1030 g/cm² to a vertical particle, so almost every primary cosmic ray initiates a shower somewhere in the upper troposphere or stratosphere. The first interaction produces a handful of charged and neutral pions, kaons, and a few other secondaries.

The charged pions and kaons live long enough to either decay in flight or interact again with the atmosphere depending on their energy. At low energies most decay before re-interacting; at higher energies the interaction probability dominates and the meson loses energy. The decay chain is

${\pi }^{+}\to {\mu }^{+}+{\nu }_{\mu }{K}^{+}\to {\mu }^{+}+{\nu }_{\mu }{K}^{+}\to {\pi }^0+e^{+}+{\nu }_e$

with similar chains for the negative counterparts. The muon then decays:

${\mu }^{+}\to e^{+}+{\nu }_e+{\bar{\nu }}_{\mu }$

Each parent meson produces one muon neutrino in its initial decay; each muon, if it decays before reaching the ground, produces one more muon neutrino and one electron neutrino. The result at the bottom of the atmosphere is a roughly 2:1 ratio of muon-flavor to electron-flavor neutrinos at energies below about 1 GeV, transitioning to a higher ratio as the muons increasingly reach the ground without decaying at higher energies.

What goes into the model

Computing the flux that emerges from this chain requires several ingredients.

The first is the primary cosmic-ray spectrum. Direct measurements by AMS-02 on the International Space Station, PAMELA, BESS, CREAM and others have measured the proton and helium spectra over many decades in energy with high precision. The combined data anchors the input at the 5-10% level for energies up to about 100 GeV and somewhat worse at higher energies, where the data come increasingly from indirect measurements via cosmic-ray-air-shower observations on the ground.

The second is the geomagnetic-field cutoff. Low-rigidity primary cosmic rays are deflected by the Earth’s magnetic field and only reach the top of the atmosphere if their rigidity exceeds the local cutoff. The cutoff varies with geomagnetic latitude — high at the equator, low at the poles — and modulates the low-energy end of the primary spectrum location-by-location. Honda et al. compute the cutoff using a detailed magnetic-field model.

The third is the atmospheric density profile. The shower development depends sensitively on how matter is distributed with altitude. Standard atmospheric models such as US Standard Atmosphere 1976 or seasonal-averaged profiles from NRLMSISE-00 are used; the variation with season and latitude is included as a separate systematic.

The fourth is the hadronic interaction model. The pion and kaon production cross-sections in proton-nitrogen and proton-oxygen collisions at the relevant energies are constrained by accelerator measurements but extrapolated to the highest atmospheric energies. The choice of hadronic-interaction model (DPMJET-III, SIBYLL, QGSJET) introduces a 10-15% systematic in the resulting neutrino flux above about 100 GeV.

The fifth is the solar-cycle modulation. Low-energy primary cosmic rays are deflected by the solar wind, with the suppression varying over the 11-year solar cycle. The effect is significant below about 10 GeV and is parameterised by the solar-modulation potential, which is tabulated as a function of time.

The output and how it is used

The HKKM tables give the neutrino flux $d\Phi /dE/d\cos \theta$ for each of the four signed neutrino species (${\nu }_e,{\bar{\nu }}_e,{\nu }_{\mu },{\bar{\nu }}_{\mu }$), at each of several major underground sites — Super-K’s Kamioka, IceCube’s South Pole, INO’s location in India, Gran Sasso, and SNOLAB — and as a function of solar modulation. Energy coverage spans from 100 MeV to 10 TeV; zenith angle is covered from straight-down to straight-up.

A typical Super-K analysis uses the HKKM tables as the input flux, propagates it through the standard PMNS three-flavor oscillation framework with the atmospheric and solar mass-squared differences and mixing angles, applies the cross-section model for each interaction channel, and folds in the detector response to produce predicted event distributions. The data are then compared and fit jointly for oscillation parameters and nuisance parameters that encode the flux uncertainties.

The model is updated periodically as new cosmic-ray data become available. The most recent major HKKM release in 2015 included the AMS-02 proton spectrum and refined the geomagnetic cutoff treatment. A 2022 update used improved hadronic-interaction models constrained by NA61/SHINE pion-production data.

Where the uncertainty sits

The HKKM uncertainty budget breaks down roughly as follows.

At low energies (sub-GeV), the dominant uncertainty is in the geomagnetic-cutoff calculation and the solar modulation, both of which affect the absolute normalisation by a few per cent. The hadronic model uncertainty is small here because the relevant proton-nitrogen cross-sections are well measured at low energies.

At GeV-scale energies, where Super-K does most of its oscillation work, the absolute flux uncertainty is around 7-10% and the flux ratio uncertainties (ν_μ/ν_e and ν/ν̄) are at the 2-3% level. The ratios cancel many common uncertainties and are what Super-K’s oscillation extraction actually relies on most heavily.

At TeV energies and above, where IceCube operates, the dominant uncertainty becomes the hadronic-interaction model — specifically the kaon-versus-pion production ratio, since kaons increasingly dominate the high-energy flux because they are heavier and decay faster than pions before re-interacting. The uncertainty rises to 15-25% on the absolute flux. IceCube’s atmospheric muon-neutrino measurement consequently has substantial systematic error from flux modelling at the highest energies.

A separate, sometimes-relevant uncertainty is the prompt atmospheric flux. Charmed mesons produced in the same shower decay essentially instantly, producing neutrinos that have a much flatter energy spectrum than the conventional pion-and-kaon flux. The prompt component is expected to dominate above about 100 TeV but has not been definitively measured yet; theoretical predictions span a factor of three or so. IceCube and KM3NeT analyses include the prompt flux as a free or partially-constrained component.

Why the model still matters

Modern oscillation analyses largely use the HKKM model in a constrained-fit framework, in which the flux uncertainties are treated as nuisance parameters with priors set by the model’s quoted uncertainties, and the data are allowed to pull the flux parameters along with the oscillation parameters. The data-driven approach reduces the impact of any specific flux mismodelling — if the data point at a slightly different flux shape, the fit can accommodate it.

But this approach degrades into pure measurement when the flux uncertainty is larger than the experimental error. Super-K’s current analyses are at the per-cent level on Δm²_32 and θ_23, comparable to the HKKM flux-ratio uncertainties, so the model is the practical limit. Further reducing the flux uncertainty — through improved primary-cosmic-ray measurements, refined hadronic-interaction modelling constrained by accelerator data, and direct measurements of relevant muon and kaon production cross-sections at the LHC fixed-target programs — is one of the active areas of work for the next generation.

For IceCube and the next-generation underwater telescopes, the same is true at TeV-PeV energies, where the prompt flux and the kaon/pion ratio are the dominant systematics. Direct measurements at LHCb of forward-region charm production are part of how this will be tightened.

Summary

The HKKM model — Honda, Kajita, Kasahara and Midorikawa — is the standard parameterisation of the atmospheric neutrino flux as a function of energy, zenith angle and geomagnetic location. It takes as inputs the directly measured cosmic-ray primary spectrum, geomagnetic-cutoff calculations, atmospheric profiles, hadronic-interaction Monte Carlo, and solar-modulation tabulations, and outputs the all-flavor neutrino flux at each major underground site. The model uncertainties — around 7-10% on absolute flux at GeV energies, rising to 15-25% at TeV scales, and a few per cent on flux ratios — are now at the level of the experimental errors in Super-K and IceCube’s precision atmospheric measurements. Improving the model through better primary-cosmic-ray data, refined hadronic-interaction modelling, and direct accelerator measurements is a quiet but essential part of pushing the precision of atmospheric-neutrino oscillation results further in the coming decade.