The hadronic interaction model Sibyll-2.3c and inclusive lepton fluxes

Muons and neutrinos from cosmic ray interactions in the atmosphere originate from decays of mesons in air-showers. Sibyll-2.3c aims to give a precise description of hadronic interactions in the relevant phase space for conventional and prompt leptons in light of new accelerator data, including that from the LHC. Sibyll is designed primarily as an event generator for use in simulation of extensive air showers. Because it has been tuned for forward physics as well as the central region, it can also be used to calculate inclusive fluxes. The purpose of this paper is to describe the use of Sibyll-2.3c for calculation of fluxes of atmospheric leptons.


I. INTRODUCTION
The main theme of this paper is the connection between hadronic interactions at high energies and the inclusive spectra of atmospheric leptons.The coupled transport equations that relate the lepton spectra to the primary cosmic-ray spectrum depend on the properties of the hadronic interactions, implemented here with Sibyll-2.3c.A brief introduction to these transport equations is given in II.The numerical methods are described in section III.Section IV establishes the connection between particle physics observables and the regions of phase space that are important for inclusive lepton observables.A key observation is that, because of the steep primary cosmic-ray spectrum, it is the forward fragmentation region of hadronic interactions that is of special importance for inclusive lepton spectra.In the second part of the paper (Section V) we describe how Sibyll-2.3cdeals with the forward fragmentation region including production of charm.In Section VI we summarize the impact of Sibyll-2.3con inclusive lepton spectra.We compare with the corresponding results of its predecessor, Sibyll2.1, and with several other event generators in current use.We try, as far as possible, to relate observed differences to specific features of hadronic interactions with the idea that precise measurements of atmospheric lepton spectra have the potential to constrain features of forward production of hadrons at high energy.

II. PHYSICS OF ATMOSPHERIC MUONS AND NEUTRINOS
Cosmic rays interacting in the atmosphere produce secondary hadrons whose decay products give rise to a spectrum of atmospheric muons and neutrinos.The primary cosmic-ray energy spectrum is approximately a broken power-law and its nuclear mass composition varies as a function of energy.As a special form of the Boltzmann transport equations, the coupled cascade equations describe the evolution of particle fluxes in a dense or gaseous medium.The average number of interactions of a particle with air nuclei is a function of the slant depth or grammage where ρ is the density of the atmosphere and h 0 the altitude above the surface.At high energies inelastic hadronic interactions dominate and result in secondary particle production.Stable and longer lived particles can again interact and produce sub-cascades similar to the initial one but at reduced energy.Unstable particles decay into other hadrons or leptons or re-interact, depending on their energy and lifetime.The cascade evolution stops as the hadrons fall below the threshold for inelastic interactions, leaving only stable hadrons and atmospheric leptons in the cascade.The number of nucleons and mesons increases up to a certain depth or altitude and then attenuates.The production of muons and neutrinos is proportional to the number of decaying mesons.Lepton production therefore decreases at lower altitudes as the atmospheric density increases.For the same reason, the flux of leptons from decay of pions and kaons increases as zenith angle increases.Low energy muons decay preferentially at large inclinations, but at energies above tens to hundreds of GeV most of them reach the ground.

A. Coupled cascade equations
The equations describe the evolution of the differential flux, defined as the differential of the particle flux φ with respect to energy per unit area, unit solid angle, and time The coupled cascade equations dΦ h (E, X) are the transport equations representing the interplay between the two dominating high energy processes, interactions and decays.The energy losses from ionization and multiple scattering ( dE/dX ≈ 2 MeV cm 2 /g) impact muon and electron neutrino fluxes below tens of GeV for near-vertical or below few TeV for near-horizontal directions.The energy dependence of the interaction lengths and decay lengths are shown in Fig. 2.1.For short-lived particles λ dec << λ int , so the secondaries preserve the energy spectrum of their mother species.At the other extreme hadronic interactions dominate, the flux of mesons is attenuated and becomes a power steeper [50].The transition between these two regimes (cross-over between solid and dashed lines in Fig. 2.1) is called critical energy ε h and is a function of the density or the altitude.For typical conditions the values are approximately 115 GeV for π ± , 850 GeV for K ± , 10 PeV for D-mesons and > 10 13 GeV for unflavored η mesons.
2.1.Decay lengths for relevant mesons (solid) compared to the interaction length for kaons calculated at 8 km altitude a.s.l.
Traditional semi-analytical solutions (see for in depth discussion [42,50,68]) have the form The energy dependence of the lepton flux follows the cosmic ray nucleon flux and becomes a power steeper above the critical energy, re-scaled by the cosine of the zenith angle θ.The Z factors are (primary cosmic ray) spectrum weighted moments, representing hadronic interactions and particle decays The values of the spectral index γ are between 2.7 and 3.0 (see Sect.III B).The weight emphasizes the very forward part of the particle spectrum with x Lab 0.2, or in other words, particles with very small scattering angles.The standard approximation is scaling, i.e. the secondary particle spectrum is only a function of x Lab = E secondary /E nucleon and independent of the nucleon's energy.In practice and in current measurements (Fig. 5.5) this approximation is valid in the very forward phase space.But at smaller x Lab and high energies scaling is known to be violated due to multiple parton interactions in one collision.
For very short lived hadrons, such as charmed or unflavored mesons the second term in the denominator of Eq. ( 6) is negligible.This particular fraction of the lepton flux thus follows the spectral index of the primary spectrum up to very high energies.Due to the attenuationless immediate decay of these mesons, the resulting lepton flux is called prompt.The number of prompt muons and neutrinos is small because the mother mesons are rarely produced (small Z N,h,γ ) and because the leptonic decay branching ratios are small (Z h→ ,γ ).All other inclusive leptons (mainly from decay of pions and kaons) constitute the conventional component, which is suppressed at very high energies > 100 TeV so that the small prompt flux eventually dominates the total rate of atmospheric muons and neutrinos.

III. CALCULATION METHOD
We carry out the computation of inclusive atmospheric lepton fluxes with a numerical method described in detail in [41,[43][44][45][46].This method is more powerful and precise than the semi-analytical solutions, especially if the cosmic-ray flux is not a single power law or when scaling violations have to be taken into account.Monte-Carlo calculations [42,60] achieve comparable accuracy but they become inefficient at energies above several hundreds GeV.Biasing techniques can reduce this inefficiency that is related to the absorption of mesons in the atmosphere, however, no such technique is available at present to bias hadronic interaction models for the generation of mesons carrying a large fraction of the projectiles momentum.

A. Matrix cascade equations
The open-source software Matrix Cascade Equations (MCEq) 1 [43] solves the system of discrete coupled cascade equations The index h represents one of the ∼ 65 particle species and the energy index i runs over an energy grid, subdivided into 8 logarithmically spaced bins per decade of energy across 11 orders of magnitude (1 GeV -100 EeV).The coefficients c (E k )→h(Ei) represent inclusive secondary particle energy spectra in the target laboratory frame, and d (E k )→h(Ei) the corresponding decay spectra dN (dec) l(E l )→h(E) /dE from Eq. ( 3).The solver performs a simultaneous integration of a coupled system of N species × N E−bins ∼ 8000 ordinary differential equations (ODEs) in parallel.For particles which suffer continuous energy losses, the ODEs become partial differential equations.This Fokker-Planck part of equation is solved through a finite differences operator of order-3 or higher.The resulting expression in matrix form is The matrices C and D contain the coefficients c and d arranged in a way so as to represent the coupling sums (source terms) from Eq. ( 8) above.∇ µ is the a first derivative operator, µ the mean energy loss or the stopping power of muons in dry air, with its energy dependence fully accounted for and arranged on the energy grid.The decay and interaction coefficients are obtained from Monte-Carlo simulations of particle interactions with air and of free decays.Interactions are simulated with Sibyll 2.3c and other cosmic ray interaction models, decays with pythia 8 [88,89].The computation is significantly accelerated by using sparse matrices and a method to reduce the stiffness [43].Depending on the choice of models and the zenith angle, the solver needs between 0.1s and a few seconds on a single high-performance x86 core.MCEq is a relatively mature software that is gaining popularity among the neutrino communities and it has been used for several practical applications, e.g.[8,30,39].

B. Primary cosmic ray flux
For the discussion in this paper we choose one representative realistic model (called H3a [51]) that contains the important features of the cosmic ray flux including the knee and the ankle.The origin of these features is understood in terms of transitions between classes of cosmic accelerators, such as supernova remnants, pulsars, gamma-ray bursts or active galactic nuclei.The flux of each mass component φ i is modeled as a broken powerlaw, where each source class j accelerates nuclei to a maximal cut-off rigidity R c Our calculation method uses the superposition approximation, where the flux of interacting nuclei is expressed in terms of the individual nucleons (all-nucleon flux) or more precisely an energy dependent proton and neutron flux that can be obtained by substituting A 2 i with A i Z i or A i N, respectively.Fig. 3.1 summarizes the different properties of this model.Upper panel: Flux of nucleons in (m2 s sr GeV 1.6 ) −1 , converted from the flux of cosmic ray nuclei as predicted by the H3a model.Middle panel: spectral index, log-derivative of the upper curve.Lower panel: Fraction of neutrons in the nucleon flux.The spectral softening at around 1 PeV is an effect of the knee and the hardening at ∼ 100 PeV of the ankle of cosmic rays.

IV. CONNECTION BETWEEN INCLUSIVE LEPTONS AND HADRONIC INTERACTIONS
A detailed study that relates the relevant properties of hadronic interactions to the spectra of atmospheric leptons is [85].Here we revisit this topic using the new Sibyll-2.3cevent generator and the efficient numerical scheme, extending it to higher energies and to prompt leptons.All Figures and computations are made with Sibyll-2.3cif not otherwise noted.

A. Distribution of cosmic ray energies
The energy spectrum of the primary nucleons that take part in the production of leptons at a fixed energy is shown in Fig. 4.1.This probability density function (PDF) distribution can be obtained by calculating the contribution from each primary energy bin to the total inclusive flux.The distributions are replotted as a function of the energy fraction of the primary nucleon in Fig. 4.2.The primary cosmic ray nucleon energies peak at 10 × E lepton with a long tail extending to the highest energies, meaning that there is a non negligible probability that the primary cosmic ray can carry significantly more energy than an observed lepton.Since muons mostly originate from pion decays, they preserve a larger fraction of the momentum of the parent meson on average and therefore have a most probable primary energy is somewhat lower than for neutrinos (around 8E µ ).This is a purely kinematical effect of the large muon/pion mass ratio.
The shape of the distribution is significantly affected by the choice of the primary flux model, in particular by the position of the knee and the spectral indices of the all-nucleon flux.Further, it also depends on the type and the longitudinal spectrum of the mother meson.For this reason the shape of the primary energy distribution is energy dependent as illustrated in Fig. 4.2.
At energies well below a PeV, the muons mostly originate from pions, and the shape is almost universal (blue lines in the upper panel of Fig. 4.2), since dσ π ± /dx Lab scales and the cosmic ray spectral index is constant.At higher energies, the spectral index of the primary flux becomes steeper since the cosmic ray spectrum is probed above the knee.Further, heavy flavor production becomes significant and the prompt flux dominates, explaining the large differences in the orange and red curves.For muon neutrinos (middle panel of Fig. 4.2) the shape changes across the entire energy range.One of the reasons is that at lower energies the dominant mother particles are pions, at intermediate energies kaons and at the highest energies mostly D mesons.The different production and decay properties result in a larger variation for the peak cosmic ray energy between 5E νµ and 20E νµ .We can conclude that for conventional muons the relation to the primary nucleon energy scales, but not if prompt fluxes are taken into account.For muon neutrinos the scaling assumption is not accurate.Electron neutrinos originate mostly from charged and neutral kaon decays, which have similar production and kinematics.The distributions look similar to the muon case, peaking at 10E νµ at intermediate energies.
It is interesting to note that the corresponding centerof-mass (c.m. ) energies, relevant for atmospheric lepton production up to multi-PeV energies, are within reach of current particle colliders.However, it is very challenging to create and operate a detector technology that is close enough to the beam to cover the relevant (as we discuss later) longitudinal momentum range with x Lab or x F > 0.2 2 .Another obstacle is that nuclear interactions at TeV energies are currently probed in lead-lead or protonlead collisions at the LHC, but lighter ions in the mass range of air molecules are accessible only in simulations.

B. Relevant hadrons for inclusive lepton production
We continue discussing the connection between hadronic interactions and inclusive leptons by looking at the different hadron species that give rise to a sizable production of inclusive leptons.Most atmospheric leptons originate from weak and partly from electromagnetic decays of the most abundant mesons, i.e. charged pions and kaons.Two aspects of particle production are relevant here; the longitudinal production spectrum, such as the p T integrated differential cross section, and, the energy distribution among the decay products and their inclusive branching ratios.
The hadronization routines in Sibyll can essentially produce all relevant hadrons and resonances up to masses of the Ω ccc baryon.Inclusive p T -integrated cross sections dσ h /dx Lab are computed for each hadronic species irrespective of its life time 3 .Decays are tabulated separately by using pythia 8 [88,89].The inelastic interaction cross sections of more exotic hadrons are assumed to be equal to σ prod nucleon-Air for all baryons, σ prod π ± −Air for pions and light mesons, and σ prod K ± -Air for heavier mesons including charmed.
The various groups of mother particles that directly decay into leptons and contribute to the inclusive flux are shown in Fig. 4.3.Sub-leading contributions are summed together in the "other" groups.As the energy increases, the decays of particles become suppressed above the critical energy.Heavier and less abundant hadrons dominate at very high energy and produce prompt atmospheric leptons.As expected, the conventional muons stem from the decays of charged pions and, with a smaller contribution, from charged kaons (upper left panel in Fig. 4.3).Prompt muons have two sources, decays of charmed mesons and a component from electromagnetic decays of unflavored mesons [63].The detailed break-down of the contributors to the unflavored component is shown Fig. 4.4.A contribution at a similar level as charm comes from the process η and η → µ + µ − γ, breaking the correlation between prompt fluxes of muons and neutrinos.The crossover between conventional and prompt flux happens at  several PeV and depends on the choice of models and the zenith angle.Further sources of high energy muons that are not included in our calculation are the photoproduction of muon pairs, which is suppressed by 10 −4 wrt. the pair production cross section σ e + e − [75], and the nuclear interactions of muons.While the muon pair production can significantly contribute to inclusive fluxes at very high (PeV) energies, the nuclear interactions are only important for the low energy tail of muon bundles in air showers.

At E
100 GeV the main source of muon neutrinos (upper right panel) are semi-leptonic and 3-body decays of charged kaons, see e.g.[61] for a more detailed discussion of relevant channels.Pion and muon decays dominate below this energy.Prompt neutrinos originate from decays of charged and neutral D-mesons, where the fluxes from D ± are a factor of three higher.Since pions do not decay into electron neutrinos (lower left panel), those come mostly from decays of neutral and charged kaons.At energies below 100 GeV and for near-horizontal zenith angles the dominant fraction of electron neutrinos is from muon decays, resulting in a strong association with the muon flux.In turn, this means that the precision of the electron neutrino prediction for a few to several tens of GeV is linked to the modeling of pion production and muon energy loss and, to a lesser extent, to kaon production.
Atmospheric tau neutrinos (lower right panel) are rare [35], but we can discuss their flux for completeness.The dominant production channel of tau neutrinos is the decay of D + s → τ + + ν τ , where the subsequent decay of τ → ν τ + X is more efficient in producing a forward tau neutrino, than the decay of the meson.Therefore most of the tau neutrino flux comes from the decay of the tau lepton itself (black and blue line in lower right panel in Fig. 4.3).
Other sources of atmospheric leptons that are not taken into account in our calculation are B-hadrons.Their contribution to the prompt flux can be of the order of 10% [64,72].

C. Muon charge ratio
The inclusive muon charge ratio R µ has been perceived in the literature as the observable with one of the highest sensitivity to the hadronic physics in atmospheric cascades [14,18,51].R µ is approximately 1.25 below 1 TeV and increases to slightly above 1.35 at higher energies.The reason for this transition can be seen in Fig. Red dashed lines show the contribution to the muon spectrum by interactions of unstable particles (π ± , Λ, etc.) and solid blue the fraction of muons that is related to hadronic interactions in the first interaction length (∼ 90 g/cm 2 ) of the cosmic ray nucleon.The upper panel is computed for a single 100 PeV proton and a vertical incidence angle.
charge ratio of the pion component, to a higher energy range where kaons become more important.
Direct constraints on the production of pions or kaons in proton-air interactions from R µ directly are difficult to derive, since the spectral index and the neutron fraction in the cosmic nucleon flux introduce some degeneracies.While the charge ratio R µ (E µ ) at a fixed E µ covers a range of primary interaction energies (compare with Fig. 4.1), the forward particle spectra approximately scale (see Section V A and Fig. 5.6), alleviating this problem in the interpretation in terms of hadronic interactions.More important is the degeneracy between the shape of the particle production spectrum and the spectral index of primary cosmic rays, or even the spectrum of secondary projectiles downstream of the air shower.As we discuss below, a sizable fraction of inclusive muons stems from these secondary interactions and thus changes as a function of depth.
Eq. (7) demonstrates that the size of the pion and the kaon component is controlled by the convolution of the projectile spectrum and the particle production spectrum (see Eq. ( 7)).Since the charge ratio is not constant along x Lab (see Fig. 5.3), R µ can only constrain the weighted integral Eq. (7) and not the the spectrum of leading mesons directly.This implies that in order to access the microscopic hadronic physics through inclusive lepton measurements one necessarily needs to take multiple correlated observables into account, like the angular distributions of muons and R µ together with the inclusive muon neutrino and electron neutrino fluxes.

D. Inclusive muons vs. muon bundles
The term muon bundles refers to the muon signature that is usually observed in volumetric underground detectors, such as L3+c, MINOS or IceCube.High energy air-showers produce large numbers of muons distributed over a wide range of energies from below 1 GeV up to energies close to those of the shower-initiating cosmic ray.Most energy is concentrated in a small cone around the shower core, which is aligned with the cosmic ray direction.The overburden, rock for underground detectors or water/ice for Cherenkov neutrino telescopes, absorbs the low energy part of the muon content and a small fraction containing between one and several hundreds of muons reaches the detector.Thus, a muon bundle contains several muons originating from the same air-shower.
Inclusive muon fluxes are obtained by scoring the energy of each individual muon and integrating over time.The time integration translates into an integration over the cosmic ray spectrum at the top of the atmosphere.For these reasons, observations of muon bundles or muons in air-showers, and inclusive muon fluxes are sensitive to different features of hadronic interactions.Fig. 4.6 illustrates this difference, showing the contribution from (secondary) interactions of unstable particles, such as pions, kaons or Λ baryons.The arguments for muon bundles apply equivalently to air shower observations.The majority of low energy muons, as it would be observed by an air-shower detector, are created in interactions of secondaries further downstream of the shower.The very first high energy interactions during the first generations of the air shower influence the shower development further downstream to a lesser extent.For the inclusive fluxes, however, the integration over the steep primary cosmic-ray spectrum results in a suppression of the importance of these re-interactions of secondaries.More than half of the muons originate from the first generations of particles emerging from the highest energy interactions.

E. Hadrons that do not decay into leptons
In the previous sections, we discussed the role of hadrons with sizable leptonic branching ratios that can directly contribute to the inclusive flux.Here, we outline why an accurate description of the other particles (with rare leptonic decays) is indispensable in inclusive flux calculations, how these hadrons are related to the conventional pion and kaon components.
The unstable particles in atmospheric cascades decay into the lightest unstable species (π ± , K, etc.), which ultimately decay into leptons or nucleons before reaching the ground.The life time of particles (for instance) ρ-mesons is so short that their explicit presence in the transport equations has little impact on the development.Through the decay into two pions (BR ∼ 100% the ρ mesons feed down into the secondary spectrum of pions.Therefore, they have to be taken into account either explicitly (in the transport equations) or implicitly in the inclusive production spectrum of pions.As it has been discussed in [38,77], leading ρ meson in pion-nucleus interactions (instead of leading π 0 that feed the electromagnetic cascade) can significantly affect the development of the muon content in an extensive air shower.Fig. 4.7 shows the feed down from vector mesons, strange baryons and resonances to the inclusive yields of pions, kaons and protons in p-Air interactions.These additional channels are remarkably large, demonstrating the importance of choosing accurately the definition of stable/resolved particles when comparing and tuning event generators to accelerator data or to other models.This example clearly demonstrates why the development of interaction models requires significant effort and why there are no accessible "knobs" for "tuning" a model to, for example, pion measurements.An enhanced production of light mesons would necessarily lead to reduced production of vector mesons or baryons, creating tension with other observables.As described in Section V, Sibyll uses the Lund string fragmentation model [32] that consistently connects the production of light and heavier mesons to the available energy in color strings using a set of phenomenological probability parameters.Simply speaking, in an ideal world these parameters should be universal and can be derived from measurements at lepton colliders.However in practice this is not always true, in particular in hadronic collisions at very high energies where additional non-perturbative methods are invoked to describe hadronization [88], at least when sticking to the fragmentation of color strings as the underlying picture for hadron formation.

F. Angular distribution
The panels of Fig. 4.8 show the dependence of the differential fluxes of leptons on the zenith angle, as it would be perceived by a detector located at the surface.Traditionally, this behavior is described by a combination of the steepening of a spectral components (from µ, π or K decay) that starts above different critical energies (see Eq. ( 6)).The critical energies are approximately 1, 115 and 850 GeV for µ, π and K mothers, respectively [50].These energies represent the transition from a decay to an interaction dominated transport and in reality they depend on the atmospheric model, the production height and the zenith angle.Muon decay is an additional process that affects the angular distribution of muon and neutrino fluxes at low energies, where a deviation from the typical 1/ cos (θ) shape can be seen in the upper left panel, and even at higher energies for very extreme angles.
The strong enhancement of near-horizontal fluxes originates from the suppression of secondary interactions of mesons, since the bulk of high energy particles from the first interactions are traversing a very thin atmosphere  The signature of a "prompt" component is a flat zenith distribution, meaning that these "mother-mesons" do not interact with the atmosphere before decay.At energies above the critical energy this probability increases and depends on the density gradient along the cascade trajectory.The gradient is smaller for horizontal trajectories, thus the probability for unstable particles to decay becomes higher again.The prompt component from charmed and unflavored meson decay is flat up to the highest energies due to the short life times these hadrons.PDF of the relevant x Lab values for inclusive leptons at fixed energy.These curves are derived from the inclusive flux, i.e. integrated over the cosmic ray spectrum.An energy independent Z-factor approach would results in a universal shape.Here, the energy dependence comes from deviations from perfect Feynman scaling and from changes in the spectral index of the primary nucleon flux due to the knee and ankle.and do not accumulate enough depth.The distinctive feature of the prompt component is the flat angular distribution.For fluxes from charmed meson decays, reinteractions with the atmosphere become relevant above 5 − 10 PeV (not shown in Fig. 4.8) leading to deviations from the flat distributions.
At very high energies the flux is not always dominated by prompt leptons.As the right panels clearly show, the prompt component is dominant for vertical directions but becomes sub-dominant towards the horizon.For the detection of prompt fluxes free from the contamination with the conventional component, it would be necessary to search in electron neutrinos at several hundreds of TeV, or, in muon neutrinos above PeV energies, preferably in vertical directions.With present neutrino telescopes, such as IceCube, ANTARES and GVD, the detection is extremely challenging since the prompt flux is low, and, diluted by backgrounds from muon bundles in the potentially interesting channels (vertical down-going or electron neutrino cascades).

G. Particle production phase space
In section IV A, we outlined that atmospheric leptons probe hadronic interactions at center-of-mass energies that are accessible to recent accelerators and we showed that the most probable energy of the incident cosmic ray is E CR ∼ 10E µ or very approximately at where the 0.6 is the average energy fraction transferred to muons in pion decays.The cascade evolution implies additional factors such as secondary interactions, contributions from other decay channels the decrease of the nucleon energy in subsequent interactions etc.A PDF for x Lab values that contribute to the flux of leptons at a fixed energy is shown in Fig. 4.9.The energy dependence of these curves can be explained analogously to the energy behavior of Fig. 4.2 and comes from deviations from ideal Feynman scaling of the interaction model and the knee and ankle in the cosmic rays.Note that the x Lab values commonly refer to all hadrons including (p, n, π, K).A value of x Lab = 0.3 for an incident proton at 10 PeV, for example, implies that the scattering angle of the secondary particle is of the order of few µrad and impossible to detect at a particle collider experiment without dedicated detectors.At high energies, the hadronic interactions are only accessible through simulation, and thus constitute the dominant source of uncertainty in the calculations of inclusive fluxes.

V. HADRON PRODUCTION IN SIBYLL 2.3C
The hadron interaction model Sibyll is designed mainly for the use in cosmic ray air shower simulations.While the general features of QCD like quark confinement, multiple interactions and jet production are included in the model, particular features that are relevant for the development of air showers, like diffraction dissociation and forward particle flow are implemented in more detail.
The interaction model in Sibyll is based on the twocomponent dual parton model with soft and hard minijets [36].It also includes low-and high mass diffraction and a model for the excitation of beam remnants [21,83].The hadronization model is based on the Lund string fragmentation model [32,87].Hard scattering is distinguished from soft scattering by a cutoff in transverse momentum.The cross section for hard scattering is calculated to LO in QCD at the scale defined by the cutoff in p T .Saturation is included by means of an energy dependent increase of the p T scale.Contributions from quarks of all flavors and gluons are included with their full kinematics as determined by the parton distribution functions.However, in the subsequent fragmentation of the scattered partons, no distinction between the different flavors is made.The string (color-flow) configuration of all the parton interactions are treated as gluon gluon scattering (see Fig. 5.7 for illustration).The soft interaction cross section is modeled with a parameterization based on the Regge field theory [37].
Two aspects of hadron interactions are improved in the latest version of Sibyll motivated by the discussion above.These are the production of leading particles and the production of charmed hadrons.The new treatment of leading particles of the remnant model is discussed in Sect.V A) and the charm model is covered in Sect.V B.

A. Leading particles
While the particle production in forward phase space (x Lab > 0.2) plays an important role for the fluxes of leptons in the atmosphere (see Sect.IV G), it is even more important for the charge and flavor ratios (see Sect.IV C).The reason is that the weight applied on the longitudinal spectrum emphasizes the forward phase space (Eq.( 7)), where the flavor asymmetry is strongest.This forward asymmetry is also known as leading particle effect and it is related to the valence quarks dominating the momentum distributions of hadrons in soft interactions.
In Sibyll the valence quarks are assumed to undergo a soft interaction (two-string model), where each hadron is split into two valence partons (baryon: quark and di-quark, meson: quark and antiquark) and two color strings are formed between the partons of the two hadrons.The basic version of this non-perturbative configuration corresponds to the upper baryon in Fig. 5.8.The nature of the leading baryon is determined by the flavor of the quark (from the q q pair) with which it recombines.(The lower baryon in Fig. 5.8 illustrates remnant splitting channel to be discussed in the next paragraph.)In the standard fragmentation of a projectile baryon, the fraction of the momentum carried by the valence quark is assumed to follow and the di-quark (quark or antiquark for mesons) is assigned the remaining momentum x diq = 1 − x q .In Sibyll-2.1 leading particles emerge from the strings due to the large momentum fractions of the di-quark in combination with a hard fragmentation function used only for the hadrons produced at the string ends that are related to the beam particle.
In Sibyll-2.3cthe hard fragmentation is replaced by a beam remnant model [38,69,90], where the valence quarks are separated from the sea partons that fragment as an independent system.The configuration of strings in the presence of a remnant is illustrated in the lower baryon line of Fig. 5.8.The fraction of momentum assigned to the remnant is taken from f r (x) = x 1.5 and the excitation mass is distributed ∝ M −3 r where the lower limit is the hadron mass and the maximal mass is M 2 r s −1 = 0.02.The energy required for the excitation is transferred from the other hadron.
To account for the possible absorption of the valence quarks in the scattering process the remnant is formed with the constant probability of 60%, suppressed by the number of soft and hard parton interactions and, in case of nuclear interactions, with the number of nucleons involved P r = P Nw+0.2 (n soft +n hard ) r,0 , with P r,0 = 0.6.Through this mechanism the leading proton distribution obtains the characteristic dip in the transition from hadron to nuclear targets, while the meson spectra are not affected.At the same time the charge ratio of leading pions can be described more accurately as the fragmentation of the proton remnant through nucleon resonances.The small strings preserve the isospin of the proton and favor the production of π + (see Fig. 5.1).For kaons the charge ratio is affected even more strongly (Fig. 5.2) as the nucleon can only transition into a hyperon and a positive charged kaon, e.g.
This model yields a viable explanation for the effect of associated production, in which both Λ and K + exhibit a leading particle effect.The models, in which the strings span between the valence quarks and ss pairs from the sea do not explain such a strong forward enhancement, since kinematically the strange hadrons form centrally.In Fig. 5.3 the resulting charge ratios for pions and kaons are compared with measurements in NA49 [24,25,29].The new model describes the pion charge ratio well.For kaons the description has improved but towards large x F the ratio is still overestimated.
Unfortunately there are no data on leading meson production and charge ratios available at √ s > 50 GeV to directly test the model.Indirectly the model is constrained by the spectrum of leading neutral pions measured at the very forward spectrometer LHCf [17].While the model reproduces the spectral shape in x F , the transverse momentum distribution is not described as well (Fig. 5.4).
LHCb [1], the detector that provides particle identification with the largest rapidity coverage at the LHC, only covers a small region of longitudinal phase space x F 0.1 and therefore can not constrain the leading charge ratios.In a proposed fixed target configuration for LHCb [67,91], the charge ratios in the forward region could be determined for √ s = 104 GeV.By comparing with Fig. 4.2, this would correspond to the charge and flavor ratios of muons and neutrinos in the TeV to PeV range.One of the central assumptions in analytical calculations of the atmospheric fluxes is the scaling behavior (energy independence) of the particle production spectra at large Feynman-x with energy, the so-called Feynman scaling [47].According to measurements the scaling hypothesis holds at least up to 7 TeV [17,79].Fig. 5.5 shows the scaling behavior of Sibyll in comparison with the data, which is compatible within the errors of the data.
In the calculation of the atmospheric fluxes the production spectra enter in a weighted form (e.g.Eq. ( 7)), emphasizing the very forward region.At the same time the primary flux covers center-of-mass energies from 10 GeV to 100 TeV.The production spectra for Sibyll with a weight for a typical cosmic ray spectrum are shown in Fig. 5.6.From the presence of scaling violations in the central region due to multi-parton interactions a slight softening of the spectra with energy is expected [78].
With the new remnant model Sibyll initially showed a hardening of the meson spectra with interaction energy 4 .This behavior was found to originate from the splitting of leading di-quarks, originally introduced to improve the leading flavor ratios [82].With Sibyll-2.3cscaling is approximately fulfilled.
The effects of these changes on the atmospheric flux heavy quarks due to the large mass (m cc ≈ O(2 GeV)) is well above the soft scale.Production of charmed hadrons can therefore be expected to be dominated by hard (perturbative) processes as in Fig. 5.7.On the other hand, measurements of leading charm production (x F > 0.6) at low energy [23,53] indicate that there is a soft (nonperturbative) component [48,49].Associated production of charm, (e.g.production of Λ c + D 0 ) is illustrated in Fig. 5.8.The model of charm quark production in Sibyll-2.3cuses the family connection between strange and charmed hadrons, exchanging s with c quarks in the fragmentation step [22].The total rate of charm production in this model is set by the probability P s→c for the replacement of a strange quark by a charm quark.At energies beyond a TeV, when the mass of the charm quark becomes negligible in comparison with the scale of the parton interactions, the difference between the light flavors and charm vanishes and the energy dependence is given by the minijet cross section.In this energy regime a fixed charm to strange rate P s→c is appropriate.In the threshold region charm production increases much more rapidly with the c.m. energy.To account for this threshold the charm rate in minijets decreases as where √ ŝ is the c.m. energy of the frame of the partons and m eff = 10 GeV is the effective charm mass.The fragmentation function used for charm quarks is the SLAC/Peterson parameterization with = 2.0 [80].The transverse momentum of the charm pair in the string is taken from an exponential transverse mass distribution with energy dependent mean m T (s) = p T,0 + (0.3 GeV) log 10 ( √ s/ 30 GeV) , where p T,0 is 0.3 GeV for charmed mesons and 0.5 GeV for charmed baryons.The comparison of the inclusive cc production cross section of the model with measurements in a wide range of energies is presented in Fig. 5.9.Apart from the full inclusive cross section, the figure also shows the cross section in central (ALICE), next-to-central (LHCb) and forward phase space.The model correctly describes the growth with energy in the most central phase space, but not as well in the next-to-central region.The different behavior in the threshold region and at large energies can be seen by comparing charm production in the model with the scaled minijet cross section σ QCD in Fig. 5.9.
While this phenomenological model of charm production is sufficient to describe the total charm yield, it has difficulties to describe differential spectra, like the transverse momentum spectra of charmed hadrons measured by LHCb [2,3] and ALICE [11], as it neglects the details of hard scattering.To account for the dominant contribution from hard scattering without altering the minijet model at the basis of Sibyll, the production of charm quarks is focused towards the string ends by invoking a separate P Minijet s→c for the gluon splitting (g → qq), as illustrated in Fig. 5.7.The nature of the production of charmed hadrons from soft interactions is not well understood.The spectra of leading D-mesons in p-, π-, and K-nucleus interactions measured at the E769 experiment, together with the asymmetry between Λ c and its antiparticle observed in pp scattering by the SELEX collaboration, indicate the presence of a non-perturbative, soft component [27,53].Irrespective of the exact origin of the leading charm production, the hard minijet component, representing perturbative interactions in Sibyll is not sufficient to describe the low energy data.In Sibyll soft interactions include soft minijets, the scattering of the valence quarks, the formation of beam remnants and diffraction dissoci- ation.The soft minijets in the model derive from soft gluon scattering which has a steep longitudinal momentum distribution, resulting again in mostly central production.In contrast, the momentum spectrum for valence quarks and the remnant is much harder, so the hadron spectra at large x F are determined by these processes.Fig. 5.10 shows how the different processes form the D-meson spectrum in pion-nucleus interactions at for the leading charm production in the soft processes are sketched in Fig. 5.8.
The production of D-mesons in the central region |y| < 0.5 as a function of the transverse momentum is shown in Fig. 5.11, compared with the measurement from ALICE at 7 TeV c.m. energy [11].Since the central region is dominated by hard scattering, the charm rate for the hard minijets (P Minijet s→c ) is determined from these data.The difficulty in describing the growth rate of the total yield in the LHCb phase space in Fig. 5.9 also appears in the differential spectra shown in Fig. 5. 13.The yield at 7 TeV is overestimated, and the shape of the p T spectra deviates at higher p T values and at large rapidities.The same trend continues at 13 TeV, although the model predicts the correct overall yield.These discrepancies are probably connected to the approximation of the perturbative charm production with the hard scattering cross section and the limitations of the simple minijet  model.The parameters of soft charm production are determined by the low energy pion-nucleon and proton-proton data shown in Fig. 5.12 and Fig. 5.10.It is interesting to note that in order to describe the forward production at E769, charm production in the central strings and in soft gluons is sufficient.No charm production in the hadronization of the remnant is required.The numerical values of the parameters of the charm model are summarized in Tab.I.An overview of the available measurements that have been used for the determination of free model parameters is in Table II.

Comparison with other models & discussion
Fig. 5.14 compares the transverse momentum spectrum of neutral D-mesons in Sibyll and a more fundamental next-to-leading order (NLO) QCD calculation based on POWHEG in the LHCb phase space [55].The spectra match nicely within the theory uncertainty, indicating that our phenomenological charm model is suf-ficient to describe the transverse momenta for perturbatively produced charm.As mentioned previously, the hard scale in the model is not equivalent to the scale in pQCD calculations.However, in the (still) central phase space covered by LHCb, the comparison is valid between pQCD calculations and the full prediction from Sibyll.
For the inclusive lepton fluxes, the particle production in longitudinal phase space is paramount.However, calculations that extend into the region of large Feynman-x / rapidity are difficult as they include scatterings smallx (y ∼ |x 1 − x 2 |, x 1 , x 2 are the momentum fractions of the partons) and typically require additional assumptions.Fig. 5.15 shows the x F spectrum of D-mesons at 13 TeV, comparing Sibyll and a small-x color-dipole calculation including saturation effects (GBW) [57,58,71].In the figure the GBW calculation is re-scaled to match our spectrum at x F ≈ 0.4.Similar to the POWHEG calculation (Fig. 5.14), the pert.* component in Sibyll and the GBW calculation match well, in particular when the SLAC/Peterson fragmentation function is applied.In our model, the additional non-perturbative components (valence quark interactions, remnant and diffractive processes) start to dominate the spectrum spectrum around x F = 0.7.Note that a weight of x 2 F is applied in the figure, which is the appropriate weight to study prompt leptons which originate from interaction energies above the cosmic ray knee.Fig. 5. 16 shows a comparison to more recent NLO QCD calculations, PROSA 2017 [54] and GMVF [33].The PROSA calculation takes into account charm production from hard scattering with leading logarithmic corrections from parton showers.Non-perturbative effects from hadronization and leading particle effects are included through PYTHIA 8 [88].In the GMVF calculation hadronization is included via fragmentation functions.In contrast to GBW, PROSA and GMVF do not include a specific model for the small-x region of the parton distribution functions, instead they are fixed by and extrapolated from HERA and LHCb data [92].The transition from proton-proton to proton-air interactions is simplified through a superposition model σ cc,p−Air = 14.5 × σ cc,pp .Within the theory uncertainty (see bands in the figure) which contains contributions from the variation of the renormalization and factorization scales, the x Lab -spectra from PROSA and Sibyll are compatible.While the difference in the normalization between GMVF and Sibyll is larger than for PROSA, the shape of the spectrum is more similar.Note that due to the additional scattering centers in the target nucleus, it is expected that the production spectra are attenuated at large x Lab .Since the superposition ansatz is used for PROSA and GMVF, the difference with Sibyll increases for protonair.Finally, Fig. 5.17 shows the spectrum of D-mesons as a function of the energy fraction in the lab.frame, also weighted by x 2  Lab , and broken down into different production processes.Their relative contribution to the spectrum weighted moment (Z-factor) is printed as a percentage value.While the soft processes in Sibyll (valence, diffractive and remnant) are important to describe the shape of the longitudinal spectra at low energy, their contribution to the Z-factor at high energy is only of the order of 10 %.Minijet charm production, in contrast, sums up to 89 %, resulting in a seemingly well determined the Z-factor.Despite the fact that the model only approximates perturbative charm production and shows slight deviations in the next-to-central phase space, the constraint on the Z-factor will be not quite as strong.Quantitatively, while 89 % of charm is produced in perturbative processes in the model, just 3 % are covered and constrained by collider measurements (LHCb).Under the assumption that the phase space extrapolation of the perturbative processes is well understood, this leaves about 10 % of charm production (valence, diffractive and remnant processes) unconstrained by the LHC.
With current LHC detectors there is little hope to directly constrain the model at higher energies, since the LHCb detector is already hitting the technological limit pT-spectrum of neutral D-mesons at √ s = 13 TeV.Compared are Sibyll and a NLO QCD calculation (POWHEG) [55].The band of the theory calculation corresponds to charm quark mass and factorization scale uncertainty.The central phase space covered by the LHCb acceptance is dominated by perturbative processes, such that the full prediction by Sibyll can be compared to the pQCD calculation.
concerning particle identification at small angles.Theoretically the model may be constrained by improving the perturbative model, in order to reduce the theoretical uncertainties.Additional input is expected from (next generation) large volume neutrino telescopes that can be sensitive to the prompt neutrino flux through a measurement of certain ratios (see Sects.VI B 3 and VI B 4).
At lower energies on the other hand the fixed-target configuration of LHCb [67,91] may provide valuable input.While the boost from c.m. to the lab.frame worsens the situation on the beam side, the target region will be shifted into the acceptance of LHCb.Measuring the production of D-mesons in fixed-target proton- Sibyll is compared with a calculation in perturbative QCD using extensions to small-x (GBW) [57,58,71].GBW model is scaled to match Sibyll at xF ∼ 0.4.While the yields are different, the shape of the spectrum predicted for the perturbative component is comparable.Comparison of the energy spectra of mesons between Sibyll (red lines) and NLO QCD calculations (PROSA 2017 [54] (black lines) and GMVF [33] (purple lines)).Energy is evaluated in the lab.frame.The lower panel shows the difference relative to Sibyll.The grey band represents the theoretical uncertainty in the PROSA calculation (factorization and renormalization scale).Calculations for pp interactions are shown by solid lines and those for p-air with are dashed or dotted.The NLO calculations for an Air target are scaled up according to 14.5×pp.
proton, or ideally oxygen-proton, collisions with LHC beams (6.5 TeV) could determine the x L -spectrum in the intermediate energy range of √ s ∼ 110 GeV.For comparison, the current highest energy measurement of the longitudinal spectrum of D-mesons is at 40 GeV.

VI. IMPACT OF THE IMPROVED SIBYLL-2.3C ON INCLUSIVE FLUX CALCULATIONS
In this last section, we challenge the new model against its predecessor Sibyll-2.1 and computations using other recent hadronic interactions models.The model EPOS-LHC [81] is currently very successful in describing cosmic ray air-shower observations and minimum-bias collider data at the same time.A representative model using the Quark-Gluon-String framework (an analogous approach to the Dual Parton Model) is QGSJet-II-04 [76].Both models have been revisited after the launch of the LHC and have been adjusted to similarly recent data as Sibyll-2.3c.We can easily swap the interaction model in our MCEq calculations and keep all other parameters equal to the Sibyll-2.3ccalculations when doing this.
In addition, we compare to the very detailed predictions of atmospheric neutrinos from the two state-of-theart calculations; HKKMS 2015 [60] is based on a fully 3-dimensional geometry (3D) at energies below 30 GeV and on a 1D Monte-Carlo calculation at higher energies.The other reference is known as the Bartol calculation [31] and consists of a 3D part at energies below 10 GeV and on a 1D Monte Carlo above that.The history of the HKKM and the Bartol calculations started over 20 years ago and played an important role in the era of the Super-Kamiokande detector and the discovery of neutrino oscillations.
Both calculations use similar technology and the parameterization of the cosmic ray flux from [52].In HKKMS, hadronic interactions are modeled with an effective particle Monte Carlo that is based on inclusive parameterizations of secondary particle spectra from the DPMJET-III event generator [84].Since DPMJET fails to describe atmospheric muon fluxes to the required precision, some corrections derived from muon measurements have been applied [62].These corrections do not have a microscopic explanation and depend on the choice of the primary cosmic ray parameterization.
The Bartol calculation uses the Target-2.1aMonte Carlo particle generator, an effective method that generates secondary particles according to parameterizations of fixed target data, conserving some physical quantities like energy and momentum.Its advantage is that in the energy range covered by fixed target data, the resulting spectra will converge to these measurements for a large number of Monte Carlo samples.However, the lack of a detailed microphysical model limits the extrapolation capabilities into a phase space without measurements.An important example of such an extrapolation problem is the associated production of K + in the process p + air → Λ + K + .The authors of [31] seem to have made an optimistic choice for this particular process that is somewhat in tension with the current (microscopic) hadronic interaction models.Since there are no fixed target results on K ± production at projectile momenta above 400 GeV, this choice can not be constrained by data and impacts the uncertainties of the present calculations.
In addition, we compare to other 1D muon [66] and neutrino flux [86] calculations that rely on numerical methods [73,74].There, hadronic interactions are parameterized inclusively using tabulated secondary spectra.This procedure allows to either use effective models that directly parameterize the measurements or secondary spectra from event generators, very similar to the possibilities in MCEq.
In the following sub-sections we explicitly do not aim to discuss the uncertainties of the calculations and leave this topic to a follow-up publication.For this reason we avoid making comparisons with inclusive flux measurements, in particular because the fluxes notably depend on the choice of the cosmic ray spectrum.The reference models were extensively compared to data and can act as a reference point.The H3a nucleon flux (Sect.III B is used as the primary cosmic ray flux throughout all calculations.It is important to keep in mind that the MCEq based calculations are 1D calculations without accounting for the geomagnetic cutoff, solar modulation effects or Earth propagation for up-going fluxes.These approximations result in up-/down-and azimuth symmetric fluxes by design and do not depend on the choice of a "detector" location.As discussed in the literature concerning the 3D modeling, these approximations are good at energies E lepton > 5 GeV.The characterization of fluxes below 5 GeV has been a very active topic in the past (see for example the review [52] or the references in [31,60]) and goes beyond the scope of this work.A. Fluxes

Muons
The break-down of the different hadron species decaying into atmospheric muons and neutrinos is shown in Fig. 6.1 for different choices of the hadronic interaction model.In the energy range where conventional leptons dominate, the recent models Sibyll-2.3c,EPOS-LHC and QGSJet have almost identical behavior.In Sibyll-2.1,however, the kaon component is more abundant.
In Fig. 6.2 we compare muon flux predictions using the different hadronic models in MCEq and another numerical calculation.The spread is within 10% for the post-LHC interaction models and increases to 20% when including Sibyll-2.1.As the muon charge ratio in the right panel suggest, this has to do with an enhanced production of K + that originates from a program artifact in the old version.These K + are copiously overproduced when diffraction occurs on nuclear targets.A correction of this behavior leads to a flat and rather small charge ratio.As expected from the explanations in Section IV C, Sibyll-2.3cand EPOS-LHC predict an increase of the charge ratio at higher energies.The "wavy" behavior of the Kochanov et al. calculation is related to its primary spectrum and not the hadronic model, that assumes scaling at high energies.The MCEq fluxes for different interaction models are calculated for the zenith angle θ = 0 • and the H3a primary spectrum.For the Kochanov et al. [66] curves the authors used the KM parameterization for hadronic interactions and the Zatsepin-Sokolskaya primary spectrum.
In summary, the muon fluxes seem well constrained since they mostly depend on the modeling of the secondary pion production, where no such associated production channel as for K + exists.The charge ratio is more sensitive to the details of forward kaons, as discussed in [51] (see also [18]).While it seems constrained within 15% by the model predictions, this range exceeds typical experimental errors of a few percent.Also, the cosmic ray composition impacts the charge ratio with a similar magnitude.

Neutrinos
The muon neutrino fluxes (left panel of Fig. 6.3) calculated with the recent interaction models are very similar.As mentioned above, the K + issue in Sibyll-2.1 is responsible for the deviation of ∼ 40% from the other models and is indeed unphysical.The disagreement of the neutrino spectral index between our calculations and HKKMS or Bartol is caused in part by the different choice of the primary flux, but it is also impacted by the way kaon production and decay are treated.As a reminder, in MCEq decays are simulated with the pythia-8 Monte Carlo [88] using methods that preserve high accuracy throughout all steps of the calculations.We find, for example, that ours and the calculation by Sinegovskaya et al. [86] (using exactly the same set of models) agree within a few percent for muon neutrinos, but in electron neutrinos do not match (right panel of Fig. 6.3).For Bartol the higher number of electron neutrinos may be the result of a higher kaon abundance in associated K + production in Target-2.1a.In the case of HKKMS it could be connected to the "tuning" to the muon charge ratio.However, we can not verify these conjectures at the level of hadronic interactions.In the case of the angular distribution of inclusive leptons, it may be that differences in geometry and computational efficiency enter.

Prompt fluxes from atmospheric charm
In Fig. 6.4 the prompt lepton flux from Sibyll-2.3c is compared to some of the recent prompt flux calculations by [34,54,56].An advantage of the charm model implemented in a Monte Carlo event generator is the integration into air-shower codes like CORSIKA [59], where it can be used to generate exclusive air-shower or muon bundle events containing the decay products of charmed hadrons and of unflavored mesons.As discussed in Sect.V B 2, the phenomenological approach to heavy flavor production in Sibyll-2.3yields charm cross sections comparable to other contemporary perturbative QCD calculations.The other inputs of the prompt flux computation, such as the proton-air cross section or the elasticity, have some influence, as well.All of the available models are compatible with the LHC measurements and the IceCube bound from [7].As expected from Figs. 5.17 and 5.16 the "perturbative" production from hard processes accounts for the largest fraction of the flux.The remaining contributions from diffraction, remnant excitation and fragmentation account only for a very small part and affect almost exclusively the description of the low energy data, as explained in Sect.V B.
The prompt electron neutrinos dominate over the conventional at much lower energies.For up-going neutrinos Leptons from decay of B mesons do not constitute more than 10% of the prompt flux [71].
The detection of the prompt flux with instruments like IceCube through an excess of the flux is extremely challenging due to the astrophysical neutrino "background" that overshoots the prompt flux by almost an order of magnitude.As discussed later in Sects.VI B 3 and VI B 4, the more sensitive observables are the flavor and the angular ratios, which would be affected by the excess of the prompt over the conventional electron neutrinos drawn as hatched area in Fig. 6.4.There are caveats, though.The spectral index of astrophysical neutrinos is under investigation and it is currently compatible with values between −2 and −3 [7].The limit by IceCube can be thought of as single power-law extrapolation of the current best fit of the through-going astrophysical muon neutrino flux at higher energies [7].In case future studies confirm either a broken power-law (soft at lower, hard at high energies), or, a generally soft astrophysical flux, the hatched area will shrink, making the prompt flux impossible to measure with neutrinos, at least in the standard (1:1:1) flavor scenario.The positive side effect of such a scenario is that the prompt flux would become a negligible background for neutrino astronomy.An alternative method to measure the prompt flux might involve atmospheric muons.As discussed in Sect.IV B, the prompt muon flux is partly independent of the prompt neutrino flux due to the extra component from unflavored meson decays and from electromagnetic pair production.

B. Angular dependence
The dependence of fluxes and flavor ratios on zenith directions is a crucial observable for the measurement of fundamental neutrino properties.In atmospheric neutrino oscillation studies with volumetric detectors, the zenith angle defines the baseline distance over which neutrinos can oscillate.This technique has been applied to obtain evidence for neutrino oscillations with the Super-Kamiokande experiment [65].Recent experiments, such as IceCube/DeepCore, evolved this method and make use of larger detector volumes and different energy bands to determine oscillation parameters to a precision competing with dedicated accelerator setups [9].At much higher energies, the pattern in the zenith-energy plane gives access to physics beyond standard model (BSM) scenarios [8].Future experiments [4,16] will increasingly rely on accurate predictions of angular distributions and require access to the underlying uncertainties of the physical models.In the following sections we scrutinize the calculations obtained from the combination MCEq + Sibyll-2.3cagainst the available reference models.[34], GRRST [56], PROSA [54] and GM-VNFS [33] are computed at NLO precision under different assumptions.They are accompanied by large error bands which are omitted for clarity.Within the errors all models are compatible with each other.The bound by IceCube [7] of 1.06x is expressed in units of normalization of the model from [40], re-computed for the H3a primary flux.The red hatched patch represents the excess of the prompt flux that leads to distinct signatures in the νµ/νe and the vertical-to-horizontal ratio.

Neutrino fluxes
The zenith angle distributions (see Fig. 6.5) of fluxes agree between the calculations within a few to ten percent with some exceptions.Excluding Sibyll-2.1, the dependence on the hadronic interaction model is weak.Muon neutrino fluxes are particularly sensitive to the K + abundance at high energies, since muon neutrinos originate primarily from decay of charged kaons in the TeV-PeV range.
The differences between the Monte Carlo calculations (HKKMS and Bartol) and MCEq do not seem to follow a conclusive pattern pointing to several contributing factors.For muon neutrinos deviations are generally small and at 5 GeV, some of the effects neglected in MCEq might play a role, for instance the geomagnetic cutoff, the onset of 3D effects or the approximation of unpo-larized muon decay might explain the 5% deviation.In electron neutrinos there is some difference at the horizon (note the increased y-axis scale in the lowest right panels).Technically, electron neutrino spectra are harder to compute with Monte Carlo methods, since their abundance in the cascade is a factor ∼ 10 lower.However, the results should in principle behave similarly, since at the energies that show the largest deviations both reference calculations use 1D schemes, in which the methods applied to improve statistics in 3D (such as the virtual detector [31,62]) do not apply.The geometry in MCEq is fully spherical (or polar in case of azimuth symmetry) and does not impose any technical difficulties even for the most extreme horizontal cascades.A possible explanation for the steeper zenith distribution in HKKMS could be a higher abundance of neutral kaons that only affects the electron but not the muon neutrinos and is therefore difficult to spot in other distributions.In HKKMS, the corrections to the secondary spectra of neutral kaons were related to that of charged kaons (to which the muon charge ratio is sensitive) via approximate isospin relations, while in Sibyll and the other interaction models, neutral kaons are a result of the partonic configurations during the interaction and the subsequent hadronization step.

Muon fluxes
The impact of the hadronic model on the angular distribution of muon fluxes is demonstrated in Fig. 6.6.At lower energies, where the muon flux is dominated almost exclusively by the pion component, the models behave similarly.However, at very high energies the angular distributions do not match.This is related to the prompt muon flux from unflavored and charmed meson decays (compare with Fig. 4.3).In addition, some contribution from muon pair production can be expected that is not accounted for in the present calculation.We checked that the angular distributions at high energy exactly match for the conventional fluxes.The previous attempts by the IceCube Collaboration to measure the atmospheric muon fluxes were negatively impacted by a mismatch in the zenith distribution [6,10].While a part of the problems may originate from experimental uncertainties, it would be worth to revisit these measurements with the post-LHC interaction models.As outlined below and in Sect.VI A 3 for neutrinos, the presence of a prompt component affects in particular the angular distributions making them flatter than the conventional-only scenario.

Neutrino ratios
The Fig. 6.7 gives a detailed overview on the behavior of the neutrino ratios.Ratios trace particularly well hadronic interactions, since the dependence on cosmic ray flux mostly cancels out.As illustrated before in  The muon neutrino/anti-neutrino ratio at higher energies in the upper panels separates the calculations, or more precisely the hadronic models of these calculations, into two classes.Those with a very enhanced forward K + production and those without a particular emphasis on this channel.The Bartol and the Sibyll-2.1 curves show a similar trend related to the high abundance of K + , while the EPOS-LHC model has contrary trend albeit less significant.The flat angular dependence of the Bartol curve at 5 GeV (most upper right) seems to be related to the high K + abundance, but without access to the individual hadronic components of this calculation it is hard to to disentangle the exact origin.
The muon-neutrino/anti-neutrino ratio inferred from the energy-dependence of the muon charge ratio in [51] is closer to the Sibyll-2.1 value, rising from 1.5 at 10 GeV to 2.2 in the TeV range and above.This is a consequence of the fact that, because of the two-body decay kinematics, most muon neutrinos come from decays of charged kaons rather than pions above 100 GeV.Fitting the increase in the measured muon charge ratio, after accounting for the neutron fraction in the primary cosmic-ray beam, normalizes the kaon contribution to the muon flux.The resulting muon charge ratio obtained in this way increases from 1.28 at 10 GeV to 1.41 at 10 TeV, some- what higher than the corresponding ratio from Sibyll-2.3c in Fig. 6.2.The corresponding muon-neutrino/antineutrino ratio is amplified by the kinematic effect, so that its difference from Sibyll-2.3c in Fig. 6.7 is greater.The electron neutrino/anti-neutrino ratio (middle panels of Fig. 6.7) processes similar discrepancies among the models.The abundance of K + dictates the ratio at energies above a hundred GeV.The up-turn of the HKKMS calculation above a TeV might be a relic of the corrections applied to fit the muon charge ratio, and this behavior can indeed be more realistic than the flatter distribution predicted by the interaction models in MCEq.During the development of Sibyll-2.3c,we aimed to have an accurate microscopical description for the muon charge ratio, but despite a significant improvement the result is not perfect.On the other hand, our extrapolations are based on a self-consistent model and are, therefore, better suited for extrapolations to higher energies.
The flavor ratio in the lower panels is less sensitive to variations in the secondary hadron production.At low energies, where electron neutrinos originate from decaying muons, the calculations agree well since muon neutrinos and muons are both coming from pion decay.In this case the muon to electron neutrino ratio is fixed by the decay kinematics of muons.At higher energies this ratio depends more on the hadronic model since electron neutrinos have an independent production channel through decays of short and long states of neutral kaons.
The prompt flux in Sibyll-2.3cgives rise to a ratio of one at the highest energies, since the branching ratios of charmed mesons are similar for muon and electron neutrino flavors.This impacts the expected track-cascade (muon-line/electron-like events) ratio in neutrino telescopes.As the lower left panel clearly demonstrates, the flavor ratio moves towards one since the prompt muon and electron neutrino fluxes are equal.The deviation from a conventional-only hypothesis emerges at energies as low as 10 TeV and it is not very dependent on the hadronic model.Close to 100 TeV this difference is striking and must yields a sensitive observable in next generation neutrino telescopes (with larger effective areas for the cascade channel).One caveat is the presence of the astrophysical flux that is currently compatible with the same muon to electron neutrino ratio of one [5].

Vertical-to-horizontal ratio
For completeness, we discuss the vertical-to-horizontal ratio (this calculation is up-down symmetric) using the same definition as in [31] (within cos θ < 0.4 around the vertical and horizontal directions).This ratio is sensitive to the differences between 3D and 1D calculations and it is shown down to the lowest energies for both reference calculations and MCEq in Fig. 6.8.
For muon neutrinos the ratios agree within a few percent between MCEq and the Bartol calculation, which switches to full-3D at 10 GeV.At energies below 30 GeV (where the HKKMS calculations witches to 3D) there is an observable shift between MCEq and HKKMS that stays below 10%.The old Sibyll-2.1 notably disagrees at medium energies due to charged kaons and the resulting error in the production height that is different for the muon neutrinos from pion decays.Apart from this, the  Neutrino-antineutrino and flavor ratios.The left panels show zenith-averaged distributions.The right panels demonstrate the relative zenith angle dependence of each ratio for a fixed neutrino energy.These curves are normalized such that at cos θ = 0.05 the value is one.Note that the axes scales for the right panels were allowed to float.dependence on the details of the hadronic model is weak.
In electron neutrinos the differences among the hadronic models are much larger, reflecting the higher uncertainty in the production mechanisms of charged and neutral kaons that are still fundamentally not understood in particle physics, unfortunately.Due to the lack of access to the pion and kaon components of the reference calculations the exact origin of these different behaviors is difficult to trace.
The angular distribution of electron neutrinos at high energies is significantly affected by the prompt flux.The right panel shows that (even when using large zenith bins of cos θ < 0.4), the expected flux of the vertical-tohorizontal ratio is a very sensitive observable.As mentioned above, the detection of the prompt neutrino flux with a future neutrino telescope in the cascade channel would not require excessively high energies that otherwise would impact the statistical errors.The two signatures, the angular distribution of cascades and the above mentioned track/cascade ratio, are sensitive to the flux excess in the (upper) triangle enclosed by the dotted black, the solid gray and the red upper limit line in Fig. 6.4.A successful determination of the prompt flux will remain a tough experimental challenge and almost certainly require an upgraded detector and a better characterization of the astrophysical flux.

VII. SUMMARY
This work is about the connection between hadronic interactions and the inclusive fluxes of muons and neutrinos in the Earth's atmosphere.The numerical solution of the transport equations with MCEq provides a study of this connection up to the highest energies at very high precision, being essentially free from statistical simulation uncertainties.
We characterized the distributions of cosmic ray energies that can produce leptons at certain energies at ground, taking into account the effect of a non-powerlaw primary spectrum affected by the knee and the ankle of cosmic rays.Essentially the entire energy range is accessible at particle colliders, but in the center of mass frame.However, the collider kinematics and the lack of proton -light ion imposes limitations on the applicability of these data to our case.We identify all types of hadrons that contribute to the inclusive fluxes and how the interplay between the production cross section and decay time impacts the zenith distributions.
The atmospheric muons, which can be measured with high precision, behave in several ways differently from the atmospheric neutrinos.The latter can only be assessed inclusively, i.e. integrated over the primary cosmic ray energy, while muons can be studied also in exclusive events like air-showers or muon bundles.Of particular importance at high energies is the prompt flux that requires a model for the production of heavy flavor hadrons.It is an integral part of the new Sibyll-2.3cmodel and it is discussed in greater detail.The prompt flux of muons is different from the prompt neutrino flux and we characterized these contributions from unflavored mesons.The muon charge ratio is known to be sensitive to the forward particle production.To constrain the forward pion and kaon production, one needs to account both for the primary spectrum of nucleons (separating protons and neutrons) and for the shape of the particles' longitudinal energy spectra.Because the fluxes in this work are calculated with a single model of the primary spectrum, a full account of the lepton ratios needs further work.
We show that the longitudinal spectra at high x F or x Lab > 0.2 are paramount for inclusive fluxes and introduce a new non-perturbative process in the Sibyll-2.3cinteraction model that gives additional degrees of freedom to better reproduce the leading particle effect, a forward flavor asymmetry that originates from the high momentum fraction carried by the valence quarks of the projectile.This mechanism in the new version of Sibyll is the remnant excitation model, in which the valence quarks are separated from the sea partons and allowed to fragment independently.The remnant model gives a significantly improved description of fixed-target data, resulting in a better, albeit not yet perfect, prediction of the muon charge ratio.We made sure that the new version is approximately compatible with Feynman scaling in the forward phase space as observed in the data and which is a central element of the Dual Parton Model.
The new model for the production of heavy flavors in Sibyll is based on the family relation between strange and charmed quarks.At different stages of the event generation, charmed quarks can be produced through the replacement s → c with certain transition probabilities that have been determined by comparison to a large variety of data sets on production of charmed hadrons.The large mass of the charm quark is beyond the non-perturbative scale, and therefore, charm is predominantly produced in perturbative (hard) processes.We find that the augmented minijet model is sufficient to describe the total yields of charm at all accessible energies from fixed-target experiments up to the LHC.The differential (p T ) spectra are tolerably described, but show some tension towards forward LHCb rapidities.We find that the hard component is insufficient to describe charm production at fixed-target energies at forward x F , which requires accounting for processes such as associated production.At high energies, where the charm production is relevant for atmospheric neutrino fluxes, the dominant contribution comes from the perturbative component and this scenario is in agreement with other contemporary calculations of the prompt flux.In extensive comparisons with NLO calculations, we generally find an agreement between our simplified approach and the more sophisticated methods within their uncertainties.
Finally, we benchmark the combination of Sibyll-2.3cand MCEq against other reference calculations including full 3-dimensional Monte Carlo calculations.For MCEq we also employ other interaction models to better disentangle the impact of hadronic interactions from the cascade physics or the calculation method.We generally find a very good agreement when using our methods that require a tiny fraction of computational time.Most differences arise from the modeling of forward kaon production that is the most uncertain component in the prediction of atmospheric fluxes and flavor ratios.However, there are additional features in the angular distributions close to the horizon that do not seem to come from differences in hadronic interactions and more likely stem from the calculation methods.The Sibyll-2.1 model is shown to overproduce K + with a notable impact on many observables.Therefore we discourage users from employing this model in future simulations and use instead the new version or one of the other interaction models.The new charm model predicts a prompt flux that is somewhat higher than the central expectations of the other current models (within errors) and it is also compatible with the experimental limit by IceCube.We discuss prospects for measuring prompt neutrinos at current or future neutrino telescopes and outline a number of distinct signatures that can be assessed through the cascade channel at moderate energies between 10 -100 TeV.The impacted variables are the muon-to-electron neutrino ratio and the vertical-to-horizontal ratio that are both sensitive to the angular distribution and the track-to-cascade ratio in volumetric detectors.
FIG. 3.1.Upper panel: Flux of nucleons in (m 2 s sr GeV 1.6 ) −1 , converted from the flux of cosmic ray nuclei as predicted by the H3a model.Middle panel: spectral index, log-derivative of the upper curve.Lower panel: Fraction of neutrons in the nucleon flux.The spectral softening at around 1 PeV is an effect of the knee and the hardening at ∼ 100 PeV of the ankle of cosmic rays.
FIG. 4.6.Red dashed lines show the contribution to the muon spectrum by interactions of unstable particles (π ± , Λ, etc.) and solid blue the fraction of muons that is related to hadronic interactions in the first interaction length (∼ 90 g/cm 2 ) of the cosmic ray nucleon.The upper panel is computed for a single 100 PeV proton and a vertical incidence angle.

FIG. 4 . 8 .
FIG. 4.8.Zenith angle dependence of the hadronic components of the atmospheric flux in (cm 2 s sr GeV) −1 .The signature of a "prompt" component is a flat zenith distribution, meaning that these "mother-mesons" do not interact with the atmosphere before decay.At energies above the critical energy this probability increases and depends on the density gradient along the cascade trajectory.The gradient is smaller for horizontal trajectories, thus the probability for unstable particles to decay becomes higher again.The prompt component from charmed and unflavored meson decay is flat up to the highest energies due to the short life times these hadrons.
FIG. 4.9.PDF of the relevant x Lab values for inclusive leptons at fixed energy.These curves are derived from the inclusive flux, i.e. integrated over the cosmic ray spectrum.An energy independent Z-factor approach would results in a universal shape.Here, the energy dependence comes from deviations from perfect Feynman scaling and from changes in the spectral index of the primary nucleon flux due to the knee and ankle.

FIG. 5 . 6 .
FIG. 5.6.Scaling behavior of the energy distribution of secondary mesons in Sibyll-2.3c.For perfect scaling the lines for the different laboratory frame projectile energies should be on top of each other.The lowest energy line deviates due to threshold.At very small x Lab scaling violations due to multi-parton interactions are expected and visible at higher energies.For π − and K − the scaling is better than for their positively charged counterparts, which are more affected by the remnant model.

FIG. 5 . 8 .
FIG. 5.8.String configuration with remnant excitation for the lower baryon.Note that the flavor content of the remnant can be different from the initial hadron in the case of remnant excitation.

5 FIG. 5 . 9 .FIG. 5 . 10 .FIG. 5 . 11 .FIG. 5 .
FIG.5.9.Inclusive cross section of the production of charmed quarks.Data below 100 GeV are from fixed-target experiments[20,27,28,70,93].Intermediate and high energy data are from the RHIC and the LHC[2,3,[11][12][13]15]. Apart from the full cross section (magenta line), cross sections for different phase space cuts of the LHC experiments are shown.The blue line corresponds to the central phase space covered in ALICE (|y| < 0.5), while the green line corresponds to the more forward coverage of LHCb (2.5 < y < 4.5).In yellow the contribution from phase space not covered by any experiment is shown.The dashed line represents the cross section for hard parton scattering scaled to represent charm production in the model.
FIG. 5.13.pT-spectrum of neutral D-mesons in LHCb at 7 and 13 TeV c.m. energy [2, 3] compared to Sibyll 2.3c FIG. 5.14.pT-spectrum of neutral D-mesons at √ s = 13 TeV.Compared are Sibyll and a NLO QCD calculation (POWHEG)[55].The band of the theory calculation corresponds to charm quark mass and factorization scale uncertainty.The central phase space covered by the LHCb acceptance is dominated by perturbative processes, such that the full prediction by Sibyll can be compared to the pQCD calculation.

FIG. 5 . 15 .
FIG. 5.15.Longitudinal momentum spectrum of D-mesons at √ s = 13 TeV.Sibyll is compared with a calculation in perturbative QCD using extensions to small-x (GBW)[57,58,71].GBW model is scaled to match Sibyll at xF ∼ 0.4.While the yields are different, the shape of the spectrum predicted for the perturbative component is comparable.
FIG. 5.16.Comparison of the energy spectra of mesons between Sibyll (red lines) and NLO QCD calculations (PROSA 2017[54] (black lines) and GMVF[33] (purple lines)).Energy is evaluated in the lab.frame.The lower panel shows the difference relative to Sibyll.The grey band represents the theoretical uncertainty in the PROSA calculation (factorization and renormalization scale).Calculations for pp interactions are shown by solid lines and those for p-air with are dashed or dotted.The NLO calculations for an Air target are scaled up according to 14.5×pp.

FIG. 5 . 17 .
FIG.5.17.Weighted energy spectrum of D-mesons at √ s = 13 TeV.The energy of the particles is taken in the lab.frame and expressed relative to the beam energy.Contributions by different phase space intervals are shown.The contribution of the phase space covered by LHCb in events that also trigger LHCb is around 3 %.

FIG. 6 . 1 .
FIG. 6.1.Fractions of hadrons decaying into atmospheric leptons.The muon flux (upper panel) is calculated with the indicated interaction models for θ = 0 • .The neutrino fluxes (lower panels) show the results for zenith averaged fluxes.

FIG. 6 . 2 .
FIG.6.2.The atmospheric muon flux and the charge ratio.The MCEq fluxes for different interaction models are calculated for the zenith angle θ = 0 • and the H3a primary spectrum.For the Kochanov et al.[66] curves the authors used the KM parameterization for hadronic interactions and the Zatsepin-Sokolskaya primary spectrum.

FIG. 6 . 3 .
FIG.6.3.Atmospheric muon and neutrino fluxes averaged over the zenith angle.The HKKMS and Bartol curves are from computations for Kamioka site and Solar Flux minimum.At energies above a few GeV the dependence on the detector location and solar modulation diminishes.The curves by Sinegovskaya et al. are computed for the similar choice of models as in our case, namely an H3a primary spectrum and QGSJet-II-03 interaction model.
FIG.6.5.Azimuth-averaged zenith distributions at fixed neutrino energies.In the upper panels, each individual curve is normalized to one at cos θ = 0.55 and offset through a multiplication with 2.5.The lower panels show the model ratios, normalized at cos θ = 0.55 to the Sibyll-2.3cvalue.

Fig. 4 .
Fig. 4.8 the angular spectrum encodes the presence of different hadronic species and it is instructive to involve that figure in the present discussion.

FIG. 6 . 6 .
FIG. 6.6.Azimuth-averaged zenith distributions of the total (conv.+ prompt) atmospheric muon flux at high energies.The angle cos θ = 0 corresponds to horizontal directions and cos θ = 1 to vertical down-going.The left panels show the muon fluxes computed with different hadronic models.Their ratios to Sibyll-2.3care located on the right hand side and normalized to one at cos θ = 0.55.Note the increased scale of the lowest right panel.
FIG. 6.7.Neutrino-antineutrino and flavor ratios.The left panels show zenith-averaged distributions.The right panels demonstrate the relative zenith angle dependence of each ratio for a fixed neutrino energy.These curves are normalized such that at cos θ = 0.05 the value is one.Note that the axes scales for the right panels were allowed to float.

FIG. 6 . 8 .
FIG.6.8.Down to horizontal ratio.It is defined as the ratio of fluxes integrated in the angular bin cos θ < 0.4 around the two extreme directions.The comparison is shown down to the lowest energies, since it is stongly impacted by 3D calculations at energies below a few GeV.Large differences at high energy originate from the presence of prompt neutrinos in the Sibyll-2.3ccalculations.
FIG.4.1.Probability density functions (PDFs) of the primary nucleon energies corresponding to inclusive leptons a a fixed energy (colors).The solid, dashed and dotted lines refer to the individual lepton species.Atmospheric neutrinos up to 1 PeV probe center of mass collisions (top axis) that are within reach of current collider experiments (the LHC reaches 13 TeV). e

TABLE I .
Table of the free model parameters that control the probabilities for charmed quarks production in different processes.As described in more detail in the text, the effective rates can be attenuated by additional factors.

TABLE II .
Experiments that collected data on charm production including the corresponding projectile-target configuration and the accessible longitudinal phase space.These data have been used for model development and parameter estimation.
. 6.4.Prompt atmospheric muon and electron neutrino fluxes averaged over zenith angles.All fluxes are computed for the H3a cosmic ray flux model and the most differences arise from the charm cross section calculations.Sibyll-2.3c(pQCD) is the component attributed to the perturbative * component from Figs. 5.17 and 5.16.The pQCD calculations BERSS FIG