QGSJet-II-04: Hadronic Interaction Model

• QGSJet-II-04 is a high-energy hadronic interaction model that employs a Reggeon Field Theory framework with multiple-Pomeron exchange for simulating extensive air showers from TeV to EeV energies.
• It integrates post-LHC retuning and enhanced Pomeron screening to accurately predict key observables such as muon densities, lateral distributions, and energy calibrations in cosmic-ray experiments.
• The model is widely implemented via Monte Carlo tools like CORSIKA, serving as a benchmark for detector-response simulations in observatories including IceCube, HAWC, and Yakutsk.

QGSJet-II-04 is a high-energy hadronic interaction model used primarily in Monte Carlo descriptions of extensive air showers (EAS), atmospheric muons, and related detector-response calculations. In the studies summarized here, it is implemented mainly through CORSIKA and supplies the high-energy hadron–air interaction kernel for analyses spanning roughly the TeV-to-EeV domain. Its standard characterization is a non-linear Reggeon Field Theory framework with multiple-Pomeron exchange, enhanced Pomeron screening, a semihard component for minijet production, and post-LHC retuning of cross sections and forward particle production. The model is therefore both a microscopic interaction generator and a practical inference engine: it is used to predict muon densities, attenuation lengths, multimuon multiplicity distributions, lateral distribution functions, and surface-detector energy calibrations (Fomin et al., 2016, Collaboration, 2024).

1. Theoretical basis and post-LHC retuning

QGSJET-II-04 is described in the cited studies as a Reggeon-Field-Theory or Gribov–Regge framework in which hadronic collisions are represented through multiple Pomeron exchanges in an eikonal formalism. Soft and semihard processes are treated in a unified way, non-linear screening is incorporated through enhanced Pomeron diagrams, and diffraction is included through the model’s standard treatment of diffractive channels. Several summaries further note a continuous transition from soft to hard interactions and explicit treatment of semihard minijet production (Fomin et al., 2016, Collaboration, 2024, Glushkov et al., 2024).

Representative parameterizations quoted in applications include an inelastic cross section scaling of the form

$$\sigma_{\mathrm{inel}}(s)\simeq \sigma_0\Bigl(\frac{s}{s_0}\Bigr)^\Delta,$$

with $\Delta\sim 0.1\mbox{–}0.15$ in the ALICE summary, and a proton–air energy dependence at PeV energies consistent with

$$\sigma_{p\text{-air}}(E)=\sigma_0\Bigl(\frac{E}{E_0}\Bigr)^\alpha, \qquad \sigma_0=265~\mathrm{mb},\quad E_0=0.2~\mathrm{TeV},\quad \alpha=0.079\pm0.010,$$

corresponding to an $8\%\pm1\%$ rise per decade and $$\sigma_{p\text{-air}}(1\,\mathrm{PeV})=350\pm15~\mathrm{mb}$$ in the Tien Shan analysis (Collaboration, 2024, Nesterova, 2015).

A central post-LHC modification emphasized in the ALICE study is the retuning to early LHC data, including total and elastic $pp$ cross sections and forward particle spectra. That retuning included a new pion-exchange mechanism in forward neutral hadron production, which increased the $\rho^0$ yield and hence the EAS muon content by $\approx20\%$ relative to previous QGSJET versions (Collaboration, 2024). The CTA comparison likewise characterizes QGSJET-II-04 as a retuned successor to QGSJET-II-03, with slightly softer forward spectra of leading baryons and neutral pions at TeV energies (Ohishi et al., 2021).

2. Computational implementation in air-shower simulations

In published applications, QGSJET-II-04 is most often interfaced through CORSIKA, with the low-energy hadronic sector supplied by FLUKA, GHEISHA, or UrQMD depending on the experiment. The thresholds summarized in the data vary by configuration: some studies use QGSJET-II-04 above $E_{\rm lab}\gtrsim80$ GeV, ALICE uses UrQMD below 100 GeV and QGSJET-II-04 above, and the KASCADE-Grande muon study uses FLUKA below 200 GeV and QGSJET-II-04 above (Glushkov et al., 2024, Collaboration, 2024, Fomin et al., 2016, Arteaga-Velazquez et al., 2013).

Detector modeling is typically not limited to particle transport in the atmosphere. Yakutsk folds the shower output through full simulations of the surface and underground scintillation detectors with the real detector geometry and thresholds; ALICE transports muons through 80 m w.e. of rock and the underground hall via GEANT3; IceTop/IceCube uses a Geant4-based detector simulation for both surface tanks and in-ice DOMs; HAWC passes CORSIKA output through a full GEANT4 model of the 300 water Cherenkov detectors (Glushkov et al., 2024, Collaboration, 2024, Verpoest, 2023, Arteaga-Velázquez, 21 Sep 2025).

The same model output can enter very different reconstruction pipelines. In the Yakutsk muon-correlation analysis, QGSJET-II-04 provides benchmark muon and surface densities for proton and iron primaries at a fixed reference energy. In IceCube, the model defines the simulated relation between primary energy and the multiplicity of muons with $E_\mu>500$ GeV. In HAWC, QGSJET-II-04 is treated as a fixed generator whose outputs populate the bidimensional response matrix $\Delta\sim 0.1\mbox{–}0.15$0 used in Gold-method unfolding (Glushkov et al., 2024, Verpoest, 2023, Arteaga-Velázquez, 21 Sep 2025).

This broad implementation pattern suggests that QGSJET-II-04 functions less as a standalone theoretical object than as infrastructure for inverse problems in cosmic-ray physics: composition inference, energy-scale setting, attenuation measurements, and detector-background estimation all depend on its forward model.

3. Model observables and analysis relations

QGSJET-II-04 enters air-shower analyses through a restricted set of observables that are strongly model dependent. The most common are the muon lateral distribution, the total or local muon content, the attenuation of muons with slant depth, and the surface-detector signal at a reference core distance.

In the Yakutsk muon-correlation method, one measures $\Delta\sim 0.1\mbox{–}0.15$1 and $\Delta\sim 0.1\mbox{–}0.15$2 and uses QGSJET-II-04 tables of $\Delta\sim 0.1\mbox{–}0.15$3 and $\Delta\sim 0.1\mbox{–}0.15$4 at a fixed reference energy $\Delta\sim 0.1\mbox{–}0.15$5 eV. The model prediction for an observed shower is then obtained through

$\Delta\sim 0.1\mbox{–}0.15$6

This avoids rerunning full Monte Carlo for each event while preserving the proton and iron benchmark lines used in composition classification (Glushkov et al., 2024).

The same Yakutsk analysis defines a mass-sensitive estimator

$\Delta\sim 0.1\mbox{–}0.15$7

and uses the ratio

$\Delta\sim 0.1\mbox{–}0.15$8

to separate proton-like, iron-like, muon-rich, and muon-deficient event groups (Glushkov et al., 2024).

For atmospheric development studies, the model is tested through the muon-attenuation length $\Delta\sim 0.1\mbox{–}0.15$9, defined by

$$\sigma_{p\text{-air}}(E)=\sigma_0\Bigl(\frac{E}{E_0}\Bigr)^\alpha, \qquad \sigma_0=265~\mathrm{mb},\quad E_0=0.2~\mathrm{TeV},\quad \alpha=0.079\pm0.010,$$0

or equivalently

$$\sigma_{p\text{-air}}(E)=\sigma_0\Bigl(\frac{E}{E_0}\Bigr)^\alpha, \qquad \sigma_0=265~\mathrm{mb},\quad E_0=0.2~\mathrm{TeV},\quad \alpha=0.079\pm0.010,$$1

This observable is used in both KASCADE-Grande and LHAASO to compare the measured zenith-angle dependence of muon content against QGSJET-II-04 predictions (Collaboration, 2024, Collaboration et al., 2018).

For surface-detector calibration, QGSJET-II-04 is used to derive energy–signal mappings. In the Yakutsk/Telescope Array comparison, the attenuation-corrected signal density at 800 m is converted into the SD energy estimator through

$$\sigma_{p\text{-air}}(E)=\sigma_0\Bigl(\frac{E}{E_0}\Bigr)^\alpha, \qquad \sigma_0=265~\mathrm{mb},\quad E_0=0.2~\mathrm{TeV},\quad \alpha=0.079\pm0.010,$$2

with $$\sigma_{p\text{-air}}(E)=\sigma_0\Bigl(\frac{E}{E_0}\Bigr)^\alpha, \qquad \sigma_0=265~\mathrm{mb},\quad E_0=0.2~\mathrm{TeV},\quad \alpha=0.079\pm0.010,$$3 and $$\sigma_{p\text{-air}}(E)=\sigma_0\Bigl(\frac{E}{E_0}\Bigr)^\alpha, \qquad \sigma_0=265~\mathrm{mb},\quad E_0=0.2~\mathrm{TeV},\quad \alpha=0.079\pm0.010,$$4 eV for the TA SD calibration. The same study reports the hybrid relation

$$\sigma_{p\text{-air}}(E)=\sigma_0\Bigl(\frac{E}{E_0}\Bigr)^\alpha, \qquad \sigma_0=265~\mathrm{mb},\quad E_0=0.2~\mathrm{TeV},\quad \alpha=0.079\pm0.010,$$5

for the Telescope Array fluorescence and surface energy scales (Glushkov et al., 2024).

4. Experimental performance across energy regimes

Published comparisons do not yield a single global verdict on QGSJET-II-04. Instead, the model performs well for some observables and regimes, while exhibiting systematic discrepancies for others.

Context	QGSJET-II-04 result	Comparison with data
Tien Shan, $$\sigma_{p\text{-air}}(E)=\sigma_0\Bigl(\frac{E}{E_0}\Bigr)^\alpha, \qquad \sigma_0=265~\mathrm{mb},\quad E_0=0.2~\mathrm{TeV},\quad \alpha=0.079\pm0.010,$$6–$$\sigma_{p\text{-air}}(E)=\sigma_0\Bigl(\frac{E}{E_0}\Bigr)^\alpha, \qquad \sigma_0=265~\mathrm{mb},\quad E_0=0.2~\mathrm{TeV},\quad \alpha=0.079\pm0.010,$$7 PeV	$$\sigma_{p\text{-air}}(E)=\sigma_0\Bigl(\frac{E}{E_0}\Bigr)^\alpha, \qquad \sigma_0=265~\mathrm{mb},\quad E_0=0.2~\mathrm{TeV},\quad \alpha=0.079\pm0.010,$$8 mb, $$\sigma_{p\text{-air}}(E)=\sigma_0\Bigl(\frac{E}{E_0}\Bigr)^\alpha, \qquad \sigma_0=265~\mathrm{mb},\quad E_0=0.2~\mathrm{TeV},\quad \alpha=0.079\pm0.010,$$9 rise per decade (Nesterova, 2015)	Reported as better matched than earlier QGSJET-II versions
EAS-MSU, $8\%\pm1\%$0–$8\%\pm1\%$1 eV, $8\%\pm1\%$2 GeV	$8\%\pm1\%$3 with 43% p + 57% Fe (Fomin et al., 2016)	No significant muon excess in inner region
ALICE multimuons, $8\%\pm1\%$4 eV	Iron sample gives MC/Data $8\%\pm1\%$5 for $8\%\pm1\%$6; proton gives $8\%\pm1\%$7 (Collaboration, 2024)	Only model reported to reproduce MMD reasonably well under heavy composition
Yakutsk UHE sample, $8\%\pm1\%$8 eV	$8\%\pm1\%$9, $$\sigma_{p\text{-air}}(1\,\mathrm{PeV})=350\pm15~\mathrm{mb}$$0 (Glushkov et al., 2024)	Enables separation of proton-like, iron-like, muon-rich, and muon-deficient groups
IceTop/IceCube, 2.5–100 PeV, $$\sigma_{p\text{-air}}(1\,\mathrm{PeV})=350\pm15~\mathrm{mb}$$1 GeV	Simulated and measured $$\sigma_{p\text{-air}}(1\,\mathrm{PeV})=350\pm15~\mathrm{mb}$$2 agree within 15–20% systematics (Verpoest, 2023)	Absolute TeV-muon multiplicity is consistent
KASCADE-Grande, $$\sigma_{p\text{-air}}(1\,\mathrm{PeV})=350\pm15~\mathrm{mb}$$3–$$\sigma_{p\text{-air}}(1\,\mathrm{PeV})=350\pm15~\mathrm{mb}$$4 eV	$$\sigma_{p\text{-air}}(1\,\mathrm{PeV})=350\pm15~\mathrm{mb}$$5 g/cm$$\sigma_{p\text{-air}}(1\,\mathrm{PeV})=350\pm15~\mathrm{mb}$$6 (Collaboration et al., 2018)	Below $$\sigma_{p\text{-air}}(1\,\mathrm{PeV})=350\pm15~\mathrm{mb}$$7 g/cm$$\sigma_{p\text{-air}}(1\,\mathrm{PeV})=350\pm15~\mathrm{mb}$$8
LHAASO, 0.3–30 PeV	$$\sigma_{p\text{-air}}(1\,\mathrm{PeV})=350\pm15~\mathrm{mb}$$9 predictions are systematically longer than measured (Collaboration, 2024)	Data favor EPOS-LHC over QGSJET-II-04
SUGAR, $pp$0–$pp$1 eV	LDF too flat; data are higher near core and lower at large $pp$2 (Kalmykov et al., 2022)	Shape mismatch rather than pure normalization mismatch

The model’s strongest reported successes in the data block are its compatibility with EAS-MSU inner-core muon densities once a surface-derived composition is imposed, and its simultaneous description of the ALICE multimuon multiplicity distribution and the high-multiplicity tail under a heavy-composition hypothesis (Fomin et al., 2016, Collaboration, 2024). The Yakutsk muon-correlation analysis further shows that QGSJET-II-04 can act as a stable benchmark for event-by-event classification at ultra-high energies, with the full 127-event sample yielding $pp$3, consistent with a predominantly light composition but with heavy, muon-rich, and muon-poor admixtures (Glushkov et al., 2024).

The principal weaknesses reported here concern the atmospheric evolution and spatial profile of the muon component. KASCADE-Grande and LHAASO both find that QGSJET-II-04 predicts steeper or deeper muon attenuation than observed, while SUGAR finds a lateral distribution that falls off too slowly with core distance (Collaboration et al., 2018, Collaboration, 2024, Kalmykov et al., 2022).

5. Muon-related tensions and the “muon puzzle”

The model’s record on muons is explicitly non-uniform. In some configurations QGSJET-II-04 removes or greatly reduces previously reported discrepancies, whereas in others it appears to shift rather than eliminate them.

At EAS-MSU, where the analysis is restricted to the inner shower region, $pp$4, and to muons above $pp$5 GeV, the measured and simulated $pp$6 distributions overlap within statistical errors and the fitted scaling factor is $pp$7; the paper concludes that no ad hoc muon rescaling is required in that regime (Fomin et al., 2016). By contrast, KASCADE-Grande reports that the measured attenuation of the muon content in the atmosphere is lower than predicted, with QGSJET-II-04 still underestimating $pp$8 despite being closer than older pre-LHC models (Collaboration et al., 2018). LHAASO extends that discrepancy across 0.3–30 PeV, with QGSJET-II-04 overpredicting the attenuation length by about $pp$9 g/cm$\rho^0$0 at 0.3 PeV and about $\rho^0$1 g/cm$\rho^0$2 at 30 PeV, corresponding to approximately 25% and 78% respectively (Collaboration, 2024).

The tension is not only longitudinal. SUGAR finds that, after normalizing to the same total muon number, QGSJET-II-04 places relatively too few muons near the core and too many at large distances. The reported residuals are up to $\rho^0$3–$\rho^0$4 for $\rho^0$5 m and roughly $\rho^0$6–$\rho^0$7 for $\rho^0$8 m, with a chance probability $\rho^0$9 for such a coherent slope difference (Kalmykov et al., 2022). IceCube identifies a different inconsistency: under QGSJET-II-04, TeV-muon multiplicities agree with simulation-based expectations, but the composition implied by those TeV muons is significantly heavier than the composition inferred from GeV-muon lateral densities measured by IceTop, with the two $\approx20\%$0 values differing by up to $\approx20\%$1 in $\approx20\%$2, or about $\approx20\%$3–$\approx20\%$4 in the 5–50 PeV range (Verpoest, 2023).

A further constraint comes from the vertical atmospheric muon spectrum. In that comparison, QGSJET-II-04 predicts a vertical muon intensity $\approx20\%$5 larger than data by a factor rising from about 1.4 at $\approx20\%$6 GeV to about 1.7 at $\approx20\%$7 GeV, and is therefore interpreted as producing too many very-forward secondary $\approx20\%$8 and $\approx20\%$9 mesons (Dedenko et al., 2015). This does not contradict the ALICE result, but it does indicate that improvement in one muon observable does not guarantee consistency across all muon thresholds, radial regions, and phase-space domains. A plausible implication is that the remaining discrepancies are differential: they concern energy partition, forward production, and transverse development rather than a single universal muon normalization.

6. Calibration, extensions, and related variants

Beyond shower-composition studies, QGSJET-II-04 serves as a calibration model for detector systems and as a reference generator in broader astroparticle pipelines. In the Yakutsk/Telescope Array calibration study, QGSJET-II-04 proton simulations are used to connect the particle density at 800 m to the SD energy scale, to compare Yakutsk and TA vertical-shower calibrations, and to interpret the hybrid result $E_{\rm lab}\gtrsim80$0 for TA (Glushkov et al., 2024). In the Yakutsk surface-detector composition analysis, QGSJET-II-04 is reported to give one of the best agreements with measured charged-particle lateral steepness, supporting an inferred decrease of $E_{\rm lab}\gtrsim80$1 from about 3 at $E_{\rm lab}\gtrsim80$2 eV to about 1.5 at $E_{\rm lab}\gtrsim80$3 eV (Glushkov et al., 2014).

In gamma-ray instrumentation, QGSJET-II-04 is close to QGSJET-II-03 in its practical impact on CTA sensitivity estimates. The CTA study attributes this to very similar $E_{\rm lab}\gtrsim80$4 production spectra: the QGSJET-II-04 spectrum is about 13% softer than QGSJET-II-03 at $E_{\rm lab}\gtrsim80$5, remains within about 10% over $E_{\rm lab}\gtrsim80$6–0.8, yields an average total residual background ratio $E_{\rm lab}\gtrsim80$7 across 1–30 TeV, and changes the 50 h point-source sensitivity by only $E_{\rm lab}\gtrsim80$8 (Ohishi et al., 2021).

A closely related derivative is QGSJET-II-04m, used in the AAfrag interpolation package for inclusive secondary production in $E_{\rm lab}\gtrsim80$9, $E_\mu>500$0, $E_\mu>500$1, and $E_\mu>500$2 collisions. QGSJET-II-04m is described as a successor with improved low-energy secondary production, retuned hadronization, and validation against LHCf, LHCb, and NA61 data. AAfrag exposes tabulated differential yields for photons, neutrinos, leptons, and antinucleons through interpolation routines rather than through full air-shower simulations (Kachelriess et al., 2019). This suggests an important bifurcation in the model family: QGSJET-II-04 remains the standard EAS interaction engine in many detector simulations, while QGSJET-II-04m extends the framework toward precision secondary-production calculations relevant to Galactic cosmic-ray transport.

Overall, QGSJET-II-04 occupies a technically central but empirically conditional position in high-energy cosmic-ray phenomenology. It is well enough constrained to underpin detector calibration, response-matrix construction, and event-by-event mass-sensitive benchmarks, yet current measurements of attenuation lengths, muon spectra, and lateral distributions show that its description of muon production and transport remains incomplete in several energy and phase-space regimes (Glushkov et al., 2024, Arteaga-Velázquez, 21 Sep 2025, Collaboration, 2024).