Hadronic interaction model sibyll $2.3\mathrm{c}$ and inclusive lepton fluxes

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

Figure 1

Decay lengths for relevant mesons (solid) compared to the interaction length for kaons calculated at 8 km altitude a.s.l.

Figure 2

Upper panel: Flux of nucleons in $({\mathrm{m}}^2\mathrm{s}\mathrm{sr}{\mathrm{GeV}}^{1.6}{)}^{-1}$, converted from the flux of cosmic ray nuclei as predicted by the H3a model. Middle panel: spectral index, log-derivative of the upper curve. Lower panel: Fraction of neutrons in the nucleon flux. The spectral softening at around 1 PeV is an effect of the knee and the hardening at $\sim 100\mathrm{PeV}$ of the ankle of cosmic rays.

Figure 3

Probability density functions (PDFs) of the primary nucleon energies corresponding to inclusive leptons a a fixed energy (colors). The solid, dashed and dotted lines refer to the individual lepton species. Atmospheric neutrinos up to 1 PeV probe center of mass collisions (top axis) that are within reach of current collider experiments (the LHC reaches 13 TeV).

Figure 4

Like Fig. 3 but as a function of fractional energy $E_{\text{cosmic ray}}/E_{\text{lepton}}$, for muons (upper panel), muon neutrinos (middle panel) and electron neutrinos (lower panel).

Figure 5

Contribution from decays of various particles to the atmospheric ${\mu }^{+}+{\mu }^{-}$ (top left), ${\nu }_{\mu }+{\overline{\nu }}_{\mu }$ (top right), ${\nu }_e+{\overline{\nu }}_e$ (bottom left) and ${\nu }_{\tau }+{\overline{\nu }}_{\tau }$ (bottom right) flux in sibyll-2.3c and H3a primary model at $\theta =60°$.

Figure 6

Breakdown of the unflavored component of the prompt muon flux.

Figure 7

Muons charge ratio and effective components from pion and kaon decays. (Upper panel) The solid line is the total muon charge ratio from sibyll and the areas the individual weights from the pion and kaon component. (Lower panel) The colored areas are a measure of the effective weight of the pion and the kaon component. The variations at 1 PeV come from the change of the spectral index at the knee, resulting in a different weight (or mean) of the $x_{\text{Lab}}$ distribution of secondaries.

Figure 8

Red dashed lines show the contribution to the muon spectrum by interactions of unstable particles (${\pi }^{\pm }$, $\mathrm{\Lambda }$, etc.) and solid blue the fraction of muons that is related to hadronic interactions in the first interaction length ($\sim 90\mathrm{g}/{\mathrm{cm}}^2$) of the cosmic ray nucleon. The upper panel is computed for a single 100 PeV proton and a vertical incidence angle.

Figure 9

Stacked plot of the feed-down contributions of short-lived hadrons to the inclusive pion, kaon, and proton spectra, computed with sibyll at $E_{\mathrm{lab}}=1\mathrm{TeV}$ for proton-carbon interactions. The hadrons produced directly in string fragmentation or resonance decay are normalized to one (dashed line). Hadrons that contribute below a few percent are omitted.

Figure 10

Zenith angle dependence of the hadronic components of the atmospheric flux in $({\mathrm{cm}}^2\mathrm{s}\mathrm{sr}\mathrm{GeV}{)}^{-1}$. The signature of a “prompt” component is a flat zenith distribution, meaning that these “mother-mesons” do not interact with the atmosphere before decay. At energies above the critical energy $\epsilon$ this probability increases and depends on the density gradient along the cascade trajectory. The gradient is smaller for horizontal trajectories, thus the probability for unstable particles to decay becomes higher again. The prompt component from charmed and unflavored meson decay is flat up to the highest energies due to the short life times these hadrons.

Figure 11

PDF of the relevant $x_{\text{Lab}}$ values for inclusive leptons at fixed energy. These curves are derived from the inclusive flux, i.e., integrated over the cosmic ray spectrum. An energy independent Z-factor approach would results in a universal shape. Here, the energy dependence comes from deviations from perfect Feynman scaling and from changes in the spectral index of the primary nucleon flux due to the knee and ankle.

Figure 12

NA49 pions in sibyll-2.1 (left) and sibyll-2.3c (right) [35, 36] for proton-proton (solid) and proton-carbon (dashed) interactions. In the fragmentation region ${\pi }^{+}$ ($\mathrm{u}\overline{\mathrm{d}}$) dominate over ${\pi }^{-}$ ($\mathrm{d}\overline{\mathrm{u}}$) because of the quark content of the proton (uud). The remnant model in sibyll-2.3c is tuned to reproduce these data in the forward region.

Figure 13

Longitudinal momentum spectrum of charged kaons measured in pp fixed-target collisions by NA49 [37], compared with sibyll.

Figure 14

NA49 measurement of the charge ratio of pions and kaons in pp interactions as a function of Feynman-$x$ [35, 36, 37], compared with sibyll.

Figure 15

Longitudinal momentum spectrum of the of neutral pions at 7 TeV. Data are from the LHCf spectrometer, covering $8.9<\eta <10.2$ in two different $p_T$ bins.

Figure 16

Yield of neutral pions as a function of the rapidity loss $y_{\mathrm{loss}}=y_{\mathrm{beam}}-y$ at different c.m. energies [38, 39]. If Feynman scaling holds in the fragmentation region, the spectrum is universal.

Figure 17

Scaling behavior of the energy distribution of secondary mesons in sibyll-2.3c. For perfect scaling the lines for the different laboratory frame projectile energies should be on top of each other. The lowest energy line deviates due to threshold. At very small $x_{\text{Lab}}$ scaling violations due to multiparton interactions are expected and visible at higher energies. For ${\pi }^{-}$ and $K^{-}$ the scaling is better than for their positively charged counterparts, which are more affected by the remnant model.

Figure 18

String configuration for minijets in sibyll.

Figure 19

String configuration with remnant excitation. Note that the flavor content of the remnant can be different from the initial hadron in the case of remnant excitation.

Figure 20

Inclusive cross section of the production of charmed quarks. Data below 100 GeV are from fixed-target experiments [52, 53, 54, 55, 56]. Intermediate and high energy data are from the RHIC and the LHC [57, 58, 59, 60, 61, 62]. Apart from the full cross section (magenta line), cross sections for different phase space cuts of the LHC experiments are shown. The blue line corresponds to the central phase space covered in ALICE ($|y|<0.5$), while the green line corresponds to the more forward coverage of LHCb ($2.5<y<4.5$). In yellow the contribution from phase space not covered by any experiment is shown. The dashed line represents the cross section for hard parton scattering scaled to represent charm production in the model.

Figure 21

Feynman-$x$ spectrum of charmed mesons in $\pi$-nucleus interactions. Data are from E769 [63], taken with beam momentum of $250\mathrm{GeV}/\mathrm{c}$. Dashed lines are contributions from individual processes in sibyll. The “remnant” contribution arises from the exchange of a quark with the soft strings (see Fig. 19).

Figure 22

Spectrum of the transverse momentum of D-mesons. Measurement from ALICE is taken in the central phase space $|y|<0.5$ at $\sqrt{s}=7\mathrm{TeV}$ [60]. Relative abundance of D-mesons, vector resonances (${\mathrm{D}}^{\star }$) and D-mesons containing a strange quark (${\mathrm{D}}_{\mathrm{s}}$) is determined by the parameters in the strange quark sector. The line is the prediction from sibyll-2.3c.

Figure 23

Feynman-$x$ spectrum of charmed mesons in pp interactions. Data are from LEBC experiment [54, 64], recorded at beam momenta of 400 and $800\mathrm{GeV}/\mathrm{c}$. The lines are the predictions from sibyll-2.3c.

Figure 24

$p_{\mathrm{T}}$-spectrum of neutral D-mesons in LHCb at 7 and 13 TeV c.m. energy [61, 62] compared to sibyll 2.3c

Figure 25

$p_{\mathrm{T}}$-spectrum of neutral D-mesons at $\sqrt{s}=13\mathrm{TeV}$. Compared are sibyll and a NLO QCD calculation (POWHEG) [65]. The band of the theory calculation corresponds to charm quark mass and factorization scale uncertainty. The central phase space covered by the LHCb acceptance is dominated by perturbative processes, such that the full prediction by sibyll can be compared to the pQCD calculation.

Figure 26

Longitudinal momentum spectrum of D-mesons at $\sqrt{s}=13\mathrm{TeV}$. sibyll is compared with a calculation in perturbative QCD using extensions to small-$x$ (GBW) [66, 67, 68]. GBW model is scaled to match sibyll at $x_{\mathrm{F}}\sim 0.4$. While the yields are different, the shape of the spectrum predicted for the perturbative component is comparable.

Figure 27

Comparison of the energy spectra of charmed mesons between sibyll (red lines) and NLO QCD calculations [PROSA 2017 [69] (black lines) and GMVF [70] (purple lines)]. Energy is evaluated in the lab. frame. The lower panel shows the difference relative to sibyll. The grey band represents the theoretical uncertainty in the PROSA calculation (factorization and renormalization scale). Calculations for pp interactions are shown by solid lines and those for p-air with are dashed or dotted. The NLO calculations for an Air target are scaled up according to $14.5\times \mathrm{pp}$.

Figure 28

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

Figure 29

Fractions of hadrons decaying into atmospheric leptons. The muon flux (upper panel) is calculated with the indicated interaction models for $\theta =0°$. The neutrino fluxes (lower panels) show the results for zenith averaged fluxes.

Figure 30

The atmospheric muon flux and the charge ratio. The mceq fluxes for different interaction models are calculated for the zenith angle $\theta =0°$ and the H3a primary spectrum. For the Kochanov et al. [79] curves the authors used the KM parametrization for hadronic interactions and the Zatsepin-Sokolskaya primary spectrum.

Figure 31

Atmospheric muon and neutrino fluxes averaged over the zenith angle. The HKKMS and Bartol curves are from computations for Kamioka site and Solar Flux minimum. At energies above a few GeV the dependence on the detector location and solar modulation diminishes. The curves by Sinegovskaya et al. are computed for the similar choice of models as in our case, namely an H3a primary spectrum and qgsjet-ii-03 interaction model.

Figure 32

Prompt atmospheric muon and electron neutrino fluxes averaged over zenith angles. All fluxes are computed for the H3a cosmic ray flux model and the most differences arise from the charm cross section calculations. sibyll-2.3c (pQCD) is the component attributed to the ${\text{perturbative}}^{*}$ component from Figs. 28 and 27. The pQCD calculations BERSS [84], GRRST [65, 87], PROSA [69] and GM-VNFS [70] are computed at NLO precision under different assumptions for vertical angles. For $\theta =0°$ the shape of the black curves would coincide with the shape of the green curve. They are accompanied by large error bands which are omitted for clarity. Within the errors all models are compatible with each other. The bound by IceCube [86] of 1.06x is expressed in units of normalization of the model from [88], recomputed for the H3a primary flux. The red hatched patch represents the excess of the prompt flux that leads to distinct signatures in the ${\nu }_{\mu }/{\nu }_e$ and the vertical-to-horizontal ratio.

Figure 33

Atmospheric muon and neutrino fluxes averaged over the zenith angle. The figure is similar to Fig. 31, but instead of the H3a primary flux the nucleon flux from the Global Spline Fit model (GSF) [89] has been used. The differences in the shape and spectral index of the zenith-averaged neutrino flux are significantly reduced in the comparison between mceq and the HKKMS or Bartol calculations.

Figure 34

Figure identical to Figs. 31 and 33 but with the primary flux model from [76] that is used in the HKKMS and Bartol calculations.

Figure 35

Azimuth-averaged zenith distributions at fixed neutrino energies. In the upper panels, each individual curve is normalized to one at $\cos \theta =0.55$ and offset through a multiplication with 2.5. The lower panels show the model ratios, normalized at $\cos \theta =0.55$ to the sibyll-2.3c value.

Figure 36

Azimuth-averaged zenith distributions of the total ($\mathrm{conv}.+\mathrm{prompt}$) atmospheric muon flux at high energies. The angle $\cos \theta =0$ corresponds to horizontal directions and $\cos \theta =1$ to vertical down-going. The left panels show the muon fluxes computed with different hadronic models. Their ratios to sibyll-2.3c are located on the right-hand side and normalized to one at $\cos \theta =0.55$. Note the increased scale of the lowest right panel.

Figure 37

Neutrino-antineutrino and flavor ratios. The left panels show zenith-averaged distributions. The right panels demonstrate the relative zenith angle dependence of each ratio for a fixed neutrino energy. These curves are normalized such that at $\cos \theta =0.05$ the value is one. Note that the axes scales for the right panels were allowed to float.

Figure 38

Down to horizontal ratio. It is defined as the ratio of fluxes integrated in the angular bin $\cos \theta <0.4$ around the two extreme directions. The comparison is shown down to the lowest energies, since it is strongly impacted by 3D calculations at energies below a few GeV. Large differences at high energy originate from the presence of prompt neutrinos in the sibyll-2.3c calculations.