SIBYLL 2.3d: High-Energy Air Shower Model

• SIBYLL 2.3d is a post-LHC hadronic interaction model that simulates cosmic-ray air showers by updating particle production, cross sections, and charm inclusion based on collider measurements.
• The model is integrated into simulation frameworks like CORSIKA 8 and CONEX, providing detailed predictions of key observables such as shower maximum (Xmax) and muon content (Nμ) for composition studies.
• Phenomenological modifications in SIBYLL 2.3d, such as enhanced baryon-pair production and leading-particle treatment, aim to reconcile simulated muon numbers with experimental data.

Searching arXiv for the primary SIBYLL 2.3d paper and closely related validation/application papers. SIBYLL 2.3d is a hadronic interaction model for cosmic-ray air-shower simulations. It is the modern SIBYLL version described in “Hadronic interaction model Sibyll 2.3d and extensive air showers” (Riehn et al., 2019), and it is used as a high-energy interaction generator in contemporary simulation frameworks and data analyses spanning extensive air showers, muon studies, composition inference, and detector-specific response modeling. In the literature represented here, SIBYLL 2.3d is characterized as a post-LHC model whose updates affect particle production, cross sections, leading-particle treatment, baryon-pair production, and charm production, with consequences for standard shower observables such as $X_{\rm max}$ and $N_\mu$ (Riehn et al., 2019).

1. Definition and position within hadronic air-shower modeling

SIBYLL 2.3d is a high-energy hadronic interaction model used inside air-shower simulation chains. In CORSIKA-based workflows, it appears as the high-energy model above a low-energy transition handled by FLUKA in several analyses, with explicit thresholds of $E<200$ GeV in KASCADE-Grande analyses (Kang et al., 2023), and a default high/low transition “fixed at ∼ 80 GeV as default, but can be modified” in CORSIKA 8 (Gaudu, 2024). Within CORSIKA 8, it is one of the standard high-energy models alongside QGSJet-II.04, EPOS-LHC, and Pythia 8.3, while decays are supported through SIBYLL 2.3d and Pythia 8 (Gaudu, 2024).

The model’s role is to provide hadronic cross sections, multiplicities, and secondary-particle spectra for hadron–air interactions that drive the hadronic cascade in extensive air showers. In that sense, it occupies the same functional layer as other post-LHC hadronic models, but with its own particle-production systematics and corresponding implications for muon production, shower maximum, and composition-sensitive observables (Kang et al., 29 Aug 2025). In CORSIKA 8, this role is embedded in a modular C++ process architecture, where SIBYLL 2.3d can be selected as a plug-in high-energy interaction model and also as a decay provider (Huege et al., 2023).

The dedicated SIBYLL 2.3d paper states that the new version preserves the core ideas of the model while changing the handling of baryon pairs and leading particles, adding charm production, and updating the extrapolation to ultra-high energy using high-precision measurements of the total and inelastic cross sections with the forward detectors at the LHC; minimum-bias measurements of particle spectra and multiplicities support the tuning of fragmentation parameters (Riehn et al., 2019). This situates SIBYLL 2.3d as a post-LHC event generator intended specifically for extensive air-shower applications.

2. Model content and high-energy phenomenology

The most explicit summary of SIBYLL 2.3d’s internal evolution comes from its dedicated description: the model introduces a new treatment of baryon-pair production and leading particles, includes production of charmed hadrons, and uses LHC forward-detector constraints on total and inelastic cross sections, together with minimum-bias particle spectra and multiplicities, to tune fragmentation parameters (Riehn et al., 2019). In related summaries, SIBYLL is described as being based on the dual parton and minijet models (Kang et al., 2021), and as a model in which soft and semihard production, diffraction, and multiparticle production are central to the event generator (Collaboration, 2024).

Several downstream analyses identify specific phenomenological consequences of the 2.3d update. In KASCADE-Grande-focused work, SIBYLL 2.3d is described as having been developed by improving the behavior of the $\pi^{\pm}/\pi^{0}$ ratio in different hadronization mechanisms, with the aim of increasing the hadronic, muon-producing channel relative to the electromagnetic channel (Kang et al., 2022). The same source states that, relative to SIBYLL 2.1, the muon number increases by more than 20%, and that SIBYLL 2.3d predicts muon numbers about 5% higher than EPOS-LHC and QGSJET-II-04 (Kang et al., 2022). This is consistent with the broader view that SIBYLL 2.3d is comparatively muon-rich among post-LHC interaction models.

The model is also treated as intermediate among post-LHC descriptions of shower development in some contexts. In the KASCADE-Grande post-LHC spectrum analysis, EPOS-LHC-R is described as having a deeper $X_{\max}$ and higher muon count than EPOS-LHC, even “deeper than the predictions from the SIBYLL 2.3d model,” which places SIBYLL 2.3d between different post-LHC options in terms of shower development and muon content (Kang et al., 29 Aug 2025). A plausible implication is that SIBYLL 2.3d frequently acts as a reference model not only because it is modern and collider-constrained, but also because it occupies a useful midpoint in several observable spaces.

3. Implementation in simulation frameworks

SIBYLL 2.3d is deeply integrated into contemporary shower-simulation software. In “The particle-shower simulation code CORSIKA 8” (Huege et al., 2023) and its later status update (Gaudu, 2024), it is explicitly listed as one of the supported high-energy hadronic interaction models in a modular C++ framework that has reached a “physics-complete” state. In this environment, SIBYLL 2.3d is part of a process list that can be assembled flexibly, and it can handle both high-energy hadronic interactions and decay processes (Huege et al., 2023).

The default high/low-energy split in CORSIKA 8 is “fixed at ∼ 80 GeV as default, but can be modified” (Gaudu, 2024). This means that if SIBYLL 2.3d is selected as the high-energy model, it governs interactions above that threshold, while FLUKA handles the low-energy regime. The same CORSIKA 8 status paper also notes that “Decay processes are supported through Sibyll 2.3d and Pythia 8,” so the model’s integration extends beyond collision generation to unstable-hadron handling (Gaudu, 2024).

The earlier CORSIKA 8 overview further notes that Sibyll 2.3d can also contribute to photohadronic treatment in electromagnetic cascades, in combination with SOPHIA, and that it can be used in multiple-media and cross-media simulation contexts (Huege et al., 2023). This does not mean the CORSIKA papers unpack the internal physics of SIBYLL 2.3d in detail; rather, they define its position in the software stack and demonstrate that it is a standard component of modern air-shower workflows.

Outside CORSIKA 8, SIBYLL 2.3d is also used in CONEX-based studies. In the deep-learning study on longitudinal shower profiles, the network is trained entirely on EPOS-LHC showers and evaluated on a 1,000,000-shower CONEX+SIBYLL-2.3d library, specifically to quantify hadronic-model dependence (Wang et al., 19 Aug 2025). This usage underscores that SIBYLL 2.3d functions not only as a production model for baseline simulations but also as a robust out-of-domain test model for model-systematics studies.

4. Validation against shower observables

A recurring theme across the literature is that SIBYLL 2.3d is validated through standard air-shower observables rather than through internal parameter exposition. In CORSIKA 8 validation, SIBYLL 2.3d participates in comparisons of proton-induced showers where the agreement between CORSIKA 8 and CORSIKA 7 for electrons/positrons, photons, muons, and hadrons is within 10% (Huege et al., 2023). The later CORSIKA 8 status paper further specifies that, in lateral distributions at ground for proton-induced showers, the muon content agrees at the 5–10% level near the shower axis, while discrepancies at larger distances remain to be explored (Gaudu, 2024). These comparisons include SIBYLL 2.3d together with EPOS-LHC and QGSJet-II.04.

At the same time, the CORSIKA 8 overview notes that with Sibyll a larger spread between muon predictions across various codes has been observed previously (Huege et al., 2023). This is not presented as a defect unique to SIBYLL 2.3d, but as a reminder that muon-related observables remain an especially sensitive diagnostic of hadronic-model implementation and shower-code coupling.

Experimental comparisons reveal a mixed but informative picture. In ALICE underground multimuon data, SIBYLL 2.3d underpredicts the muon multiplicity distribution in a large multiplicity interval by more than 20% under heavy-composition assumptions, while QGSJET-II-04 reproduces the distribution more reasonably and EPOS-LHC underpredicts it by more than 30% (Collaboration, 2024). However, for high muon multiplicity events with $N_\mu>100$, the rate obtained with SIBYLL 2.3d is compatible with the data when the composition is assumed to be dominated by heavy elements (Collaboration, 2024). This dual outcome is characteristic: SIBYLL 2.3d can succeed for some integral or extreme-event observables while remaining under tension for differential muon distributions.

In LHAASO’s measurement of the attenuation length of the muon content from 0.3 to 30 PeV, SIBYLL 2.3d predicts attenuation lengths significantly longer than the data, while EPOS-LHC is relatively closer to the measurement (Collaboration, 2024). The paper explicitly states that LHAASO data favor EPOS-LHC over QGSJET-II-04 and SIBYLL 2.3d (Collaboration, 2024). This indicates that SIBYLL 2.3d’s muon-development systematics are not uniformly preferred across experimental environments and energy ranges.

5. Composition and spectrum inference

SIBYLL 2.3d is widely used in composition and spectrum reconstruction because its predictions feed directly into energy calibration, mass-group separation, response matrices, and unfolding. In KASCADE-Grande analyses, the high-energy interaction model is combined with FLUKA below 200 GeV, and energy is reconstructed from charged-particle shower size via

$$\log_{10}\left(\frac{E}{\mathrm{GeV}}\right)=a\cdot\log_{10}(N_{ch})+b.$$

For SIBYLL 2.3d, model-specific coefficients are derived separately for light and heavy primaries in different analyses (Kang et al., 2023, Kang et al., 2022).

In the 2025 KASCADE-Grande post-LHC spectrum analysis, SIBYLL 2.3d yields a heavy-component knee at

$\log_{10}(E_k/\mathrm{GeV})=7.78\pm0.02,$

with spectral indices

$\gamma_1=2.73\pm0.01,\qquad \gamma_2=3.24\pm0.03,$

and $\Delta\gamma=0.51$, with $N_\mu$0 (Kang et al., 29 Aug 2025). The same study describes a light-component hardening above about $N_\mu$1, an ankle-like feature seen across post-LHC models including SIBYLL 2.3d (Kang et al., 29 Aug 2025). The all-particle spectrum depends only weakly on model choice in that analysis, and SIBYLL 2.3d lies very close to QGSJet-II-04 and EPOS-LHC in flux and shape (Kang et al., 29 Aug 2025).

Earlier KASCADE-Grande work emphasizes that SIBYLL 2.3d gives the lightest inferred composition among the considered models because it predicts the lowest flux of heavy primaries, a consequence of its comparatively high muon content (Kang et al., 2022, Kang et al., 2021). Yet the same works stress that the main spectral structures—heavy knee, light hardening, concave behavior around $N_\mu$2 eV—remain robust (Kang et al., 2022, Kang et al., 2021).

At ultra-high energy, SIBYLL 2.3d also appears centrally in “heavy-metal” scenarios. In one such analysis, SIBYLL 2.3d is one of the two reference models used to interpret Pierre Auger Observatory $N_\mu$3 data under a pure-iron-above-$N_\mu$4 eV assumption (Vícha et al., 12 Feb 2025). The model’s $N_\mu$5 scale is shifted by

$N_\mu$6

to achieve consistency with the data (Vícha et al., 12 Feb 2025). A later data-driven formulation of the heavy-metal scenario reuses this global $N_\mu$7 recalibration of SIBYLL 2.3d and interprets it as leading to a heavier UHECR composition than conventionally assumed (Vícha et al., 31 Jul 2025). These studies do not prove that SIBYLL 2.3d is uniquely correct; rather, they show that modest global recalibration of its $N_\mu$8 scale can support a self-consistent heavy-composition interpretation.

6. Muons, model tensions, and phenomenological modifications

The most persistent controversy surrounding SIBYLL 2.3d concerns muons. On the one hand, several analyses describe it as relatively muon-rich compared with older SIBYLL versions and even compared with some contemporary models (Kang et al., 2022, Kang et al., 2021). On the other hand, Auger hybrid analyses still require larger hadronic signals at ground than Sibyll 2.3d predicts by default (Collaboration et al., 2024). In that study, the best fit requires the hadronic component of $N_\mu$9 to be rescaled by about 15–18% and $E<200$0 to be shifted deeper by about 10 g/cm$E<200$1 for SIBYLL 2.3d (Collaboration et al., 2024). This suggests that even a comparatively muon-rich post-LHC model remains insufficient for a fully satisfactory joint description of longitudinal and ground observables.

That tension motivates explicit phenomenological deformations of SIBYLL 2.3d. In “Sibyll$E<200$2: ad-hoc modifications for an improved description of muon data in extensive air showers” (Riehn et al., 2023) and the later “Sibyll$E<200$3” paper (Riehn et al., 2024), the baseline model is modified after event generation to enhance $E<200$4, baryon–antibaryon pair, or kaon production. These variants remain within accelerator bounds comparable to Sibyll 2.3d and can increase the muon count in extensive air showers by up to 35%, while minimally affecting $E<200$5 (Riehn et al., 2024). The studies argue that enhanced $E<200$6 production and mixed $E<200$7+baryon modifications can bring the model closer to Auger muon data (Riehn et al., 2023, Riehn et al., 2024).

A crucial implication, stated explicitly in the modified-model work, is that increasing muon production reduces the difference in predicted total muon numbers for proton and iron showers (Riehn et al., 2023). This means that using such variants may improve agreement with observed muon content but can weaken mass discrimination in total $E<200$8. This is not a contradiction; it is a reminder that fixing one hadronic observable can compress composition sensitivity in another.

A related but distinct line of work modifies effective hadronic characteristics of Sibyll 2.3d—cross section, elasticity, and multiplicity—to see whether reasonable macro-parameter changes can reconcile the muon excess and $E<200$9 behavior (Blazek et al., 2023). That study finds strong correlations among the resulting shifts in observables and concludes that simple cross-section, elasticity, and multiplicity rescalings do not naturally reproduce both more muons and deeper $\pi^{\pm}/\pi^{0}$0 simultaneously. A plausible implication is that the missing physics is more microscopic than a global rescaling of a few bulk parameters.

7. Broader applications and model dependence in modern analyses

SIBYLL 2.3d is used well beyond classical spectrum reconstruction. In IceCube’s improved gamma–hadron separation for PeV photon searches, it is the central hadronic interaction model in CORSIKA simulations used to optimize discriminants and estimate photon survival fractions (Schröder et al., 8 Jul 2025). The paper states that gamma–hadron separation is better than $\pi^{\pm}/\pi^{0}$1, in the sense that measured air-shower events are suppressed more than 1000 times more strongly than photon-induced showers of the same energy simulated with Sibyll 2.3d (Schröder et al., 8 Jul 2025). The hadronic background is anchored in data rather than simulation, so SIBYLL 2.3d primarily controls the photon-efficiency side of the analysis.

In machine-learning studies of composition, SIBYLL 2.3d is often used as the out-of-domain test model. In the longitudinal-profile neural-network study, an EPOS-LHC-trained network applied to SIBYLL-2.3d showers still achieves a proton–iron merit factor of about 2.0, versus about 2.2 on EPOS, but predicts heavier nuclei too lightly on the absolute $\pi^{\pm}/\pi^{0}$2 scale (Wang et al., 19 Aug 2025). This shows that full-profile composition-sensitive features are relatively robust across hadronic models, while absolute mass calibration remains model-dependent.

In knee-region composition studies based on observables designed to suppress model dependence, SIBYLL 2.3d behaves similarly to EPOS-LHC and QGSjet-II-04. The $\pi^{\pm}/\pi^{0}$3 azimuthal-fluctuation parameter is reported to have minimal dependence on the specific hadronic interaction model considered, including SIBYLL 2.3d (Arsene, 2023). Likewise, a PCA-based KASCADE analysis using $\pi^{\pm}/\pi^{0}$4, $\pi^{\pm}/\pi^{0}$5, $\pi^{\pm}/\pi^{0}$6, and lateral age finds that, although the input observables remain model-dependent, the resulting mass-composition reconstruction significantly reduces dependence on the hadronic interaction model used in simulation, with SIBYLL 2.3d yielding results consistent with EPOS-LHC and QGSjet-II-04 within uncertainties (Arsene, 2 Jul 2025).

Taken together, these applications suggest that SIBYLL 2.3d is not merely a default shower model but a reference point for systematic studies: a model modern enough to be operationally standard, yet distinct enough in its muon and shower-development predictions to make it valuable for cross-model robustness tests.