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

###### Abstract

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.3c deals with the forward fragmentation region including production of charm. In Section VI we summarize the impact of Sibyll-2.3c on 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

(1) |

where is the density of the atmosphere and 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.

### ii.1 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

(2) |

The coupled cascade equations

(3) |

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 () 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

(4) |

and decay lengths

(5) |

are shown in Fig. 2.1. For short-lived particles , 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 gaisser_book (). The transition between these two regimes (cross-over between solid and dashed lines in Fig. 2.1) is called critical energy and is a function of the density or the altitude. For typical conditions the values are approximately GeV for , GeV for , PeV for D-mesons and GeV for unflavored mesons.

Traditional semi-analytical solutions (see for in depth discussion gaisser_book (); Lipari:1993hd (); Fedynitch:2012fs ()) have the form

(6) |

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 factors are (primary cosmic ray) spectrum weighted moments, representing hadronic interactions and particle decays

(7) |

The values of the spectral index are between 2.7 and 3.0 (see Sect. III.2). The weight emphasizes the very forward part of the particle spectrum with , 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 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 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 attenuation-less 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 ) and because the leptonic decay branching ratios are small (). All other inclusive leptons (mainly from decay of pions and kaons) constitute the conventional component, which is suppressed at very high energies 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 cpc_paper (); Fedynitch:2015zbe (); Fedynitch:icrc2013 (); Fedynitch:2015zma (); Fedynitch:2016nup (). 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 Fedynitch:2012fs (); Honda:2015fha () 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.

### iii.1 Matrix cascade equations

The open-source software Matrix Cascade Equations (MCEq)^{1}^{1}1https://github.com/afedynitch/MCEq cpc_paper () solves the system of discrete coupled cascade equations

(8) |

The index represents one of the particle species and the energy index runs over an energy grid, subdivided into 8 logarithmically spaced bins per decade of energy across 11 orders of magnitude (GeV - EeV). The coefficients represent inclusive secondary particle energy spectra in the target laboratory frame, and the corresponding decay spectra from Eq. (3). The solver performs a simultaneous integration of a coupled system of 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

(9) |

The matrices and contain the coefficients and 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 pythia (); Sjostrand:2015cx (). The computation is significantly accelerated by using sparse matrices and a method to reduce the stiffness cpc_paper (). 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.TheIceCube:2016oqi (); Arguelles:2017eao (); Edsjo:2017kjk ().

### iii.2 Primary cosmic ray flux

For the discussion in this paper we choose one representative realistic model (called H3a Gaisser:2012em ()) 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 is modeled as a broken power-law, where each source class accelerates nuclei to a maximal cut-off rigidity

(10) |

Our calculation method uses the superposition approximation, where the flux of interacting nuclei is expressed in terms of the individual nucleons (all-nucleon flux)

(11) |

or more precisely an energy dependent proton and neutron flux that can be obtained by substituting with or , respectively.

Fig. 3.1 summarizes the different properties of this model.

## 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 Sanuki:2006yd (). Here we revisit this topic using the new Sibyll-2.3c event generator and the efficient numerical scheme, extending it to higher energies and to prompt leptons. All Figures and computations are made with Sibyll-2.3c if not otherwise noted.

### iv.1 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 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 ). 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 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 and . 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 at intermediate energies.

It is interesting to note that the corresponding center-of-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 or ^{2}^{2}2Feynman- . Another obstacle is that nuclear interactions at TeV energies are currently probed in lead-lead or proton-lead collisions at the LHC, but lighter ions in the mass range of air molecules are accessible only in simulations.

### iv.2 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 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 baryon. Inclusive -integrated cross sections are computed for each hadronic species irrespective of its life time ^{3}^{3}3 is defined in the laboratory frame as . Decays are tabulated separately by using pythia 8 pythia (); Sjostrand:2015cx (). The inelastic interaction cross sections of more exotic hadrons are assumed to be equal to for all baryons, for pions and light mesons, and 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 illana_eta_2010 (). 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 , breaking the correlation between prompt fluxes of muons and neutrinos. The cross-over 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 photo-production of muon pairs, which is suppressed by wrt. the pair production cross section Nelson:1985ec (), 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 GeV the main source of muon neutrinos (upper right panel) are semi-leptonic and 3-body decays of charged kaons, see e.g. Honda_1D_1995 () 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 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 Bulmahn:2010pg (), but we can discuss their flux for completeness. The dominant production channel of tau neutrinos is the decay of , where the subsequent decay of 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% 2003AcPPB..34.3273M (); 2009JCAP…09..008I ().

### iv.3 Muon charge ratio

The inclusive muon charge ratio has been perceived in the literature as the observable with one of the highest sensitivity to the hadronic physics in atmospheric cascades Gaisser:2012em (); Adamson:2007ww (); Agafonova:2014mzx (). 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. 4.5 as a transition from an energy range, dominated by the 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 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 at a fixed covers a range of primary interaction energies (compare with Fig. 4.1), the forward particle spectra approximately scale (see Section V.1 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 (see Fig. 5.3), 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 together with the inclusive muon neutrino and electron neutrino fluxes.

### iv.4 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.

### iv.5 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 ( 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 Drescher:2007hc (); Ostapchenko:2013pia (), leading meson in pion-nucleus interactions (instead of leading 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 Bengtsson:1987kr () 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 Sjostrand:2015cx (), at least when sticking to the fragmentation of color strings as the underlying picture for hadron formation.

### iv.6 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 gaisser_book (). 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 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 and do not accumulate enough depth. The distinctive feature of the prompt component is the flat angular distribution. For fluxes from charmed meson decays, re-interactions with the atmosphere become relevant above 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).

### iv.7 Particle production phase space

In section IV.1, 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 or very approximately at

(12) |

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 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 values commonly refer to all hadrons including (p, n, , K).

A value of 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 two-component dual parton model with soft and hard minijets Capella:1992yb (). It also includes low- and high mass diffraction and a model for the excitation of beam remnants Riehn:2015oba (); Ahn:2009wx (). The hadronization model is based on the Lund string fragmentation model Bengtsson:1987kr (); Sjostrand:1987xj (). 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 . Saturation is included by means of an energy dependent increase of the 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 Donnachie:1992ny ().

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.1) and the charm model is covered in Sect. V.2.

### v.1 Leading particles

While the particle production in forward phase space () plays an important role for the fluxes of leptons in the atmosphere (see Sect. IV.7), it is even more important for the charge and flavor ratios (see Sect. IV.3). 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 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

(13) |

and the di-quark (quark or antiquark for mesons) is assigned the remaining momentum . 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.3c the hard fragmentation is replaced by a beam remnant model Drescher:2007hc (); Liu:2003we (); Sjostrand:2004pf (), 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 and the excitation mass is distributed where the lower limit is the hadron mass and the maximal mass is . 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

with . 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 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 Alt:2005zq (); Alt:2006fr (); Anticic:2010yg (). The new model describes the pion charge ratio well. For kaons the description has improved but towards large the ratio is still overestimated.

Unfortunately there are no data on leading meson production and charge ratios available at 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 Adriani:2015iwv (). While the model reproduces the spectral shape in , the transverse momentum distribution is not described as well (Fig. 5.4).

LHCb Aaij:2012ut (), the detector that provides particle identification with the largest rapidity coverage at the LHC, only covers a small region of longitudinal phase space and therefore can not constrain the leading charge ratios. In a proposed fixed target configuration for LHCb Lansberg:2012wj (); Ulrich:2015uwb (), the charge ratios in the forward region could be determined for 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- with energy, the so-called Feynman scaling Feynman:1969ej (). According to measurements the scaling hypothesis holds at least up to TeV Pare:1989mr (); Adriani:2015iwv (). 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 GeV to 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 Ostapchenko:2016ytp (). With the new remnant model Sibyll initially showed a hardening of the meson spectra with interaction energy^{4}^{4}4intermediate version Sibyll 2.3. This behavior was found to originate from the splitting of leading di-quarks, originally introduced to improve the leading flavor ratios Riehn:2017mfm (). With Sibyll-2.3c scaling is approximately fulfilled.

The effects of these changes on the atmospheric flux of muons and the charge ratio of muons are discussed in Sect. IV.3.

### v.2 Charm production

#### v.2.1 Charm model

Despite their low production yield, charmed hadrons play an important role for the prompt flux of high energy leptons in the atmosphere (Sect. IV.2 and Sect. VI.1.3). We have therefore included them in the model. The production of heavy flavors in hadron interactions is qualitatively different from the production of light flavors. While the light u,d and s quarks are abundantly produced in soft fragmentation processes, the production of heavy quarks due to the large mass ( 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 () at low energy Aitala:1996hf (); Garcia:2001xj () indicate that there is a soft (non-perturbative) component frixione1994charm (); Frixione:1997ma (). Associated production of charm, (e.g. production of ) is illustrated in Fig. 5.8.

The model of charm quark production in Sibyll-2.3c uses the family connection between strange and charmed hadrons, exchanging s with c quarks in the fragmentation step Ahn:2011wt (). The total rate of charm production in this model is set by the probability 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 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 GeV is the effective charm mass. The fragmentation function used for charm quarks is the SLAC/Peterson parameterization with Peterson:1982ak (). The transverse momentum of the charm pair in the string is taken from an exponential transverse mass distribution with energy dependent mean

where is GeV for charmed mesons and GeV for charmed baryons.

The comparison of the inclusive 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 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 Aaij:2013mga (); Aaij:2015bpa () and ALICE ALICE:2011aa (), 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 for the gluon splitting (), 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 and its antiparticle observed in pp scattering by the SELEX collaboration, indicate the presence of a non-perturbative, soft component Alves:1996rz (); Garcia:2001xj (). 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 dissociation. 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 are determined by these processes. Fig. 5.10 shows how the different processes form the D-meson spectrum in pion-nucleus interactions at the E769 experiment. While minijets below the hard scale in the framework of the model are referred to as soft, they are still within the range of perturbative calculations. Correspondingly contributions from soft and hard minijets in the model are shown together as . Truly non-perturbative contributions are D mesons from the fragmentation of the valence quarks, the diffractive contribution and the remnant (all labeled accordingly). In addition to the rate of charm from string fragmentation and in minijets, the model generates charm in the splitting of the soft gluons () and as an intrinsic charm content in the remnant (). The sources for the leading charm production in the soft processes are sketched in Fig. 5.8.

The production of D-mesons in the central region as a function of the transverse momentum is shown in Fig. 5.11, compared with the measurement from ALICE at TeV c.m. energy ALICE:2011aa (). Since the central region is dominated by hard scattering, the charm rate for the hard minijets () 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 TeV is overestimated, and the shape of the spectra deviates at higher values and at large rapidities. The same trend continues at 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. 1. An overview of the available measurements that have been used for the determination of free model parameters is in Table 2.

parameter | value |
---|---|

perturbative | |

0.08 | |

non-perturbative | |

0.004 | |

0.002 | |

0.0 | |

0.004 |

Name | (GeV) | (GeV) | spectrum | coverage | Beam config. | Ref. |
---|---|---|---|---|---|---|

E-769 | yes | p-Nuc | Alves:1996rz (); Alves:1996qz () | |||

EHS | yes | p-p | AguilarBenitez:1988sb (); AguilarBenitez:1987rc () | |||

MPS | yes | p-p | Ammar:1988ta () | |||

HERA-B | no | p-Nuc | Zoccoli:2005yn () | |||

STAR | TeV | no | p-p | Adamczyk:2012af () | ||

PHENIX | TeV | no | p-p | Adare:2010de () | ||

ALICE | PeV | TeV | no | p-p | Abelev:2012vra () | |

PeV | TeV | no | p-p | ALICE:2011aa () | ||

LHCb | PeV | TeV | no | p-p | Aaij:2013mga () | |

PeV | TeV | no | p-p | Aaij:2015bpa () |

#### v.2.2 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 Gauld:2015yia (). The spectra match nicely within the theory uncertainty, indicating that our phenomenological charm model is sufficient 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- / rapidity are difficult as they include scatterings small- (, , are the momentum fractions of the partons) and typically require additional assumptions. Fig. 5.15 shows the spectrum of D-mesons at TeV, comparing Sibyll and a small- color-dipole calculation including saturation effects (GBW) Martin:2003us (); GolecBiernat:1998js (); GolecBiernat:1999qd (). In the figure the GBW calculation is re-scaled to match our spectrum at . Similar to the POWHEG calculation (Fig. 5.14), the 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 . Note that a weight of 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 Garzelli:2016xmx () and GMVF Benzke:2017yjn (). 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 Sjostrand:2015cx (). 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- region of the parton distribution functions, instead they are fixed by and extrapolated from HERA and LHCb data Zenaiev:2015rfa (). The transition from proton-proton to proton-air interactions is simplified through a superposition model . Within the theory uncertainty (see bands in the figure) which contains contributions from the variation of the renormalization and factorization scales, the -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 . Since the superposition ansatz is used for PROSA and GMVF, the difference with Sibyll increases for proton-air.

Finally, Fig. 5.17 shows the spectrum of D-mesons as a function of the energy fraction in the lab. frame, also weighted by , and broken down into different production processes. Their relative contribution to the spectrum weighted moment (-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 -factor at high energy is only of the order of %. Minijet charm production, in contrast, sums up to %, resulting in a seemingly well determined the -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 -factor will be not quite as strong. Quantitatively, while % of charm is produced in perturbative processes in the model, just % 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 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 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.2.3 and VI.2.4).

At lower energies on the other hand the fixed-target configuration of LHCb Lansberg:2012wj (); Ulrich:2015uwb () 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-proton, or ideally oxygen-proton, collisions with LHC beams (TeV) could determine the -spectrum in the intermediate energy range of GeV. For comparison, the current highest energy measurement of the longitudinal spectrum of D-mesons is at 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 Pierog:2013ria () 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 Ostapchenko:2010vb (). 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.3c calculations when doing this.

In addition, we compare to the very detailed predictions of atmospheric neutrinos from the two state-of-the-art calculations; HKKMS 2015 Honda:2015fha () 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 Barr:2004br () 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 Gaisser:2002jj (). 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 dpmjetIII (). Since DPMJET fails to describe atmospheric muon fluxes to the required precision, some corrections derived from muon measurements have been applied Honda:2006qj (). 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.1a Monte 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 . The authors of Barr:2004br () 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 Kochanov:2008pt () and neutrino flux Sinegovskaya:2014pia () calculations that rely on numerical methods Naumov:1984jfv (); Naumov:1993yr (). 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.2 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 GeV. The characterization of fluxes below 5 GeV has been a very active topic in the past (see for example the review Gaisser:2002jj () or the references in Barr:2004br (); Honda:2015fha ()) and goes beyond the scope of this work.

### vi.1 Fluxes

#### vi.1.1 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.3, Sibyll-2.3c and 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.

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 Gaisser:2012em () (see also Agafonova:2014mzx ()). 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.

#### vi.1.2 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 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 Sjostrand:2015cx () 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. Sinegovskaya:2014pia () (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.

#### vi.1.3 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 Bhattacharya:2015ds (); Gauld:2015uq (); Garzelli:2016xmx (). An advantage of the charm model implemented in a Monte Carlo event generator is the integration into air-shower codes like CORSIKA CORSIKA_report (), 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.2.2, the phenomenological approach to heavy flavor production in Sibyll-2.3 yields 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 Aartsen:2016xlq (). 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.2.

The prompt electron neutrinos dominate over the conventional at much lower energies. For up-going neutrinos this transition can occur as low as a few tens of TeV. For muon neutrinos this transition happens in the PeV range due to the higher conventional component. The prompt flux is identical for both neutrino flavors, since the leptonic branching ratios of D mesons are almost equal for electron and muon flavors. Atmospheric tau neutrinos are suppressed by one order of magnitude (see Fig. 4.3). Leptons from decay of B mesons do not constitute more than 10% of the prompt flux Martin:2003us ().

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.2.3 and VI.2.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 and Aartsen:2016xlq (). 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 Aartsen:2016xlq (). 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.2, 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.

### vi.2 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 Kajita:2016vhj (). 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 Aartsen:2017nmd (). At much higher energies, the pattern in the zenith-energy plane gives access to physics beyond standard model (BSM) scenarios TheIceCube:2016oqi (). Future experiments Aartsen:2014oha (); Adrian-Martinez:2016fdl () 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.3c against the available reference models.

#### vi.2.1 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 unpolarized 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 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 Barr:2004br (); Honda:2006qj ()) 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.

#### vi.2.2 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 Abbasi:2012kza (); Aartsen:2015nss (). 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.1.3 for neutrinos, the presence of a prompt component affects in particular the angular distributions making them flatter than the conventional-only scenario.

#### vi.2.3 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 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.

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 Gaisser:2012em () 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, somewhat higher than the corresponding ratio from Sibyll-2.3c in Fig. 6.2. The corresponding muon-neutrino/anti-neutrino 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.3c gives 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 Aartsen:2015ivb ().

#### vi.2.4 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 Barr:2004br () (within 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 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 ), the expected flux of the vertical-to-horizontal 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-power-law 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.3c model 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 or are paramount for inclusive fluxes and introduce a new non-perturbative process in the Sibyll-2.3c interaction 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 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 () 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 , 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.3c and 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.

###### Acknowledgements.

For their steady interest, the inspiring discussions and for the early adoption of our codes, we would like to thank our colleagues from the IceCube Collaboration, in particular Carlos A. Arguëlles, Dennis Soldin, Summer Blot, Jakob van Santen and Tyce de Young. We also thank Maria Vittoria Garzelli for valuable discussions and for in-depth comparisons between the charm models. This project has received funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme (Grant No. 646623). Work of T. Gaisser and T. Stanev is supported in part by grants from the U.S. Department of Energy (DE-SC0013880) and the U.S. National Science Foundation (PHY 1505990).## References

- (1) R. Aaij et al. Measurement of prompt hadron production ratios in collisions at 0.9 and 7 TeV. Eur. Phys. J., C72:2168, 2012, 1206.5160.
- (2) R. Aaij et al. Prompt charm production in pp collisions at sqrt(s)=7 TeV. Nucl.Phys., B871:1–20, 2013, 1302.2864.
- (3) R. Aaij et al. Measurements of prompt charm production cross-sections in collisions at TeV. JHEP, 03:159, 2016, 1510.01707.
- (4) M. G. Aartsen et al. Letter of Intent: The Precision IceCube Next Generation Upgrade (PINGU). 2014, 1401.2046.
- (5) M. G. Aartsen et al. Flavor Ratio of Astrophysical Neutrinos above 35 TeV in IceCube. Phys. Rev. Lett., 114(17):171102, 2015, 1502.03376.
- (6) M. G. Aartsen et al. Characterization of the Atmospheric Muon Flux in IceCube. Astropart. Phys., 78:1–27, 2016, 1506.07981.
- (7) M. G. Aartsen et al. Observation and Characterization of a Cosmic Muon Neutrino Flux from the Northern Hemisphere using six years of IceCube data. Astrophys. J., 833(1):3, 2016, 1607.08006.
- (8) M. G. Aartsen et al. Searches for Sterile Neutrinos with the IceCube Detector. Phys. Rev. Lett., 117(7):071801, 2016, 1605.01990.
- (9) M. G. Aartsen et al. Measurement of Atmospheric Neutrino Oscillations at 6–56 GeV with IceCube DeepCore. Phys. Rev. Lett., 120(7):071801, 2018, 1707.07081.
- (10) R. Abbasi et al. Lateral Distribution of Muons in IceCube Cosmic Ray Events. Phys. Rev., D87(1):012005, 2013, 1208.2979.
- (11) B. Abelev et al. Measurement of charm production at central rapidity in proton-proton collisions at tev. JHEP, 1201:128, 2012, 1111.1553.
- (12) B. Abelev et al. Measurement of charm production at central rapidity in proton-proton collisions at TeV. JHEP, 1207:191, 2012, 1205.4007.
- (13) L. Adamczyk et al. Measurements of and Production in Collisions at GeV. Phys.Rev., D86:072013, 2012, 1204.4244.
- (14) P. Adamson et al. Measurement of the atmospheric muon charge ratio at TeV energies with MINOS. Phys. Rev., D76:052003, 2007, 0705.3815.
- (15) A. Adare et al. Heavy Quark Production in and Energy Loss and Flow of Heavy Quarks in Au+Au Collisions at GeV. Phys.Rev., C84:044905, 2011, 1005.1627.
- (16) S. Adrian-Martinez et al. Letter of intent for KM3NeT 2.0. J. Phys., G43(8):084001, 2016, 1601.07459.
- (17) O. Adriani et al. Measurements of longitudinal and transverse momentum distributions for neutral pions in the forward-rapidity region with the LHCf detector. Phys. Rev., D94(3):032007, 2016, 1507.08764.
- (18) N. Agafonova et al. Measurement of the TeV atmospheric muon charge ratio with the complete OPERA data set. Eur. Phys. J., C74:2933, 2014, 1403.0244.
- (19) M. Aguilar-Benitez et al. Meson Production From 400 GeV/ Interactions. Phys. Lett., B189:476, 1987. [Erratum: Phys. Lett.B208,530(1988)].
- (20) M. Aguilar-Benitez et al. Charm hadron properties in 400-gev/c p p interactions. Z.Phys., C40:321, 1988.
- (21) E.-J. Ahn, R. Engel, T. K. Gaisser, P. Lipari, and T. Stanev. Cosmic ray interaction event generator sibyll 2.1. Phys. Rev. D, 80:094003, 2009.
- (22) E.-J. Ahn, R. Engel, T. K. Gaisser, P. Lipari, and T. Stanev. Sibyll with charm. 2011, 1102.5705.
- (23) E. M. Aitala et al. Asymmetries between the production of and mesons from 500-gev/c - nucleon interactions as a function of and . Phys.Lett., B371:157–162, 1996, hep-ex/9601001.
- (24) C. Alt et al. Inclusive production of charged pions in p p collisions at 158-gev/c beam momentum. Eur. Phys. J., C45:343–381, 2006.
- (25) C. Alt et al. Inclusive production of charged pions in p + c collisions at 158-gev/c beam momentum. Eur. Phys. J., C49:897–917, 2007.
- (26) G. A. Alves et al. Feynman x and transverse momentum dependence of D meson production in 250-GeV pi, K and p - nucleon interactions. Phys. Rev. Lett., 77:2392–2395, 1996.
- (27) G. A. Alves et al. Forward cross-sections for production of d+, d0, d/s, d*+ and lambda/c in 250-gev pi+-, k+-, and p nucleon interactions. Phys. Rev. Lett., 77:2388–2391, 1996.
- (28) R. Ammar et al. D-meson production in 800-gev/c p pinteractions. Phys.Rev.Lett., 61:2185–2188, 1988.
- (29) T. Anticic et al. Inclusive production of charged kaons in p+p collisions at 158 GeV/c beam momentum and a new evaluation of the energy dependence of kaon production up to collider energies. Eur.Phys.J., C68:1–73, 2010, 1004.1889.
- (30) C. A. Argüelles, G. de Wasseige, A. Fedynitch, and B. J. P. Jones. Solar Atmospheric Neutrinos and the Sensitivity Floor for Solar Dark Matter Annihilation Searches. JCAP, 1707(07):024, 2017, 1703.07798.
- (31) G. D. Barr, T. K. Gaisser, P. Lipari, S. Robbins, and T. Stanev. A Three - dimensional calculation of atmospheric neutrinos. Phys. Rev., D70:023006, 2004, astro-ph/0403630.
- (32) H.-U. Bengtsson and T. Sjostrand. The Lund Monte Carlo for Hadronic Processes: Pythia Version 4.8. Comput. Phys. Commun., 46:43, 1987.
- (33) M. Benzke, M. V. Garzelli, B. Kniehl, G. Kramer, S. Moch, and G. Sigl. Prompt neutrinos from atmospheric charm in the general-mass variable-flavor-number scheme. JHEP, 12:021, 2017, 1705.10386.
- (34) A. Bhattacharya, A. Stasto, I. Sarcevic, R. Enberg, and M. H. Reno. Perturbative charm production and the prompt atmospheric neutrino flux in light of RHIC and LHC. JHEP, 1506(6):110, Feb. 2015, 1502.01076.
- (35) A. Bulmahn and M. H. Reno. Secondary atmospheric tau neutrino production. Phys. Rev., D82:057302, 2010, 1007.4989.
- (36) A. Capella, U. Sukhatme, C.-I. Tan, and J. Tran Thanh Van. Dual parton model. Phys. Rept., 236:225–329, 1994.
- (37) A. Donnachie and P. V. Landshoff. Total cross-sections. Phys. Lett., B296:227–232, 1992.
- (38) H.-J. Drescher. Remnant Break-up and Muon Production in Cosmic Ray Air Showers. Phys. Rev., D77:056003, 2008, 0712.1517.
- (39) J. Edsjo, J. Elevant, R. Enberg, and C. Niblaeus. Neutrinos from cosmic ray interactions in the Sun. JCAP, 1706(06):033, 2017, 1704.02892.
- (40) R. Enberg, M. H. Reno, and I. Sarcevic. Prompt neutrino fluxes from atmospheric charm. Physical Review D, 78(4):43005, Aug. 2008.
- (41) A. Fedynitch. Phenomenology of atmospheric neutrinos. EPJ Web Conf., 116:11010, 2016.
- (42) A. Fedynitch, J. Becker Tjus, and P. Desiati. Influence of hadronic interaction models and the cosmic ray spectrum on the high energy atmospheric muon and neutrino flux. Phys. Rev., D86:114024, 2012, 1206.6710.
- (43) A. Fedynitch and R. Engel. MCEq - efficient computation of particle particle cascades in the atmosphere. 2018.
- (44) A. Fedynitch, R. Engel, T. K. Gaisser, F. Riehn, and T. Stanev. Atmospheric neutrinos at high energy. Braz. J. Phys., 44(5):pp.415–608, 2014.
- (45) A. Fedynitch, R. Engel, T. K. Gaisser, F. Riehn, and T. Stanev. Calculation of conventional and prompt lepton fluxes at very high energy. EPJ Web Conf., 99:08001, 2015, 1503.00544.
- (46) A. Fedynitch, R. Engel, T. K. Gaisser, F. Riehn, and T. Stanev. MCEQ - numerical code for inclusive lepton flux calculations. PoS, ICRC2015:1129, 2015.
- (47) R. P. Feynman. Very high-energy collisions of hadrons. Phys. Rev. Lett., 23:1415–1417, 1969.
- (48) S. Frixione, M. L. Mangano, P. Nason, and G. Ridolfi. Charm and bottom production: Theoretical results versus experimental data. Nuclear Physics B, 431(3):453–483, 1994.
- (49) S. Frixione, M. L. Mangano, P. Nason, and G. Ridolfi. Heavy quark production. Adv. Ser. Direct. High Energy Phys., 15:609–706, 1998, hep-ph/9702287.
- (50) T. K. Gaisser. Cosmic rays and particle physics. Cambridge Univ Pr, 1990.
- (51) T. K. Gaisser. Spectrum of cosmic-ray nucleons, kaon production, and the atmospheric muon charge ratio. Astroparticle Physics, 35(12):801–806, July 2012.
- (52) T. K. Gaisser and M. Honda. Flux of atmospheric neutrinos. Ann. Rev. Nucl. Part. Sci., 52:153–199, 2002, hep-ph/0203272.
- (53) F. G. Garcia et al. Hadronic production of Lambda(c) from 600-GeV/c pi-, Sigma- and p beams. Phys. Lett., B528:49–57, 2002, hep-ex/0109017.
- (54) M. V. Garzelli, S. Moch, O. Zenaiev, A. Cooper-Sarkar, A. Geiser, K. Lipka, R. Placakyte, and G. Sigl. Prompt neutrino fluxes in the atmosphere with PROSA parton distribution functions. JHEP, 05:004, 2017, 1611.03815.
- (55) R. Gauld, J. Rojo, L. Rottoli, and J. Talbert. Charm production in the forward region: constraints on the small-x gluon and backgrounds for neutrino astronomy. JHEP, 11:009, 2015, 1506.08025.
- (56) R. Gauld, J. Rojo, L. Rottoli, and J. Talbert. Charm production in the forward region: constraints on the small-x gluon and backgrounds for neutrino astronomy. JHEP, 11:009, 2015, 1506.08025.
- (57) K. J. Golec-Biernat and M. Wusthoff. Saturation effects in deep inelastic scattering at low Q**2 and its implications on diffraction. Phys. Rev., D59:014017, 1998, hep-ph/9807513.
- (58) K. J. Golec-Biernat and M. Wusthoff. Saturation in diffractive deep inelastic scattering. Phys. Rev., D60:114023, 1999, hep-ph/9903358.
- (59) D. Heck, J. Knapp, J. Capdevielle, G. Schatz, and T. Thouw. CORSIKA: A Monte Carlo Code to Simulate Extensive Air Showers. Technical Report FZKA 6019, Karsruhe, 1998.
- (60) M. Honda, M. S. Athar, T. Kajita, K. Kasahara, and S. Midorikawa. Atmospheric neutrino flux calculation using the NRLMSISE-00 atmospheric model. Phys. Rev., D92(2):023004, 2015, 1502.03916.
- (61) M. Honda, T. Kajita, K. Kasahara, and S. Midorikawa. Calculation of the flux of atmospheric neutrinos. Phys. Rev. D, 52(9):4985, Mar. 1995.
- (62) M. Honda, T. Kajita, K. Kasahara, S. Midorikawa, and T. Sanuki. Calculation of atmospheric neutrino flux using the interaction model calibrated with atmospheric muon data. Phys. Rev., D75:043006, 2007, astro-ph/0611418.
- (63) J. I. Illana, P. Lipari, M. Masip, and D. Meloni. Atmospheric muon and neutrino fluxes at very high energy. Astroparticle Physics, 34(9):663–673, Oct. 2010.
- (64) J. I. Illana, M. Masip, and D. Meloni. Atmospheric lepton fluxes at ultrahigh energies. J. Cosmol. Astropart. Phys., 09(0):008, Sept. 2009.
- (65) T. Kajita, E. Kearns, and M. Shiozawa. Establishing atmospheric neutrino oscillations with Super-Kamiokande. Nucl. Phys., B908:14–29, 2016.
- (66) A. A. Kochanov, T. S. Sinegovskaya, and S. I. Sinegovsky. High-energy cosmic ray fluxes in the Earth atmosphere: calculations vs experiments. Astropart. Phys., 30:219–233, 2008, 0803.2943.
- (67) J. P. Lansberg et al. A Fixed-Target ExpeRiment at the LHC (AFTER@LHC) : luminosities, target polarisation and a selection of physics studies. PoS, QNP2012:049, 2012, 1207.3507.
- (68) P. Lipari. Lepton spectra in the earth’s atmosphere. Astropart. Phys., 1:195–227, 1993.
- (69) F. M. Liu, T. Pierog, J. Aichelin, and K. Werner. Hadron production in proton proton scattering in NEXUS3. 2003, hep-ph/0307204.
- (70) C. Lourenco and H. K. Wohri. Heavy flavour hadro-production from fixed-target to collider energies. Phys.Rept., 433:127–180, 2006, hep-ph/0609101.
- (71) A. D. Martin, M. G. Ryskin, and A. M. Stasto. Prompt neutrinos from atmospheric and production and the gluon at very small x. Acta Phys. Polon., B34:3273–3304, 2003, hep-ph/0302140.
- (72) A. D. Martin, M. G. Ryskin, and A. M. Stasto. Prompt Neutrinos from Atmospheric cbar c Production and the Gluon at Very Small x. Acta Physica Polonica B, 34(6):3273, June 2003, hep-ph/0302140.
- (73) V. A. Naumov. A method for calculating geomagnetic corrections to differential energy spectra of cosmic-ray nucleons in the atmosphere. Res. Geomagn. Aeron. Solar Phys., 69:82–94, 1984.
- (74) V. A. Naumov, S. I. Sinegovsky, and E. V. Bugaev. High-energy cosmic ray muons under thick layers of matter. A new method for calculating the energy spectrum of cosmic ray muons under thick layers of matter. Phys. Atom. Nucl., 57:412, 1994, hep-ph/9301263.
- (75) W. R. Nelson, H. Hirayama, and D. W. O. Rogers. The Egs4 Code System. 1985.
- (76) S. Ostapchenko. Monte Carlo treatment of hadronic interactions in enhanced Pomeron scheme: I. QGSJET-II model. Phys. Rev., D83:014018, 2011, 1010.1869.
- (77) S. Ostapchenko. QGSJET-II: physics, recent improvements, and results for air showers. EPJ Web Conf., 52:02001, 2013.
- (78) S. Ostapchenko, M. Bleicher, T. Pierog, and K. Werner. Constraining high energy interaction mechanisms by studying forward hadron production at LHC. 2016, 1608.07791.
- (79) E. Pare et al. Inclusive Production of s and Feynman Scaling Test in the Fragmentation Region at the S p S Collider. Phys. Lett., B242:531–535, 1990.
- (80) C. Peterson, D. Schlatter, I. Schmitt, and P. M. Zerwas. Scaling Violations in Inclusive e+ e- Annihilation Spectra. Phys. Rev., D27:105, 1983.
- (81) T. Pierog, I. Karpenko, J. M. Katzy, E. Yatsenko, and K. Werner. EPOS LHC: Test of collective hadronization with data measured at the CERN Large Hadron Collider. Phys. Rev., C92(3):034906, 2015, 1306.0121.
- (82) F. Riehn, H. P. Dembinski, R. Engel, A. Fedynitch, T. K. Gaisser, and T. Stanev. The hadronic interaction model SIBYLL 2.3c and Feynman scaling. PoS, ICRC2017:301, 2017, 1709.07227.
- (83) F. Riehn, R. Engel, A. Fedynitch, T. K. Gaisser, and T. Stanev. A new version of the event generator Sibyll. 2015, 1510.00568.
- (84) S. Roesler, R. Engel, and J. Ranft. 2000, hep-ph/0012252.
- (85) T. Sanuki, M. Honda, T. Kajita, K. Kasahara, and S. Midorikawa. Study of cosmic ray interaction model based on atmospheric muons for the neutrino flux calculation. Phys. Rev., D75:043005, 2007, astro-ph/0611201.
- (86) T. S. Sinegovskaya, A. D. Morozova, and S. I. Sinegovsky. High-energy neutrino fluxes and flavor ratio in the Earth’s atmosphere. Phys. Rev., D91(6):063011, 2015, 1407.3591.
- (87) T. Sjostrand. Status of fragmentation models. Int. J. Mod. Phys., A3:751, 1988.
- (88) T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C. O. Rasmussen, and P. Z. Skands. An introduction to PYTHIA 8.2. CPC, 191:159–177, June 2015.
- (89) T. Sjöstrand, S. Mrenna, and P. Skands. PYTHIA 6.4 physics and manual. Journal of High Energy Physics, 2006(05):026–026, May 2006.
- (90) T. Sjostrand and P. Z. Skands. Multiple interactions and the structure of beam remnants. JHEP, 03:053, 2004, hep-ph/0402078.
- (91) R. M. Ulrich, C. Baus, R. Engel, A. Fedynitch, U. Kraemer, T. Pierog, and F. Riehn. The impact of a fixed-target experiment with LHC beam for astroparticle physics. PoS, ICRC2015:407, 2016.
- (92) O. Zenaiev et al. Impact of heavy-flavour production cross sections measured by the LHCb experiment on parton distribution functions at low x. Eur. Phys. J., C75(8):396, 2015, 1503.04581.
- (93) A. Zoccoli et al. Charm, beauty and charmonium production at HERA-B. Eur.Phys.J., C43:179–186, 2005.