Hadronic interaction model SIBYLL 2.3d and extensive air showers

We present a new version of the hadron interaction event generator SIBYLL . While the core ideas of the model have been preserved, the new version handles the production of baryon pairs and leading particles in a new way. In addition, production of charmed hadrons is included. Updates to the model are informed by high-precision measurements of the total and inelastic cross sections with the forward detectors at the LHC that constrain the extrapolation to ultrahigh energy. Minimum-bias measurements of particle spectra and multiplicities support the tuning of fragmentation parameters. This paper demonstrates the impact of these changes on air-shower observables such as X max and N μ , drawing comparisons with other contemporary cosmic-ray interaction models.


I. INTRODUCTION
Studying cosmic rays at energies above 100 TeV imposes a challenge since the intensity is too low for direct measurements with high-altitude balloons or spacecraft. Instead the properties of the primary cosmic-ray nucleus must be inferred indirectly from the properties of extensive air showers (EAS) that can be observed with large, ground-based detectors. At energies in excess of several tens or hundreds of PeV (so-called ultrahigh energy cosmic-rays (UHECR) ) the event rate per unit area and solid angle quickly drops, requiring ever larger and more sparsely instrumented detectors. Therefore, the interpretation of these air-shower data has necessarily to rely on detailed Monte Carlo simulations of the shower development and the experimental observables. The main challenge in these simulations is the modeling of nuclear and hadronic interactions that can occur at all possible energies ranging from the MeV up to ultrahigh energies E ∼ 10 21 eV. While interactions of hadrons with protons and nuclei are well studied up to several hundreds of GeV (in target rest frame) at fixed target detectors, at the highest energies it is necessary to rely on model extrapolations from collider experiments that measure primarily the central region. This leads to the subclass of event generators in high-energy physics called cosmic-ray interaction models. SIBYLL is one of the first microscopic event generators for EAS [1] and it is based in its core on the dual parton model (DPM) [2] and the minijet model [3][4][5][6]. Particle formation (or hadronization) is adopted from the Lund algorithms [7,8] and shares in this sense many ideas about the interactions of color strings with the popular PYTHIA event generators [9]. A summary of the principles and ideas behind SIBYLL and a review of its long history can be found in Ref. [10].
From the beginning SIBYLL aimed to describe a broad range of ppðpÞ measurements at the intersecting storage rings (ISR), the SpðpÞS at CERN and the Tevatron at Fermilab, providing the highest interaction energies available at that time; for example, the growth of the average transverse momentum with center-of-mass (c.m.) energy is adjusted according to the results of the CDF experiment at the Tevatron, UA1 at the SppS and the ISR at CERN [11][12][13]. The hard interaction cross section is calculated in the minijet model. The Glauber scattering theory [14] is applied in hadron-nucleus collisions and extended with a semisuperposition approach [15] to nucleus-nucleus collisions.
Since the previous version 2.1 [16] soft interactions and diffraction dissociation are implemented in a more sophisticated way by including multiple soft interactions and a two-channel eikonal model for diffraction, respectively. The current extension of the model is motivated by recent developments in cosmic-ray (CR) and astroparticle physics and new measurements at accelerators. At the high-energy frontier, the LHC provides for the first time constraints on extrapolation of the model to energies corresponding to cosmic rays beyond the knee. In addition, dedicated forward physics experiments (for example LHCf and CASTOR) and recent fixed target experiments (NA61) studied a larger part of the phase space that is particularly important for EAS.
There are several challenges for the present cosmic-ray interaction models. One example arises in the interpretation of EAS data in terms of CR mass composition where simulations predict a lower muon content than required to interpret the observations [17][18][19]. This challenge is specifically addressed by careful evaluation of ρ 0 and p=p production, both of which increase muon content in EAS.
Another example is the need to include production of charmed hadrons in event generators for EAS. The observation of high-energy astrophysical neutrinos above 100 TeV by IceCube [20,21] extends to the energy range where prompt muons and neutrinos from decays of charmed hadrons become larger than the conventional (light meson) channels. Eventually prompt muons and neutrinos become the main leptonic backgrounds for the astrophysical neutrino flux. Production of charm was first introduced as a modification of SIBYLL 2.1 [22]. Its implementation in SIBYLL 2.3d is based on comparison with recent accelerator data on production of charmed hadrons and fully supports the production of charm [23][24][25]. The model of the production of charm and the application of SIBYLL 2.3d to the calculation of inclusive lepton fluxes is the subject of a separate paper [25].
The objective of this paper is twofold. The first is a description of the post-LHC version SIBYLL 2.3d. 1 The changes to the microscopic interaction model with respect to the predecessor are detailed in Sec. II. The second objective is the evaluation of the impact on EAS observables. Section III contains the benchmark calculations and comparisons against other contemporary post-LHC models [28,29] including the previous SIBYLL 2.1. We conclude with a discussion in Sec. IV.

A. Basic model
The aim of the event generator SIBYLL is to account for the main features of strong interactions and hadronic particle production as needed for understanding air-shower cascades and inclusive secondary particle fluxes due to the interaction of cosmic rays in Earth's atmosphere. Therefore, the focus is on the description of particle production at small angles and on the flow of energy in the projectile direction. Rare processes, such as the production of particles or jets at large p T or electroweak processes, are either included approximately or neglected.
The model supports interactions between hadrons (mostly nucleons, pions or kaons) and light nuclei (h-A), since the targets in EAS mainly are nitrogen and oxygen. The CR flux at the top of the atmosphere contains elements up to iron, requiring a model for interactions of nuclei (A-A). Nuclear binding energies have negligible impact for high-energy interactions, allowing for the approximate construction of interactions of cosmic-ray nuclei from individual hadron-nucleon (h-N) collisions. On the target side, nucleons are combined to light nuclei on amplitude level using the Glauber model [14,30] together with the semisuperposition [15] approach. This means that the interaction of an iron nucleus (A ¼ 56), for example, with a nitrogen nucleus in air is treated as 56 separate nucleonnitrogen interactions. With the exception of inelastic screening (Sec. II E), the model extensions discussed in the following are introduced at the level of hadron-nucleon interactions.

Parton level
The total scattering amplitude that determines the interaction cross sections is defined in impact parameter space by using the eikonal approximation (see Refs. [6,31,32] and, for a pedagogical introduction, also Ref. [33]) 1 Preliminary versions of this model were released as SIBYLL 2.3 [26] and SIBYLL 2.3c [25,27]. Explanations of the changes between versions can be found in Appendix A.
where i is the unit imaginary number, ⃗ b is the impact parameter of the collision and s is the Mandelstam variable, which for the interaction between hadrons k and l is defined as s ¼ ðp k þ p l Þ 2 . The eikonal function χ is given by the sum of two terms representing soft and hard interactions χðs; ⃗ bÞ ¼ χ soft ðs; ⃗ bÞ þ χ hard ðs; ⃗ bÞ and then unitarized as in Eq. (1) (jaj ≤ 0.5). The soft and hard eikonal functions take the form with R A int ðs; ⃗ bÞd 2 ⃗ b ¼ 1 and int ¼ soft; hard. Within the parton model, there is a straightforward interpretation of Eq. (2) for hard interactions of asymptotically free partons. Then σ hard is the inclusive hard scattering cross section of partons in the interaction of hadron k with hadron l. The spatial distribution of partons available for hard interaction is encoded in the overlap function A hard ðs; ⃗ bÞ. This overlap function between hadrons k and l is given by the individual transverse profile functions of partons in the scattering hadrons, A k=l ðs; ⃗ b l Þ, and the transverse profile of the individual parton-parton interaction, A par ðs; ⃗ b par Þ: where ⃗ b k=l are the positions of the interacting partons in the hadrons k and l and ⃗ b par is the impact parameter between the partons; see Ref. [16]. For pointlike parton-parton interactions, A par would be a Dirac δ function.
A geometrical (gluon) saturation condition [16,34,35] is approximated by an energy-dependent transverse momentum p min ⊥ ðsÞ cutoff that separates soft and hard parton interactions: Values of the parameters can be found in Appendix C. Hard interactions are calculated in leading-order quantum chromodynamics (QCD) at the minimal scale p min ⊥ ðsÞ with a K factor to account for higher-order corrections. The hard interaction is assumed to be pointlike and the partons are spatially distributed inside the hadron according to the electric form factor of the proton [36]. The distribution of partons in momentum space is given by the parton distribution functions (PDFs) parameterized by Glück, Reya and Vogt [37,38].
The parametrization of the soft cross section is inspired by the Donnachie-Landshoff model [39]. The soft cross section has two components, one declining and one increasing with energy, corresponding to Reggeon and Pomeron exchange. In contrast to the hard parton interactions, the soft interactions are thought of as spatially extended; i.e., A soft ðs; ⃗ b z Þ in Eq. (4) is given by a Gaussian profile instead of Dirac's delta function. The width of the profile is energy dependent B s ðsÞ ¼ B 0 þ α 0 i ð0Þ ln ðs=s 0 Þ, with α 0 ð0Þ being a parameter known from Regge phenomenology; see, for example, [40,41]. To obtain an analytic solution for the overlap integral [Eq. (4)], the distribution of soft partons [A x;y ðs; ⃗ b x;y Þ] is defined as Gaussian, i.e., for a pp collision The effective width parameter 2B p þ B 0 is determined from a fit to cross section data and the slope of the energy dependence α 0 i ð0Þ is given by the slope of the Pomeron (Reggeon) trajectory known from soft interactions [42]. The interaction cross sections are calculated by integration of the above amplitude in impact parameter space, e.g., for the inelastic cross section The obtained values are given in Appendix C. A twochannel Good-Walker formalism is used for low-mass diffractive interactions, where the two channels correspond to the hadron's ground state and a generic excited state [43]. For simplicity, high-mass diffraction is assumed to account for 10% of the nondiffractive interactions and contributes with only a single cut. A more in depth discussion of the basic principles of the model can be found in Ref. [16]. The partial cross sections for multiple Pomeron scattering are calculated from the elastic amplitude using unitarity cuts (Abramovsky-Gribov-Kancheli cutting rules) [44]. The multiple cuts (or parton interactions) are assumed to be uncorrelated and Poisson-distributed at tree level, but at later steps of the event generation correlations can arise from e.g., energy and momentum conservation. The cross sections for multiple cuts are calculated (neglecting diffractive channels) from where N soft;hard is the number of soft or hard parton scatterings in the interaction. n int ðs; ⃗ bÞ ¼ 2χ int ðs; ⃗ bÞ is the average number of soft or hard interactions.
For runtime optimization the momenta of the partons in an event are sampled from approximate parameterizations instead of the full amplitude. The hard component (σ QCD ) is calculated at leading order assuming collinear factorization, in which the full PDFs that resolve individual quark flavors and gluons are replaced by an effective PDF for all partons of the form fðxÞ ¼ gðxÞ þ 4 9 ½qðxÞ þqðxÞ, where qðxÞ represents the combined distribution of all quark flavors [45]. Neglecting initial transverse momentum, the transverse momentum of the partons is determined by the scattering process given byt −2 , wheret is the four momentum transfer after Mandelstam.
For the soft interaction, which are assumed to include the valence quarks, the momentum fractions are taken from the distribution In case of the valence quarks, d, which leads to the suppression of large momentum fractions, is set to 3 (2) for baryons (mesons). The pole at small momentum fractions is controlled by the choice of an effective quark mass m 2 q ¼ 0.3 GeV 2 . For soft sea quarks and gluons, d ¼ 1.5 and m 2 q ¼ 0.01 GeV 2 . The conservation of energy is enforced by assigning one (the last) parton the remaining fraction. Since these distributions favor small momentum fractions, the remainder usually constitutes the largest fraction and thus emerges as leading particle. For baryons this fraction is always assigned to pairs of valence quarks, the so-called diquarks. For mesons one of the valence quarks is randomly selected as leading.
The excitation mass M D for diffractive interactions is sampled from a M −2 D distribution without distinguishing between the contributions from low-and high-mass diffraction. The minimal mass of the diffractively excited system is chosen such that the difference between the mass of the excited system and the original projectile hadron is larger than 1.5, 0.2 and 0.6 GeV for protons, pions and kaons, respectively. The upper limit for the diffractive mass universally is set to M 2 D =s ¼ 0.2. The transverse momentum in the diffractive interaction is assumed to be exponential in p 2 T with a slope [46,47].

Hadron level
The hadronization model in SIBYLL is based on the Lund string fragmentation model [8,48]. Each (nondiffractive) interaction involves the exchange of color between the hadrons. For the valence quarks a single soft gluon (two colors) is exchanged forming two color fields (strings) between the two quark-diquark pairs for baryons and quark-antiquark pair for mesons, respectively (Fig. 1). Since gluon scattering is the dominant process at high energy, all the additional hard or soft interactions are modeled as gluon-gluon scattering. Furthermore, the color flow of the gluon scattering is approximated by a closed color loop between two gluons resulting in two strings (see Fig. 2). In general, a single hadron-hadron interaction will be a complex combination of such two string configurations, where the probability density for the multiple cut (or string) topology is determined by σ N soft ;N hard [Eq. (8)].
The fraction of the string energy z assigned to the quarks in each step in the fragmentation is taken from the symmetric Lund function [49] where a ¼ 0.5 and κ string ¼ 0.8 c 2 =GeV 2 and m 2 T is the transverse mass p 2 T þ m 2 . The transverse momentum of a quark-antiquark pair of flavor i is sampled from a Gaussian distribution with the mean The Hadronic interactions with zero net quantum number exchange, and in particular no color exchange between the scattering partners, may leave one or both of the hadrons in an excited state and are referred to as low-mass diffraction. The deexcitation of this state is separated into the resonance region at the lowest masses (M D < 2 GeV), modeled with isotropic phase space decay (thermal fireball), and the continuum region where string fragmentation is used to produce the multiparticle final state. The hadron-Pomeron scattering in high-mass diffraction is approximated by π 0hadron scattering in the rest frame of the diffractive system.

Basic model characteristics
SIBYLL gives a remarkably good description of the general features of hadronic interactions. Particularly encouraging is the comparison of predictions of SIBYLL 2.1 with the results from LHC run I as demonstrated, for example, in Fig. 3 by the yield of charged particles at large scattering angles (pseudorapidity η ∼ tan θ=2). The widening of the distributions is a phase space effect and arises from the available interaction energy. At central rapidities particle production increases with energy as in Fig. 3 according to the growth of the multiple parton scattering probability. The energy dependence of the average number of soft and hard interactions in Fig. 4 shows that below 1 TeV mostly one soft scattering occurs. At higher energies, hard scatterings dominate due to the steep rise of the parton-parton cross section (see σ QCD in Fig. 6). In combination, these figures demonstrate the energy scaling of interaction cross sections, multiple interactions and particle production.
For the high-energy data in Fig. 3, the new model is underestimating the width of the pseudorapidity distribution, indicating a problem with the transition from hard (central) to soft (forward) processes. This problem is becoming more evident with the shift to PDFs in SIBYLL 2.3d that include a steeper rise of the sea quark and gluon distributions toward small x values as favored by measurements at the Hadron-electron ring accelerator (HERA). The scale of the hard scatterings is integrated out for the event generation and the PDFs are evaluated at an effective scale. In nature, the separation between soft and hard scatterings is not well defined and can be thought of as a gradual transition. In principle there should be mixed processes, usually referred to as semihard, which are currently not included in SIBYLL leading to a faster drop of multiplicity for rapidities around the hard-soft scale transition. The comparison to TOTEM measurements in this region (5 < η < 6) reveals a underestimation of the particle density of 30%-40% [54]. However, the more important quantity for EAS than the particle density is the energy flow. Measurements are available in the very forward region by LHCf [55] and at the edge of the central region by CMS and CASTOR [56,57]. The former is described reasonably well by the new model (see Fig. 14 in Sec. II C 1 below), whereas the CASTOR measurement indicates a deficit [56]. The largest part of the energy is carried by particles produced in between these regions and hence remains unobserved. Therefore it is not evident from these data that the omission of semihard processes in the model has an impact on the EAS predictions.

B. Interaction cross section
The parameters of the amplitude are determined by fitting the interaction cross section to measurements.  When the cross section fit was performed for SIBYLL 2.1, the highest-energy data points that were available were the ones obtained at the Tevatron [58-60] (see Table I). These data suffered from an unresolved ambiguity between the measurement by CDF and the other measurements (Fig. 5). The higher data point was supported by some cosmic-ray measurements at the time. Recent measurements at the LHC [61] agree well with each other and suggest a lower cross section. These higher-energy data impose stronger constraints on the extrapolation to UHECR energies constitute an important input in SIBYLL 2.3d.
Despite an overestimation of the interaction cross section, SIBYLL 2.1 gives a remarkably good description of the general features of minimum-bias data. Therefore, we aim for an evolutionary extension of the previous model, in which the hard interaction cross section is smaller. This change yields smaller total and inelastic cross sections in the TeV range and above, while at lower energies remain mostly unaffected according to Fig. 5. Hard parton scattering is calculated in perturbative QCD, generally leaving little room for alterations. The hard cross section can be reduced by increasing the transverse momentum cutoff p min T ðsÞ that defines the transition between soft and hard interactions. However, in SIBYLL the energy dependence is derived from a geometrical saturation condition [see Eq. (5)] and is, therefore, fixed.
A different possibility is the modification of the opacity profile A hard ð ⃗ bÞ. The overlap integral for two protons [the formal definition is given in Eq. (4)] in the model takes the explicit form given by where K 3 ðxÞ is a modified Bessel function of the second kind. The parameter ν h determines the width of the profile that controls the share between more peripheral and central collisions; i.e., narrow profiles lead to a reduction of peripheral collisions. Since most collisions are peripheral, a narrower profile reduces the interaction cross section.   Table I).
In the scattering of waves a refraction pattern is determined by the form of the scattering object. For hadrons, the shape of the refraction pattern in first approximation is described by the elastic slope parameter B ela , the slope of the forward peak of the differential elastic cross section: The −t is the transferred momentum squared. Decreasing the width of the proton profile results in a broadening of the refraction pattern and hence a decrease of the slope. While the interaction cross sections are better described by the narrower profile, the measurements of the elastic slope [65] do not reflect this preference (see Fig. 7). More recent, LHC-constrained parameterizations of the PDFs (e.g., CT14 [72]) instead of the older GRV98-LO [37,38] typically show a less steep rise of the gluon distribution toward small x and hence result in a smaller hard scattering cross section. This would lead to a smaller rise of σ QCD and hence a wider profile can be chosen to reduce the tension with data in B ela . As the integration of the new PDFs in the complete event generator requires the readjustment of almost all model parameters this endeavor is left to a future update.
These modifications to the proton-proton cross sections also affect the cross sections for hadron-nucleus and nucleus-nucleus collisions. The extension to mesonnucleus interactions is discussed in Sec. II F, σ p-air is presented in Fig. 25 and the interaction lengths of iron nuclei, protons, pions and kaons in air are given in Appendix B and discussed in Sec. III.

C. Leading particles
Secondary particles that carry a very large momentum fraction of the initial projectile are called leading particles.
They are of utmost importance for the longitudinal development of EAS since they transport energy more efficiently into the deeper atmosphere requiring at the same time fewer interactions. The origin of leading particles is not clearly related to one hadronic or partonic process and can be thought of as a superposition of all processes contributing to the forward phase space, often involving valence quark interactions.

Leading protons and hadron remnants
In the parton model the leading hadron is related to the partons with the largest momentum fractions, which in most cases are the valence quarks. Figures 8 and 9 FIG. 7. The elastic slope parameter in proton-proton interactions. The slope parameter is related to the width of the impact parameter profile. The decrease in the width of the hard profile between SIBYLL 2.1 and SIBYLL 2.3d means the slope parameter decreases.
. The latter naturally fits into the leading particle definition since in diffraction no quantum numbers are exchanged. The flat region below 0.9 corresponds to the leading particle in nondiffractive events. The presence of this plateau in the proton spectrum and its absence for secondary particles that do not share quantum numbers with the projectile (see antiproton spectrum in Fig. 10) identifies the valence quarks as high momentum constituents of the projectile.
In SIBYLL 2.1 the leading particles are implemented by assigning one of the valence (di)quarks a large momentum fraction. In a proton-proton interaction each proton is split into a quark-diquark pair forming a pair of strings between a quark and a diquark of the other proton, as illustrated in Fig. 11(a). The momentum fraction of the quark is sampled from a soft distribution as in Eq. (9), leaving a larger fraction to the diquark. In addition, the subsequent fragmentation of the quark-diquark string is biased toward the diquark by sampling the energy fraction in the first string break next to the diquark from f lead ðzÞ ∼ z instead of the standard Lund function [Eq. (11)]. This mechanism reproduces the observed flat proton spectra in Figs. 8-10. Interactions of hadrons at low energies [e.g., ffiffi ffi s p ¼ Oð20 GeVÞ] are dominated by soft parton scattering. In SIBYLL, most of these interactions happen between the valence quarks (see Fig. 4). The conservation of energy and baryon number for such systems introduces a strong correlation between the production of leading protons and central (x F ∼ 0) antiprotons, as both come from the hadronization of the same valence quark system. In the leading proton scenario, where a large momentum fraction is assigned to the leading string break, an antiproton produced in a later break is necessarily slow. Often its production will be energetically forbidden because the antiproton has to be produced alongside a second baryon. The opposite case, in which the leading proton is slow [f Lund ðzÞ ∼ exp ð−1=zÞ ≪ f lead ðzÞ ∼ z as z → 0], is more problematic since the antiproton can carry a large momentum fraction. Measurements of x F spectra of protons and antiprotons in Fig. 10 do not confirm the presence of antiprotons with large momenta (an additional discussion of baryon-pair production can be found in Sec. II D). By changing the momentum fraction of the leading protons the production of antiprotons with large momentum fraction cannot be avoided since the protons demonstrate a flat spectrum down to the central region.
In SIBYLL 2.3d the issues with leading baryon production are addressed with the so-called remnant formation. In this mechanism, the leading protons are produced from the remnant, while antiprotons and central protons are produced from strings that are attached to soft sea quarks [Figs. 11(b) and 11(c)]. The momentum fraction of the sea quarks is sampled from f soft q ðxÞ ¼ ð1 − xÞ 1.5 ðx 2 − m 2 q =sÞ −1=4 with m q ¼ 0.6 GeV. The momentum fraction for the remnant (system of valence quarks) is distributed like x 1.5 .  [75]. The flat distribution for protons is achieved using an ad hoc mechanism in SIBYLL 2.1. However, the central and the leading particles are produced by the same process, resulting in a hard spectrum for antiprotons. In SIBYLL 2.3d the central and the fragmentation region are related to separate processes, leading to more accurate descriptions for longitudinal baryon spectra. The energy and the momentum transferred in the remnant interaction are modeled similarly to diffractive interactions as discussed at the end of Sec. II A 1. The squared mass spectrum approximately follows dN=dM 2 r ∼ 1=M 2 r and the slope of the p T spectrum is with the parameters B 0;r ¼0.2GeV 2 =c 4 , a r ¼ 7.0 GeV 2 =c 4 and b r ¼ −2.5 GeV 2 =c 4 . In addition to the continuous spectrum, discrete excitations of resonances are included. Due to their isospin structure, the decay channels may be weighted differently than for isotropic phase space decay. For each projectile two resonances are included (e.g., see Table II). When parton densities become large at high energies and the number of parton interactions increases, it is less likely that partons remain to form a remnant. In this case the situation is more similar to the two-string approach in SIBYLL 2.1. This transition effect is taken into account by imposing a dependence on the sum of soft and hard parton interactions ðn s þ n h Þ to the remnant survival probability In nuclear interactions (even at low energies) parton densities can be large. Correspondingly, the remnant probability depends on the number of nucleon interactions N w . The relative importance of nucleon and parton multiplicity is determined by ϵ and is set to 0.2. The remnant survival probability at low energies P r;0 is 60%. The spectrum of the remnant excitation masses for proton interactions in Fig. 12 demonstrates how different hadronization mechanisms apply for different regions of the mass spectrum. For large masses (ΔM ¼M remnant −m projectile > 1GeV, where m projectile is the mass of the projectile), indicating the presence of a fast valence quark, the deexcitation is very anisotropic and particles are emitted mostly in the direction of the leading quark. In this case, the hadronization of high-mass remnants is implemented as the fragmentation of a single string. At intermediate masses (0.4 GeV < ΔM < 1 GeV), a continuum of isotropic particles is produced by phase space decay. The number of particles produced is selected from a truncated Gaussian distribution with the mean n thermal ¼ 2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ΔM=GeV p , n thermal > 2. Below the threshold for the production of particles and resonances (ΔM < 0.2 GeV), the remnant is recombined to the initial beam particle. This recombination region determines the proton distribution at intermediate and large Feynman x. Hence, the shape of the final particle spectra depends on the combination of the separate hadronization mechanisms. The adjustment of the remnant model parameters has been mainly achieved from comparisons with the leading low-energy proton data shown in Fig. 8 together with the antiprotons distribution shown in Fig. 10. In particular, the latter is much better described by the updated model. At the higher energies probed in the ZEUS experiment [74] (see Fig. 9), the contribution from the remnant in the region x F > 0.9 overlaps with the diffractive peak, resulting in an overestimation of the spectrum in SIBYLL 2.3d, while in the region of 0.6 < x F < 0.8 the spectrum is underestimated. This can be addressed in the future by adjusting the remnant and the diffractive mass distribution.
Another drawback of the model for leading particle production in SIBYLL 2.1 is the insufficient attenuation of the leading particles in the transition from proton to nuclear targets (see secondary proton spectrum in Fig. 13). While the proton spectrum is clearly affected by the number of target nucleons, this effect is much smaller for mesons (pions). The model for the reduced remnant formation probability in the presence of multiple target nucleons [Eq. (15)] in SIBYLL 2.3d reproduces this effect correctly.
The model parameters are adjusted according to low-energy data from the NA49 experiment that provides a large x F coverage. However, the remnant model affects high energies as well, resulting in a significant improvement of leading neutrons at LHCf [79] (7 TeV), as shown in Fig. 14.

Leading mesons and ρ 0 production
A second important role of leading particles in EAS is their impact on the redistribution of energy between the hadronic and the electromagnetic (EM) shower component. Any charged pion of the hadronic cascade can transform into a neutral pion in a charge exchange interaction. Through the prompt decay of the neutral pion into two photons, all the energy is then transferred to the EM component The influence of this reaction is largest for the leading particles and usually results in a decrease of the muon production that occurs at late stages of the EAS development [80]. A suppression of the pion charge exchange process has the opposite effect.
An example for such a competing reaction is the production of neutral vector mesons [ρ 0 ∶IðJ CP Þ¼1ð1 −− Þ] from a pion beam: Whereas a neutral pion decays into two photons, the conservation of spin requires a ρ 0 to decay into two charged pions.
In the Heitler-Matthews model [81] the average number of muons in an EAS initiated by a primary cosmic ray with energy E 0 is given by and critical energy E c . The change of the number of muons per decade of energy (α) thus depends on the total and charged multiplicities. It is evident that the ratio between ρ 0 and π 0 production directly affects the exponent α.
In charged pion-proton interactions the NA22 fixed target experiment found that at large momentum fractions vector mesons are more abundantly produced than neutral pions (Fig. 15) [82,83]. In the dual parton approach with standard string fragmentation, as it is used in SIBYLL and several other models, this result is unexpected and probably cannot be reproduced without invoking an additional exchange reaction. Recent measurements by the NA61 Collaboration have confirmed the leading ρ 0 enhancement in case of pion nuclear interactions [84].
The leading ρ enhancement and π 0 suppression can be reproduced in SIBYLL by adjusting the hadronization for the remnant and for diffraction dissociation. The result is shown in Fig. 15. The transition from proton to nuclear targets is entirely described by the dependence of the remnant survival probability on N w in Eq. (15). As demonstrated in Fig. 16, the softening of the leading ρ 0 spectrum in pion-carbon interactions is well reproduced by the current model. The intersection between the ρ 0 and π 0 spectra is predicted to occur at the same x F in pion-proton and pion-carbon collisions (x F ≈ 0.5). The position of this intersection is important for EAS since it determines the fraction of the energy that goes either into the EM or hadronic shower component. Until the spectrum of π 0 is measured for meson-nucleus interactions, this intersection is experimentally not fully determined. Thus the total effect of the leading ρ 0 on the number of muons in EAS remains unconstrained (this topic is further discussed in Sec. III C).

Baryon-pair production
While the importance of leading particles for the development of EAS is clear, it is not directly evident how a relatively rare process as baryon-pair production affects muon production [80,85,86]. The role the baryons play is similar to a catalyst in a chemical reaction. Any baryon produced in an air shower will undergo interactions and produce new particles; in particular, it will regenerate at least itself due to the conservation of baryon number. The interactions continue until the kinetic energy falls below the particle production threshold. Through this mechanism any additional baryon yields more pions and kaons and hence ultimately more muons. In terms of the Heitler-Matthews model, where the number of muons is given by Eq. (18), additional baryons represent an increase of the exponent α.
In SIBYLL's string model, baryonic pairs are generated through the occurrence of diquark pairs in the string splitting with a certain probability, which in SIBYLL 2.1 is the global diquark rate P diq =P q ¼ 0.04. This model works well at low energies where mostly a single gluon exchange occurs. It fails, however, in the multiminijet regime at higher energies [87,88] (see Fig. 17). Both regimes can be jointly described by choosing a different value for the diquark pair rate in events with multiple parton interactions. The constant ratio of baryons to mesons in the measurement that is shown in Fig. 18 suggests that baryon-pair production cannot depend on the number of minijets or the centrality of the interaction. In the model, the diquark probability is then FIG. 15. Feynman-x spectrum of neutral pions and their spin-1 resonance state ρ 0 in π þ -proton collisions at p lab ¼ 250 GeV=c [82,83]. The expectation from standard quark splitting (π þ ∶ud) and fragmentation is that a fixed fraction of the leading π þ transforms into neutral pions [π 0 ∶ðuū − ddÞ= ffiffi ffi 2 p ] and a smaller fraction into the resonance state ρ 0 (upper figure). Data, on the other hand, show an enhancement of the production of the resonant state and a suppression of the ground state in the region of the leading particle. The effect is reproduced in SIBYLL 2.3d (lower figure) by increasing the rate at which resonances occur in the fragmentation of diffractive processes and by including the ρ 0 as a resonance state in the remnant formation of the pion.

FIG. 16.
Feynman-x spectrum of neutral ρ mesons in pioncarbon interactions as measured in the NA61 experiment [84]. This measurement confirms the enhancement of leading ρ 0 for nuclear targets. Compared to the data obtained with a proton target (gray triangles), the carbon data (blue squares) reveal a softening of the spectrum, indicating the relevance of interactions with multiple target nucleons. The new remnant model (bottom) correctly reproduces the softening of the leading ρ 0 and predicts a suppression of the production of leading neutral pions (red curve).
where P single ¼ 0.06, P multi ¼ 0.13, P diff ¼ 0.04 and n s þ n h is the sum of the number of soft and hard interactions. This purely phenomenological model is inspired by the observation in e þ e − collision experiments where it is found that baryon-pair production in the fragmentation of quarks or gluons can be different [92].

Transverse momentum
The transverse momentum in the string fragmentation model (string p T ) is usually derived from the tunneling of the quark pairs in the string splitting, which results in a Gaussian distribution [93]. However, the observed distribution of transverse momenta in hadron collisions [94,95] more closely resembles an exponential distribution as predicted by models of "thermal" particle production [96], motivating us to distribute the string p T in SIBYLL 2.3d according to where i denotes different flavors of quarks and diquarks. The energy dependence of the average transverse mass hm T;i i is parameterized as with the parameters A T;i and m 0;i . The values are given in Table III. These values are derived from the measured p T spectra of pions, kaons and protons at low (NA49) and high energies (CMS; see Fig. 19). In addition to the string p T , the hadrons acquire their transverse momentum from the initial partonic interaction. As previously mentioned, the parton kinematics in SIBYLL 2.3d are determined from post-HERA PDFs (GRV98-LO [37,38]), which predict a steeper rise of the gluon density at small x, when compared to the old parameterization in SIBYLL 2.1. With the new parameterizations the transition between the regions dominated by soft scattering (p T < 3 GeV) and hard scattering is described better (see Fig. 20).
While the new PDFs help in describing the transition region, the rise of the average transverse momentum with energy is not described well (not shown). To account for the rapid rise with energy seen in the data (see Fig. 21), the energy dependence of the average transverse mass in Eq. (19) is set to be quadratic in logð ffiffi ffi s p Þ. The integration of post-LHC PDFs, in which the small-x gluon densities tend to be smaller than in the GRV98 parameterizations, is not expected to help with this.  [90]. The central multiplicity is sensitive to the number of parton interactions. High-multiplicity data suggest a constant rate of baryon production per minijet, whereas the region at low central multiplicities is populated by diffractive events and events with a single parton interaction. The substructure for SIBYLL 2.3d is due to the remnant model.

E. Nuclear diffraction and inelastic screening
Nuclear cross sections in SIBYLL 2.1 are calculated with the Glauber model [14,30] neglecting screening effects due to inelastic intermediate states [98] in which an excited nucleon may reinteract and return to its ground state. Also, diffraction dissociation in hadron-nucleus interactions is restricted to the incoherent component. SIBYLL 2.3d also makes use of the Glauber model but includes screening and the diffractive excitation of the beam hadron in a coherent interaction [99,100].
In analogy to diffraction dissociation in hadron-nucleon interactions [16,43], the coherent diffractive excitation of a hadron by a nucleus is implemented using a two-channel formalism with a single effective diffractive intermediate state, where the shape of the transition amplitude to the excited state is equal to the elastic amplitude. The remaining free parameter of the model is the coupling between the states λ. In the following, we will limit the discussion to proton-nucleus interactions and substitute the nucleon with a proton. With where jpi represents the proton and jp ⋆ i is the effective intermediate state or diffractive final state, the generalized amplitude for the described model of proton-proton interactions isΓ The proton-nucleus cross sections σ pA are calculated with the standard Glauber expressions using the proton-proton amplitudeΓ pp , projected onto the desired transition hpj Á Á Á jpi. The diffractive cross sections are calculated in the same way but for the projection hp ⋆ j Á Á Á jpi [99,101]. The assumed equivalence of the elastic and diffractive amplitude (Γ pp→p ⋆ p ¼ λΓ pp→pp ) implies for the energy dependence of the coupling λ where M 2 D;max is the upper limit for the excitation mass in diffraction dissociation motivated by the coherence limit [41] and s is the square of the center-of-mass energy. We assume the coupling λðsÞ to be universal for different hadrons. The cross sections in Eq. (22) are taken from FIG. 21. Average transverse momentum of charged hadrons as a function of center-of-mass energy [11,51,97]. The low-energy limit is given by the confinement of the partons to the hadron. The increase with energy is in part due to the increase in the hard scattering (jets) threshold (p min T ) and in part due to the hardening of the string-p T spectrum according to Eq. (19). While the rise in the p T cut is given by QCD and saturation, the rise of string p T is entirely phenomenological. parameterizations [102,103]. The single diffractive cross section used in proton-proton collisions and the parametrization of the coupling λðsÞ are shown in Fig. 22. The difference is due to the larger value for the upper mass limit of M 2 D;max =s ¼ 0.1 for hadron targets, whereas a lower value of M 2 D;max =s ¼ 0.02 was found to give the best description of the production cross sections in proton-carbon and neutron-carbon interactions [104][105][106]. Although the description of data in Fig. 22 does not look ideal, one shall consider that several shown data points are extrapolations of rapidity gap data from limited detector acceptance and must not represent accurately σ SD . A more accurate description of rapidity gaps [107,108] and particle production with diffractive cuts [109] has to be addressed in a different revision of the model.
The cross section for the diffractive dissociation of the projectile proton in proton-carbon interactions is shown together with the predictions from commonly used interaction models in Fig. 23. The diffractive cross in SIBYLL 2.1 section drops toward high energies, whereas the contribution from coherent diffraction in SIBYLL 2.3d compensates this trend. QGSJETII-04 [28] and EPOS-LHC [29] predict almost constant cross sections. Since the diffractive cross section is small relative to the production cross section of Oð400 mbÞ, the differences among the models are not expected to be important in EAS.

F. Meson-nucleus interactions
The extension of the model from proton-nucleon collisions (as discussed Sec. II B) to pion-and kaon-nucleon collisions is straightforward, since at the microscopic level the interactions are treated universally as scatterings of FIG. 24. Secondary particle spectra in pion-carbon interactions with p Lab ¼ 350 GeV measured by NA61 [114,115] and shown together with predictions from SIBYLL 2.3d (full line) and SIBYLL 2.1 (dashed line). Note that these newer data were not yet available during the development of the models. Some aspects of the distributions are better described by the newer model, in particular the antiprotons; SIBYLL 2.3d is far from perfect. Although SIBYLL 2.3d lacks forward kaons, the description of the charge ratio is improved, resulting in a positive impact on the atmospheric muon charge ratio. quarks and gluons. Differences, in particular at low energies, arise from the different profile functions [112], momentum distributions (PDFs) [113] and Regge couplings in the soft interaction cross section (see Appendix C).
Since the measurements [114,115] from Fig. 24 were not yet available during the development of the model, the distributions obtained with SIBYLL 2.3d and SIBYLL 2.1 are predictions. Some improvement is observed in the distributions of baryons and kaons. However, the production of central pions, forward kaons and antiprotons clearly demonstrates that the model requires more work.

III. AIR-SHOWER PREDICTIONS
Some relations between air-shower observables and specific properties of hadronic interactions have been studied in the past [116]. Here we focus on the depth of shower maximum hX max i and the number of muons N μ . The calculations are obtained with CONEX [117], using FLUKA [118,119] to simulate interactions at E kin < 80 GeV. The employed scheme is hybrid, meaning that all subshowers with less than 1% of the primary energy are treated semianalytically using numerical solutions of the average subshower. We compare the predictions from SIBYLL 2.3d with the previous SIBYLL 2. 1 and two other post-LHC models, EPOS-LHC [29] and QGSJETII-04 [28]. In addition, we calculate some of the observables with modified versions of SIBYLL 2.3d to show the impact of individual extensions introduced in Sec. II. The extensions are labeled in Table IV and will be used throughout the next sections. Tables with the predictions for hX max i, N μ and λ int can be found in Appendix B.

A. Interaction length and σ air
The simplest and most direct connection between the development of an air shower and hadronic interactions is governed by the interaction length λ int ðEÞ ¼ hm air i= σ prod ðEÞ. It determines the position of the first interaction in the atmosphere and thus directly influences the position of the shower maximum (X max ). In the Glauber model [30], the inelastic cross section in proton-air interactions, σ prod , is derived from the proton-proton cross section σ pp . A smaller σ pp , as in SIBYLL 2.3d (Sec. II B), translates into a smaller proton-air cross section. The effect on σ prod is less than proportional since σ pp is only a small contribution to the overall value that is mostly defined by the nuclear geometry. An additional small reduction of the cross section originates from inelastic screening (Sec. II E). The updated proton-air cross section results in a better compatibility with observations as can be seen in Fig. 25. The impact of the updated interaction length on hX max i is demonstrated in Fig. 29. The reduction of the cross section at high energy leads to a shift of 5-10 g=cm 2 . Interaction lengths for different primary nuclei and secondary mesons in air are listed in the Appendix.
B. hX max i and σðX max Þ The depth at which an individual shower reaches the maximum number of particles is determined by the depth of the first interaction and the subsequent development of the particle cascade. In very general terms, the development of the cascade is influenced by how the energy of the interacting particle is distributed among the secondaries, in particular by how energy is shared among electromagnetic and hadronic particles. The average shower maximum for proton-initiated showers in SIBYLL 2.3d is almost 20 g=cm 2 deeper than that in SIBYLL [120][121][122][123][124][125][126][127]. The reduction between the versions of SIBYLL comes mainly from the updated proton-proton cross section, whereas the correction due to inelastic screening is small. The most precise measurement at the highest energies by the Pierre Auger Observatory also favors a lower cross section [69,127] in agreement with the extrapolations of the LHC measurements.
protons in air. Another contribution to the difference in hX max i is the decreased inelasticity of the interactions (see Fig. 28). Figure 29 illustrates the effect of the individual modifications on the shift in hX max i. This comparison is produced by individually switching off the model extensions introduced in Sec. II and summarized in Table IV. The change in the interaction length (cyan line) is responsible for 10 g=cm 2 out of the 20 g=cm 2 difference between SIBYLL 2.1 and SIBYLL 2.3d at high energy. Coherent diffraction on the nuclei in the air (purple line) contributes another 5 g=cm 2 . The remaining 7 g=cm 2 cannot be attributed to a single feature and emerge from the combination of the model modifications.
The enhanced ρ 0 production (green line) and the improved baryon-pair production (not shown) have a small effect on hX max i. These processes mostly affect the later stages of EAS that are more important for muon production (see the next section for more details).
The overall effect of the changes in the multiparticle production between the 2.1 and 2.3d versions result in a decreased inelasticity in Fig. 28 for proton and pion interactions. Compared to SIBYLL 2.1, the inelasticity increases less steeply with energy and should have impacted the elongation rate for protons. This effect seems to have been compensated by the change in the energy dependence of the interaction length or cross section (cyan line in Fig. 29).
FIG. 26. Average depth of air-shower maxima hX max i for different models compared to recent data from the Pierre Auger Observatory [128,129] obtained with the fluorescence detectors. The model lines represent the expectations for a pure proton and iron composition, respectively. The deviation of the data from the pure composition indicates a change toward a mixed composition; i.e., cosmic rays consist of a combination of light and heavier nuclei. The modifications in SIBYLL 2.3d drive the interpretation toward heavier nuclei since the hX max i becomes deeper.   Table IV. The change of the cross section for coherent diffraction as described in Sec. II E increases the hX max i by 5 g=cm 2 . The change in λ int;p due to the smaller proton-proton cross section amounts to another 10 g=cm 2 . ρ 0 production has a negligible effect on hX max i.
The separation between proton and iron showers in hX max i at lower energies is larger in SIBYLL 2.3d (see Fig. 30), since coherent diffraction only deepens the proton showers and has no effect for nuclear projectiles. This effect is expected to have a higher impact on the measurements of the cosmic-ray composition that were previously interpreted using predictions from SIBYLL 2.1.
The width of the distribution of shower maxima σðX max Þ in Fig. 31 increased by 10 g=cm 2 between the versions, becoming the largest of all CR models. This change is dominated by the increased interaction length, as is shown Fig. 32. Note that the σðX max Þ increases only for protons, widening the distance between the pure protons and other masses. This behavior has an important impact on the theoretical interpretation of the measurements in terms of cosmic-ray sources and it has been shown that SIBYLL 2.3d produces distinctly different results compared to other contemporary interaction models [130].

Number of muons
In recent years it became evident that the muon content observed in air showers differs from the predictions of the interaction models [131]. Recently the Pierre Auger Observatory quantified this "muon excess" at ground to be at the order of 30%-60% [19]. This result is in agreement with the numbers obtained by the Telescope Array [132]. In contrast to the hX max i, the production of muons is very sensitive to hadronic particle production at all stages of the shower. It is therefore legitimate to attribute the muon excess to a combination of flaws in the modeling of hadronic interactions. Alternatively, the excess could also be seen as the signature of a new physical phenomena beyond the scales probed by current colliders [133,134].
Most muons in EAS originate from decays of hadrons, most abundantly of pions and kaons. Due to their relatively long lifetime, especially at high energy, these mesons reinteract with air molecules and initiate additional cascades, copiously creating more mesons. The large dependence of the number of muons N μ on hadronic interactions can be understood by considering that any flaw in the production spectrum of secondaries that persists across multiple generations of reinteractions has a multiplicative effect at the final stages of the shower. In fact, most muons are produced at the end of the cascade where the energies of mesons are low enough to allow a significant fraction to decay before the next interaction. This cascade process leads to a power law relation between the number of muons and the primary energy as shown in Fig. 33 and by Eq. (18). The slope corresponds to the exponent α that depends on the fraction of hadrons that effectively participate in the production of muons. The FIG. 30. Difference in the hX max i between proton-and ironinduced showers. This observable is relevant for measurements of the cosmic-ray mass composition that are based on observations of hX max i.
FIG. 31. The width of the X max distribution expected from models using a pure composition compared to data from the Pierre Auger Observatory [128,129]. The σðX max Þ plays an important role in the determination of the mixture of different mass groups at a particular energy.  Table IV. The fluctuations are most strongly affected by the change in the interaction length. Since the nuclear cross sections are not very sensitive to changes of σ pp , the impact is highest for proton primaries. This is clearly seen for the iron predictions in Fig. 31. enhanced baryon-pair and leading ρ 0 production in SIBYLL 2.3d result in a higher number of charged pions and hence a higher value of α. Relative to SIBYLL 2.1 (see Fig. 34) the new version has at least 30% more muons at PeV energies, which increases to ∼60% at the highest energies due to a steeper slope. The other post-LHC models include similar extensions and therefore show the same behavior in the muon number.
The influence of baryon-pair production and ρ production on the number of muons is shown in Fig. 35, from which the contribution from each enhancement can be seen individually. A reduction of the baryon-pair production to the level of SIBYLL 2.1 results in only 10% less muons at ground. As discussed in Sec. II C 2, the ratio between ρ 0 and π 0 is more important for muon production. This is confirmed by Fig. 35 where the difference is at the level of 25%. With such large variations to the observable number of muons induced by qualitative improvements to the physics of the model, in contrast to just parameter settings, it appears likely that the muon excess in UHECR interactions originates from the shortcomings of the current hadronic interaction models.

Muon energy spectrum
The energy spectra of muons for the post-LHC interaction models relative to SIBYLL 2.1 are shown in Fig. 36. The clear rise in the number of low-energy muons predominantly originates from the increased number of cascading hadrons due to the modified baryon-pair and ρ production. The enhancement of muons at high energies originates from decays of charmed hadrons which are an exclusive feature of SIBYLL 2.3d in current air-shower simulations. The number of these, so-called, prompt muons is very low and hence no impact is expected for air-shower observations since experimentally an energy threshold around a few PeV is required. Muons with an energy in excess of 1 TeV (100 TeV) constitute only 0.1% (3.1 × 10 −5 %) of all muons at ground for a 10 19 eV shower (see also Appendix B). For inclusive lepton fluxes this contribution has important implications as discussed in Ref. [25].
In the left panel of Fig. 36 the energy and incident angle of the primary CR resemble the typical experimental conditions of IceTop and IceCube [135,136], whereas the right panel resembles typical conditions at the Pierre Auger Observatory [18]. It is remarkable that the modelspecific features of the spectrum are present across very different primary energies.
Another observation is that the current models predict different shapes of the muon spectrum. With a combination of the surface air-shower array IceTop and the main instrumented IceCube volume deep in the Antarctic ice, the IceCube Observatory has the potential to discriminate among the interaction models by measuring the muon  Table IV).
content of a single air shower at two different energy regimes simultaneously. IceTop is sensitive to the lowenergy muons while only the muons with E μ ∼ TeV can penetrate the ice deep enough to generate the "in-ice" muon signal. The preliminary results clearly indicate that SIBYLL 2.1 has too many high-and too few low-energy muons [137]. The discrepancy is expected from the discussion of Fig. 36 above, since SIBYLL 2.1 neither describes the baryon-pair production nor the ρ production very well. The same analysis shows that SIBYLL 2.3d accurately reproduces both low-and high-energy muons. The result is, however, difficult to translate into constraints on the hadronic parameters since the (unknown) mass composition has to be simultaneously taken into account. The impact of each modification on the muon spectrum is illustrated in Fig. 37. According to the figure baryon-pair production contributes dominantly at low energies, while the contribution from ρ affects all energies.

Effect of the projectile mass on muon production
The spectra for the individual mass groups of cosmic-ray nuclei are not well known across the entire energy range of the indirect air-shower measurements [138]. The main source of this systematic uncertainty stems from ambiguities among the interpretations of EAS observables with different hadronic interaction models. At present, at ultrahigh energies the most robust method to estimate the composition relies on the electromagnetic component only. Recent attempts to use the surface detector and exploit the muon content as a sensitive variable often result in incompatible results [139].
We study the ratio of the muon energy spectra for the two extreme composition assumptions, pure protons and pure iron. The ratios in Fig. 39 demonstrate that the difference in the number of GeV muons is small between UHE protons and iron nuclei (∼20%-40%). As discussed in the previous section, similar variations are expected just from swapping the interaction model. At higher muon energies (E μ > 100 GeV) protons and iron are well separated. The shape comes from two effects: the earlier development of iron showers due to the shorter interaction length of the primary nucleus and the lower energy carried by the individual nucleons in the iron nucleus. If one would take the muon energy spectrum from iron primaries with E Fe ¼ 56E p and compare with the spectrum in proton showers at the shower maximum, they would have identical shapes.
The superposition ansatz (E 0 → E 0 =A and N A μ ¼ AN 1 μ ) in the Heitler-Matthews model of Eq. (18) yields for the composition dependence of the total muon number an additional multiplicative term ð1 − αÞ lnðAÞ. If α approaches unity, as is the case for the current model  [135]. Right: Showers at 10 EeV are simulated with a zenith angle of 67°as they are observed at the Pierre Auger Observatory [18]. The increased number of PeV muons in SIBYLL 2.3d is due to the prompt decay of charmed hadrons not present in any of the other models [23,88].
FIG. 37. Ratio of the muon energy spectrum between the versions of SIBYLL 2.3d and SIBYLL 2.1 for 10 EeV proton showers. The models labeled "off" refer to modified versions of SIBYLL 2.3d where the extensions for enhanced ρ 0 and baryon production have been switched off (see Table IV). Baryon-pair production enhances mostly the number of low-energy muons, while ρ 0 production also affects high-energy muons. extensions, the difference between protons and nuclei decreases. This expectation is confirmed by full model calculations in Fig. 38, in which the muon number varies by only 35% between proton and iron for post-LHC models, while for SIBYLL 2.1 the difference is almost 50%. However, the ratio of iron to proton spectra from different interaction models agree remarkably well (see Fig. 39).
The influence of individual model processes on the separation between proton and iron are demonstrated in Fig. 40. Both baryon-pair production and ρ production enhance low-energy muons and essentially reduce this separation through a more elongated hadronic cascade (or, in other terms, a larger α in the Heitler-Matthews model). However there are subtle differences. At 10 16 eV only enhanced ρ production is important for the difference between the primaries in TeV muons, while low-energy muons are affected by both mechanisms. At 10 19 eV, the difference between primaries is not much affected by ρ production and baryon-pair production and other changes in the model seem to play more central roles.

IV. DISCUSSION AND CONCLUSION
This paper documents the latest extensions to the hadronic interaction model SIBYLL and discusses their impact on extensive air showers. The model update is motivated through the availability of recent particle accelerator measurements, where measurements from experiments at the LHC and those from fixed-target experiments are equally important. The goal is to improve the consistency in the description of extensive air showers, in particular related to the muon content that impacts the interpretation of the mass composition of the primary cosmic rays. A tabulated overview of the changes between the SIBYLL 2.1 and SIBYLL 2.3d is available in Appendix C.
The interaction cross sections from measurements at the LHC point toward lower total and inelastic proton-proton cross sections that favor the low data points from measurements at the Tevatron. Our new fits take the measurements up to ffiffi ffi s p ¼ 13 TeV into account, reducing the extrapolation uncertainties up to ultrahigh cosmic-ray energies. The effect on the proton-air cross section is a reduction of the tension between SIBYLL and the cross section measurement derived from UHECR observations at the Pierre Auger Observatory. The spectra of identified particles, measured in central phase space at the LHC, allow us to adjust the hadronization to account for a higher baryon-pair production compared to the previous version.
Together with the updated PDFs, the high-energy data constrain the shape and energy dependence of transverse momentum distributions.
On the other hand, the fixed-target measurements in p-p, p-C, π-p and π-C beam configurations yield enough information to identify the shortcomings of the previous model version and entirely revise the leading particle production. We implement a model that makes use of the remaining hadron content in the beam remnants that can undergo further excitation and hadronization processes. This mechanism adds necessary degrees of freedom to decouple very forward particle production from central.
None of the new features requires drastic changes in the underlying principles and assumptions that were defining SIBYLL during the last decades. Microscopically, the main picture is still a combination of the dual parton and the minijet model, a fusion of perturbative QCD (hard component) and elements of the Gribov-Regge field theory (soft component).
We identified, however, a number of problems that indicate a necessity to depart from these well-explored principles in future versions. One of these problems is related to the growth of the multiplicity distribution that rises faster in the model than in data. A second problem is the narrow width of the pseudorapidity distributions that most likely is an effect of the missing contribution from semihard processes. Both aspects are related to the underlying partonic picture, and a permanent solution will require an overhaul of several old principles in the code base.
On the nuclear side, the previous Glauber-based model is extended to include screening corrections on the production cross section due to inelastic intermediate states. The updated model for diffraction dissociation now incorporates the process of coherent diffraction, in which the beam hadron transitions to an excited state without the target side nucleus losing its coherence.
Charm hadron production is added explicitly for particle astrophysics applications. In particular this affects calculations of atmospheric neutrinos at very high energies, where the flux of atmospheric leptons competes with that of astrophysical origin. The details of this topic are discussed in a separate publication [25].
Regarding air showers, several of the changes to the hadronic interaction model impact the simulations. The showers reach their maximum deeper by 20 g=cm 2 with respect to SIBYLL 2.1, mainly due to the modifications to nuclear diffraction and the updated interaction cross sections for protons and pions. The fluctuations of the X max in proton showers are almost 10 g=cm 2 larger as an effect of the increased interaction length and elasticity. Both modifications are likely to yield a notably heavier composition in the interpretation of the flux of UHECR.
The muon number in SIBYLL 2.3d drastically increases by 20%-50% relative to SIBYLL 2.1, which was previously known to yield too few muons. Compared to the other interaction models the new version has the highest number of muons but only exceeding the numbers from EPOS-LHC and QGSJETII-04 by ∼1%-5%. This change will certainly reduce the muon excess seen by the Pierre Auger Observatory and the Telescope Array but will most likely not be sufficient to remove entirely the tension between simulation and data. We demonstrated that the forward spectrum of π 0 and leading ρ mesons in π-nucleus interactions effectively modulates the total muon number and that a constraining measurement of the π 0 is one of the leading uncertainties.
We expect that the combined measurements with the IceCube and IceTop detectors at two energy regimes, and the event-by-event composition sensitivity of the upgrade of the Pierre Auger Observatory (AugerPrime) [140], will help to resolve the mysteries around the muon component in EAS.  Table V) have become publicly available throughout the development cycle and were available in the air-shower simulator CORSIKA (versions 7 and 8) [141,142] and the cascade equation code MCEq [25]. In this section, we give a brief overview of the changes and estimate the quantitative impact on hX max i and the number of muons in air showers. SIBYLL 2.1 is the basic implementation of the hadron interaction model and was outlined in Sec. II A and described in detail elsewhere [16]. The first public release of the SIBYLL 2.3 [26] model improved the compatibility with LHC measurements and astroparticle experiments as described in the main text. The model exhibited a stronger violation of Feynman scaling in the fragmentation region than supported by data [25,27] that has been addressed in SIBYLL 2.3c.
In a recent publication the behavior of the π AE to π 0 ratio in different mechanisms of hadronization and the role in muon production in air showers were discussed [143]. In SIBYLL 2.3c this ratio has a stronger than expected energy dependence (see left panel of Fig. 41), because a part of the model responsible for the leading ρ 0 (Sec. II C 2) interfered with the fragmentation of minijets.
Although this behavior increases the number of muons in air showers and reduces the tension with the observations, it was unintended and has been addressed in SIBYLL 2.3d. The maximal effect occurs in the central phase space but as shown by the distribution of charged particles over pseudorapidity in the right panel in Fig. 41, the impact is small.
In general the different versions of SIBYLL 2.3 have rather small effects on air-shower observables. The differences in hX max i for proton-induced showers is shown in the left panel in Fig. 42, which are up to 5 g=cm 2 at high energies. The muon number at 10 19 eV is ≈7% smaller in SIBYLL 2.3d than in SIBYLL 2.3c (right panel in Fig. 42). We verified that these two releases have almost identical inclusive lepton fluxes (as in [25]).  TABLE VI. Predictions for the depth of shower maximum, the fluctuations thereof and the number of muons in SIBYLL 2.3d for protonand iron-induced showers. X max is calculated by fitting a parabola to the profile of energy deposit in the atmosphere. The number of muons is taken at a depth of 2030 g=cm 2 , counting all muons with an energy exceeding 1 GeV. Showers were simulated with an inclination of 67°using CONEX hybrid simulations [117]. The Monte Carlo to cascade threshold was set to E thr =E 0 ¼ 10 −2 . Hard minijets Leading-order QCD with energy-dependent p T threshold PDF: cross section GRV-98LO [37,38] GRV-98LO PDF: sampling Eichten et al.