Cosmic Ray Simulation with PYTHIA

We present recent developments in PYTHIA for the modelling of hadronic cascades in a medium. Several improvements have been made in the Angantyr model for collisions with nuclei, especially in the limit of low collision energies, allowing it to be used throughout the hadronic cascades. Also the simplified nuclear model in the PythiaCascade} module has been updated. We find that the two models give consistent results for cosmic ray air showers initiated by both high energy protons and nuclei.

💡 Research Summary

This paper reports on recent developments in the PYTHIA 8.3 event generator that enable it to serve as a full‑featured hadronic interaction model for cosmic‑ray air‑shower simulations. The authors focus on two complementary modules: Angantyr, which provides a Glauber‑based description of nucleus–nucleus collisions, and the simpler PythiaCascade (hereafter “Cascade”) add‑on, which implements a fast, parametrised nuclear geometry. Historically, Angantyr has excelled at high‑energy (TeV) proton–nucleus and nucleus–nucleus collisions, but it lacked a reliable treatment of low‑energy (∼10 GeV) interactions, a clear separation of elastic, inelastic and diffractive components, and the ability to switch rapidly between different projectile–target combinations as required by a cascade evolving through the atmosphere. Cascade, on the other hand, could handle such rapid switching but used a very crude wounded‑nucleon count and did not model the full nuclear geometry.

The paper introduces several key upgrades. First, the total, elastic, and diffractive cross sections for a broad set of hadrons (π, K, η, D, B, Λ, Σ, Ξ, Ω and their charm/bottom counterparts) are now implemented using the Donnachie–Landshoff (DL) parametrisation together with the SaS extension. The pomeron term X is scaled according to the Additive Quark Model (proportional to the number of valence quarks, with mass‑dependent suppression for heavy flavours), while the reggeon term Y is estimated from light‑quark exchange arguments. This yields a consistent set of cross sections down to the threshold region, where non‑perturbative handling (gluon exchange with colour‑octet conversion) is applied up to √s≈10 GeV, after which the standard perturbative multi‑parton interaction (MPI) machinery takes over.

Second, a new family of parton distribution functions (PDFs), labelled SU21, has been generated for all 21 projectile species. Starting from simple valence‑quark forms at Q0²≈0.26 GeV², the PDFs are evolved with QCDNUM to the relevant scales and stored in grid files. This allows MPI calculations for pions, kaons, heavy mesons and baryons, which previously relied on crude approximations. The regularisation scale p⊥0 is kept common to all hadrons, but the framework permits future species‑dependent tuning.

Third, the authors pre‑tabulate the integrated MPI cross sections and diffractive mass spectra for each projectile over a logarithmic energy grid (starting at 10 GeV). During a cascade simulation, the required values are obtained by fast interpolation, reducing the per‑collision overhead from seconds to milliseconds. This speed‑up is essential because a realistic air‑shower involves millions of secondary interactions with constantly changing projectile types and energies.

Fourth, the Angantyr implementation has been extended to include fluctuations of the nucleon–nucleon cross section via a Gamma‑distributed interaction radius, following the Good–Walker formalism. The elastic amplitude is taken as purely imaginary, with a radius‑dependent profile T(b,σ). Elastic, single‑diffractive, double‑diffractive and central‑diffractive processes are now distinguished, each with its own parametrised mass spectrum (approximately d m²/m²) and low‑mass enhancement. This provides a more realistic treatment of low‑energy diffractive events that can dominate the early stages of a cascade.

Fifth, the forward fragmentation region—critical for air‑shower development because it determines the inelasticity of each interaction—has been refined. The authors reduce the “popcorn” mechanism for diquark break‑up and adjust the Lund symmetric fragmentation parameters a and b for diquark ends, making forward baryons harder (i.e., carrying a larger fraction of the projectile momentum). A bug in the handling of a‑dependent diquark fragmentation, discovered in PYTHIA 8.315, has been fixed, which appears to alleviate the need for further b‑tuning.

The paper then presents validation studies. Simulations of proton‑initiated and iron‑initiated air showers are performed with both the upgraded Angantyr and the Cascade model. Key observables—depth of shower maximum (Xmax), muon content, and the electromagnetic‑to‑hadronic energy ratio—are compared. The two models yield consistent results, confirming that the simplified Cascade geometry reproduces the full Angantyr predictions while being computationally cheaper. Moreover, the forward spectra generated by the new fragmentation settings show a reduced softness compared with the default PYTHIA tunes, bringing the predictions closer to measurements from the Pierre Auger Observatory.

Finally, the authors discuss integration prospects. The upgraded PYTHIA code can now serve as the hadronic‑interaction plugin for the next‑generation air‑shower framework Corsika 8, which already uses PYTHIA 8 for particle decays. Because both Corsika 8 and PYTHIA are written in C++, the interface is straightforward. The same code could also be linked to Geant4, enabling a unified description of hadronic interactions in collider detectors and atmospheric cascades. Future work will address neutrino‑ and photon‑initiated cascades, further tuning of low‑energy cross sections against fixed‑target data, and systematic studies of the impact of the new PDFs on muon production.

In summary, this work expands PYTHIA’s applicability from collider physics to cosmic‑ray physics by providing a comprehensive, low‑energy‑aware nuclear interaction model (Angantyr) together with a fast, parametrised alternative (Cascade). The developments ensure consistent treatment of elastic, inelastic, and diffractive processes across a wide range of projectiles and energies, improve forward fragmentation, and deliver the computational performance required for large‑scale air‑shower simulations.