Atmospheric Lepton Fluxes
This review of atmospheric muons and neutrinos emphasizes the high energy range relevant for backgrounds to high-energy neutrinos of astrophysical origin. After a brief historical introduction, the main distinguishing features of atmospheric and are discussed, along with the implications of the muon charge ratio for the ratio. Methods to account for effects of the knee in the primary cosmic-ray spectrum and the energy-dependence of hadronic interactions on the neutrino fluxes are discussed and illustrated in the context of recent results from IceCube. A simple numerical/analytic method is proposed for systematic investigation of uncertainties in neutrino fluxes arising from uncertainties in the primary cosmic-ray spectrum/composition and hadronic interactions.
Atmospheric leptons are of current interest in two contexts, as a beam for the study of neutrino oscillations and the mass hierarchy (energy range MeVTeV) and as the background in the search for high energy neutrinos of astrophysical origin (energy range GeVPeV). The emphasis of this paper is on the higher energy range, motivated by the discovery by IceCube [1, 2] of neutrinos of extraterrestrial origin at very high energy above the background of atmospheric neutrinos. In view of this discovery it is important to understand the atmospheric neutrino spectrum as precisely as possible, not only to understand the backgrounds, but also in order to learn how the spectrum of the astrophysical component extends to lower energy.
It is well known that the idea of using Cherenkov light in a large volume of water to detect interactions of high-energy neutrinos was suggested in 1960 by Markov , by Greisen  and by Reines . Markov was interested in using atmospheric neutrinos to measure the energy dependence of the neutrino cross section. He speculates about the existence of a vector boson intermediary and the question of whether there are two different kinds of neutrinos related to the electron and the muon. At the end of his review, “Cosmic Ray Showers,” Greisen speculated briefly about detecting extraterrestrial neutrinos. Reines’ review of “Neutrino Interactions” in the same volume devotes a section to “Cosmic and Cosmic Ray Neutrinos.” By “cosmic neutrinos” he means neutrinos produced by cosmic rays in or near their distant sources, while “cosmic-ray neutrinos” refers to atmospheric neutrinos produced by interactions of cosmic rays at Earth. After noting the potential of cosmic neutrinos as a probe of the origin of cosmic rays and noting the ignorance of what flux to expect, he writes, “The situation is somewhat simpler in the case of cosmic-ray neutrinos: they are both more predictable and of less intrinsic interest.” He goes on to estimate a detection rate of 1 atmospheric neutrino per day in 5000 m of water.
Perhaps less well known is the early (1961) paper  by Zatsepin and Kuz’min in which the essential phenomenology of atmospheric neutrinos is derived, including the important role of kaons relative to pions as parents of neutrinos. Their papers include the production of neutrinos from decay of muons, and they describe the strong angular dependence at high energy, which is a consequence of the dependence of the ratio of decay to interaction above the critical energies of the mesons. They also explain that charged kaons are more efficient producers of neutrinos than charged pions because of the higher mass of the kaon relative to the muon and the shorter lifetime of the kaon. In the 60’s and 70’s Volkova and Zatsepin published a series of papers refining the calculations, as described in the paper of Volkova , which remains a standard reference for fluxes of atmospheric neutrinos.
3 Overview of neutrinos and muons
The linear development of the hadronic cascade in the atmosphere is described by a system of equations of the form
where d is the flux of particles of type at slant depth in the atmosphere with energies in the interval to [8, 9]. The probability for a particle of type to interact in d is d, where is the corresponding interaction length (in g/cm). Similarly, d is the probability that a particle of type decays in d. The function is
where d is the number of particles of type i produced on average in the energy bin d around per collision of an incident particle of type . Depending on the boundary condition, the solution of Eq. 3 can describe either a single air shower or the inclusive flux of a particular particle type in the atmosphere. In the latter (inclusive) case, the boundary condition is and for . Here is the spectrum of nucleons at the top of the atmosphere. Correspondingly, the solutions give the inclusive rates of particles of type per unit area at depth independently of whether or not there are particles nearby. In contrast, for the air shower boundary condition () the solution in principle contains all the correlations among the particles in the shower.
3.1 Muon neutrinos
For the inclusive boundary condition with a power-law spectrum of primary nucleons and assuming scaling of the particle interactions in the forward kinematic region, a good approximation to the solution of Eq. 3 for the flux of is
where the sum is over the contributions of charged pions, charged kaons and charmed hadrons. At low energy () the contribution from decay of muons also needs to be added. The asterisk on cosine of the zenith angle indicates the correction necessary to account for the curvature of the Earth for . The quantity is the energy for meson above which re-interaction in the atmosphere becomes more likely than decay. At energies below the critical energy () for each channel, most mesons decay, so the angular distribution of the corresponding neutrinos is isotropic for an isotropic primary cosmic-ray spectrum. In addition, the neutrino spectrum is similar to the cosmic-ray spectrum. For the energy spectrum steepens, first near the vertical and at higher energy for more horizontal directions. As a consequence, at high energy the intensity of neutrinos from near the horizon is an order of magnitude larger that near the vertical, as illustrated in Fig. 1. The short lifetime of charmed hadrons corresponds to GeV, so the small fraction of prompt neutrinos will be isotropic below this energy.
The numerator of Eq. 3 contains the spectrum-weighted moments for production of secondary particles in particle collisions in the atmosphere. For example, for nucleon
and has a similar definition in terms of the cross section for a nucleon of energy to produce a secondary nucleon of energy . is the spectrum weighted moment of the decay distribution.
Eq. 5 is a general definition given in Ref.  that takes account of the energy dependence of interaction cross sections and particle production and allows for a general form of the primary spectrum. For a power-law primary spectrum with differential index , constant interaction cross sections and particle production that depends only on the ratio , Eq. 5 simplifies to the more familiar, energy-independent form,
3.2 Electron neutrinos
Electron neutrinos come primarily from three-body decays of kaons at energies high enough so that the contribution from muons is small. The intensity of atmospheric electron neutrinos is approximately 5% of muon neutrinos at high energy. The contribution of charm is the same for both and , which means that the fraction of prompt is significantly larger than the fraction of prompt . The lower normalization of from kaons also means that the contribution from muon decay is relatively more important to higher energy and lower zenith angle than for . A detailed calculation of electron neutrinos at high energy, including the contribution from the rare three-body decay of , is given in Ref. . Because of its short lifetime, the critical energy for is TeV. As a consequence, the spectrum of neutrinos from this channel is harder by one power of energy than that from and . Above TeV the contribution from the three channels is nearly equal because the lifetime is inversely proportional to the total decay width and the branching ratio is the ratio of the semi-leptonic width to the total width.
3.3 Muons and the ratio
The flux of muons in the atmosphere is similar in form and closely related to the flux of muon neutrinos, with three terms as in Eq. 3. The abundant and well-measured spectrum of atmospheric muons is therefore an important benchmark for any calculation of neutrinos. In fact, measurements of atmospheric muons  are used to tune the model for hadron production for one of the standard neutrino flux calculations . The main difference between and is the kinematics of the two-body decays of charged pions and kaons. Because the muon carries most of the energy in pion decay, and because more pions are produced than kaons, most muons come from the pions, even after accounting for the steep primary spectrum and the higher value of .
The signature of the kaon channel nevertheless shows up in an important way in the charge ratio of muons, which increases from around GeV  to at a TeV and above [15, 16]. Although most muons come from , the fraction from kaons increases above the energy at which charged pions begin to prefer interaction to decay in the atmosphere. The critical energy for pions is GeV compared to GeV for kaons. There is an increase in the muon charge ratio in this energy range because of the importance of associated production, , which makes the charge ratio higher for kaons than for pions. The fluxes of and can be calculated in a straightforward way by keeping track of positive and negative meson fluxes separately . The ratio depends both on the excess of protons in the primary spectrum and on the spectrum-weighted moments for the separate meson charge channels. Using the calculation of Ref.  to fit the proton excess in the primary spectrum and the Z-factor, the OPERA group found and at a mean muon energy of TeV . This result has implications for the ratio in the same energy region. With the value of from OPERA, the expected ratio for muon neutrinos increases from at low energy to above a TeV, as shown in Fig. 2.
4 Analytic and numerical calculations of atmospheric neutrinos
A formal expression for the flux of neutrinos is
where is the primary flux of nuclei with total energy and is the yield of neutrinos from a given primary nucleus. A Monte Carlo evaluation of this integral is the standard approach for detailed calculations of the fluxes of atmospheric muons and neutrinos [13, 18]. The results can be re-weighted to correspond to any desired description of the primary spectrum and composition, including the effects of geomagnetic cutoffs at different locations. It is straightforward to include muon energy loss and decay. The flux can be extended to high energy with good statistics by forcing all mesons to decay and recording the fractional probability that the decay would actually have occurred. In the calculation of Ref. , the yields were calculated by Monte Carlo for five representative nuclei (p, He, N, Si, Fe) at ten primary energies per decade, equally spaced in . The neutrino flux bins were filled with the appropriate fractional weights in logarithmically spaced bins of neutrino energies. For both Refs. [13, 18] the calculations extend only to TeV.
Analytic and numerical calculations are a useful alternative to the Monte Carlo approach because of the insight they provide into which aspects of the physics are most important for different features of the spectra of atmospheric leptons. A simple, fast analytic/numerical routine is also suitable for systematic exploration of uncertainties in the lepton fluxes due to lack of knowledge of the primary spectrum/composition and to the hadronic interactions. Two approaches were described at this conference, the numerical integration of the coupled matrix of equations  and the iterative scheme  of Sinegovskaya et al. . The approach developed in this paper uses energy-dependent Z-factors to account for energy dependence of hadronic interactions and non-power-law behavior of the primary spectrum.
4.1 Energy-dependent Z-factors
In their paper on neutrinos and muons from charm decay, Gondolo, Ingelman & Thunman  showed that the simple Eq. 3 valid for a power-law primary spectrum and scaling can be adapted to include energy dependence. With the generalized definition of the Z-factors as in the example of Eq. 5, the algebraic solutions (Eq. 3) correctly account both for the energy dependence of meson production and for the non-power-law behavior of the primary spectrum, provided that the energy dependences are sufficiently gradual. This scheme was used in Ref.  and Ref.  to account for the effect of the knee in the primary spectrum. Here I also give an example of the comparison of two different interaction models. Evaluating the energy-dependent Z-factors as in Eq. 5 requires a representation of the production distribution for each particle considered. Using the simple, two-parameter forms of Ref.  is sufficient, provided the parameters are tabulated at a sufficiently fine grid of energies. The production spectrum for is
with and the two energy-dependent parameters and . One way to characterize a hadronic event generator is by a set of spectrum weighted moments calculated at each beam energy for a range of power-law spectra. Then from
the energy-dependent parameters can be evaluated at each energy as
Having determined the parameters for particle production, The next step is to evaluate the spectrum-weighted moments at each energy for an arbitrary spectrum using
(For simplicity, the factor in Eq. 5 has been approximated as unity here because the interaction cross section varies slowly over the range of near that is important for the integral over the steep primary spectrum.)
For illustration, two hadronic interaction models are compared here. One is a very simple energy-independent extrapolation of the parameters from Ref. , which correspond to the Z-factors of Ref.  for . In this case, the only energy dependence in the is from the knee of the spectrum, as in Ref. . For comparison, the more realistic case corresponding to the hadronic interaction model of Ref.  are used. The energy-dependent Z-factors (Eq. 12) are shown in Fig. 3 (left) evaluated with the H3A primary spectrum of Ref. . The corresponding fluxes of and are shown in Fig. 3 (right), and the angular distributions at two energies are shown in Fig. 1. The analytic calculations here do not include neutrinos from decay of muons, and the tabulated data from Ref.  are above the calculation at low energy. In addition, the tabulated values were calculated with a slightly different energy spectrum, as noted in the discussion of Fig. 4 below.
4.2 Primary spectrum
For an analytic or numerical approximation as in Eq. 3 the primary spectrum enters the calculation as the spectrum of nucleons per GeV/nucleon. This superposition approximation neglects collective effects in nuclei, but correctly captures the kinematics of meson production in nucleon-nucleon collisions. We are interested in an energy range that spans direct measurements of the primary spectrum below 100 TeV/nucleon and air shower measurements to the knee and beyond. Fits to the spectrum attempt to extrapolate the direct measurements into the air-shower regime where the spectrum is generally presented as an all-particle spectrum in energy per nucleus. One common assumption in fitting the data is to assume that the primary spectrum depends only on magnetic rigidity . The rationale is that both propagation and acceleration depend on motion in magnetic fields. Several such representations of the nucleon spectrum are shown in Fig. 4.
The Polygonato model  is an example in which each element is assigned a power-law index based on direct measurements at low energy and then suppressed above an energy , where , the cutoff rigidity, is tuned to fit the knee of the cosmic-ray spectrum. In Ref. , parameters for a model (H3a) with the standard five nuclear groups (p, He, CNO, Si, Fe) and three populations are given. Two galactic populations have cutoffs of PV (for the knee) and PV. The cutoff for the extragalactic population depends on what composition is assumed. Other parameterizations (GST1 and GST2) based on a different fit to direct measurements and guided by measurements of average primary composition from air showers measurements are given in Ref. . The Polygonato model shown in Fig. 4 is a modified version with 5 nuclear components at low energy and populations 2 and 3 of H3a. The red dashed line is a version of GST2 modified by Lipari  to include an unrealistically high fraction of protons in order to estimate the maximum possible atmospheric neutrino flux above 100 TeV. The dotted line is the spectrum used in the calculation of the neutrino flux by Honda et al.  modified to include the knee as in H3a. This parameterization is the high helium version of a parameterization and extrapolation of direct measurements from Ref. . The 30% range in the knee region in Fig. 4 translates into a comparable uncertainty for the atmospheric neutrino flux above 30 TeV.
5 Neutrino yields and neutrino self-veto
The High Energy Starting Event (HESE) analysis in IceCube [1, 2] defines a veto region in the outer parts of the deep array. Events are selected that start in the fiducial volume inside the veto region. The analysis also sets a very high energy threshold by requiring an amount of light equivalent to approximately 30 TeV or more of energy deposition in the detector. A key feature of this analysis is that high-energy atmospheric neutrinos from above are excluded from the event sample if they are accompanied at the detector by a muon from the same shower that triggers the veto. Evaluation of the veto probability is an important case where the forced decay scheme for Monte Carlo evaluation of Eq. 7 at high energy does not work. The reason is that in this case it is essential to keep the correlation between the neutrino and the high-energy muons in the same shower which provide the veto.
The probability that a muon produced in the same decay as a reaches the deep detector with a minimum energy can be calculated analytically . This calculation applies only to muon neutrinos. In general, an atmospheric neutrino can also be accompanied by a muon from a different branch of the same air shower in which the neutrino was produced. Such a generalized self-veto can be applied to electron neutrinos and to prompt neutrinos of either flavor. A procedure for evaluating the generalized self-veto  consists of evaluating an expression similar to Eq. 7 but with an extra factor in the integrand:
The exponential factor is the Poisson probability that a shower generated by a primary nucleus with energy has no muons with , the minimum muon energy at production needed for the muon to reach the detector with sufficient residual energy to be detected. A new set of parameterizations analogous to the Elbert formula  was obtained from simulations with Sibyll 2.1  to use for the yields in Eq. 13. The yields for both neutrinos and muons have the integral form
In the original Elbert formula for muons, . In the expression for prompt neutrinos the factor is omitted. The differential yields are
Then the passing rate for atmospheric neutrinos of energy entering the detector with zenith angle and interacting inside is
An IceCube analysis that makes extensive use of the atmospheric neutrino self-veto was presented at this conference by J. van Santen, since published as Ref. . The importance of the veto for the discovery of neutrinos in IceCube is illustrated by Fig. 5, reproduced from the supplemental material of the IceCube publication  and shown at this conference in a presentation by G. Binder. The pink histograms in the figure show the angular distribution of prompt neutrinos from decay of charmed hadrons, which would likely be the dominant contribution to the atmospheric background for TeV. The dashed line shows the shape before the veto, while the solid line is the expectation after the veto. The angular axis is labeled in declination. At the South Pole the zenith angle is related to declination by . The veto suppresses the downward atmospheric background by 50% or more for zenith angles , a shape completely different from the data.
6 Rates in IceCube
Effective areas in the IceCube HESE analysis are given for three flavors of neutrinos separately for the Northern and Southern hemispheres in the supplementary material of Ref. . Expected rates for flavor are obtained by the convolution of the neutrino spectrum with the corresponding effective area,
The effective areas include all the details of the veto and absorption by the Earth except for the atmospheric neutrino self-veto, which does not apply to astrophysical neutrinos. For example, for muon neutrinos below TeV, the acceptance is higher for neutrinos from the Northern hemisphere because the atmospheric muon veto must be more severe for events from above the detector (Southern hemisphere). On the other hand, at high energy the acceptance is smaller from below because of absorption of neutrinos by the Earth.
Overall, the acceptance of the HESE event selection is smallest for because of the starting track criterion coupled with the high threshold for deposited energy in the detector and the muon veto. At high energy, much of the energy of a -induced muon that starts in the detector is deposited after the muon leaves the detector. For a flavor ratio of at Earth and for an astrophysical spectrum
the expected numbers of events in 988 days are . For comparison, the total number of events, including backgrounds in the 3-year analysis  is 37. The astrophysical spectrum in Eq. 17 is a fit with a differential spectral index of from Ref. .
It is instructive to investigate the distributions of neutrino energies that give rise to these astrophysical signals in comparison with the corresponding distributions for the atmospheric backgrounds. This comparison is shown in Fig. 6. The plots show the number of events per logarithmic interval of energy on a linear scale, so the area under each segment of a curve correctly represents the fraction of events from that range of energy. Another important point is that the plots show true neutrino energy whereas the events are characterized by energy deposited in the detector. Most of the atmospheric events have true neutrino energy less than TeV, while most of the astrophysical signal is above that energy. In calculating the rates of downward atmospheric neutrinos, the neutrino fluxes at production have been reduced by the passing fraction as a function of energy and zenith angle as described in Ref. . A consequence is that these backgrounds are expected to be smaller from the South than from the North, unlike the case for astrophysical neutrinos, where the opposite is true.
The atmospheric neutrino numbers include the contribution of prompt neutrinos calculated with the model of Enberg et al.  assuming a primary spectrum with the knee as in Ref. . The conventional rates are calculated with the neutrino fluxes of the energy-independent model discussed in Section 4 (the higher of the two fluxes in the right panel of Fig. 3). Most (7 of 8) of the expected atmospheric events are conventional (from decay of charged kaons and pions). For atmospheric electron neutrinos slightly more than half of an expected 3 events are prompt. The expected numbers of background neutrino events would be reduced proportionately if the lower extrapolated flux in the right panel of Fig. 3 were used, but the ratios would be essentially unchanged.
The two panels in Fig. 6 show -induced events on the left and the sum of and on the right. Most of the will be classified as tracks, but some (including all neutral current interactions) will be classified as cascades. On the other hand all and most will be classified as cascades. At present the observations are consistent with a flavor ratio of the astrophysical neutrinos at Earth. The atmospheric neutrino backgrounds should be mostly -induced. In principle, given a significantly larger event sample and a good knowledge of normalization of the background of atmospheric neutrinos at high energy, it would be possible to subtract the background and measure the astrophysical flavor ratio. At present, however, the events are too few and the uncertainties in the atmospheric background too large to do so.
7 Summary and Outlook
It is clearly important to acquire the best possible understanding of the atmospheric neutrino fluxes at high energy. The two main sources of uncertainty are the primary spectrum and the limited knowledge of hadronic interactions. The biggest uncertainty at high energy is the level of charm production. In a paper at this conference , work on a post-LHC model of the event generator SIBYLL was presented. The new version includes production of charm and will therefore give further insight into the level of prompt neutrinos. Charm is introduced in SIBYLL with a non-perturbative component, tuned to results of fixed target experiments, for example Refs. [37, 38] and a perturbative QCD component tuned to LHC results, for example from ALICE  and LHCb .
In Section 4 of this paper, we describe a method that can be used to explore in a systematic way the implications of uncertainties in hadronic interactions and primary spectrum for the fluxes of atmospheric neutrinos in the TeV-PeV energy range.
I thank Dr. M. Honda for providing details of the spectrum weighted moments for the hadronic interaction model used to calculate the atmospheric neutrino fluxes in Ref. . This research is supported in part by grants from the U.S. National Science Foundation, NSF-PHY-1205809 and from the U.S. Department of Energy, 12ER41808.
The original version of this paper contains an important typographical error in Eq. 11. (The pre-factor on the r.h.s. was incorrectly written as rather than .) The corresponding code used to produce the energy-dependent Z-factors was correct, so there are no changes in the results. In particular, the Z-factors and neutrino spectra in Fig. 3 are correct and unchanged. I am grateful to Felipe Campos Penha for finding this error.
-  IceCube Collaboration*, “Evidence for High-Energy Extraterrestrial Neutrinos at the IceCube Detector,” Science 342, 1242856 (2013) [arXiv:1311.5238 [astro-ph.HE]].
-  IceCube Collaboration*, “Observation of High-Energy Astrophysical Neutrinos in Three Years of IceCube Data,” Phys. Rev. Lett. 113, 101101 (2014) [arXiv:1405.5303 [astro-ph.HE]].
-  “On High energy Neutrino Physics” M.A. Markov in Proc. 1960 Annual Int. Conf. on High Energy Physics at Rochester, (ed. E.C.G. Sudarshan, J.H. Tinlot & A.C. Melissinos) pp. 578-581 (1060).
-  K. Greisen, “Cosmic ray showers,” Ann. Rev. Nucl. Part. Sci. 10, 63 (1960).
-  F. Reines, “Neutrino interactions,” Ann. Rev. Nucl. Part. Sci. 10, 1 (1960).
-  “Neutrino Production in the Atmosphere” G.T. Zatsepin & V.A. Kuz’min, Soviet Physics JETP 14 (1961) 1294.
-  L. V. Volkova, “Energy Spectra and Angular Distributions of Atmospheric Neutrinos,” Sov. J. Nucl. Phys. 31, 784 (1980) [Yad. Fiz. 31, 1510 (1980)].
-  T. K. Gaisser, “Cosmic rays and particle physics,” Cambridge, UK: Univ. Pr. (1990) 279 p
-  P. Lipari, “Lepton spectra in the earth’s atmosphere,” Astropart. Phys. 1, 195 (1993).
-  P. Gondolo, G. Ingelman and M. Thunman, “Charm production and high-energy atmospheric muon and neutrino fluxes,” Astropart. Phys. 5, 309 (1996) [hep-ph/9505417].
-  T. K. Gaisser and S. R. Klein, “A new contribution to the conventional atmospheric neutrino flux,” arXiv:1409.4924 [astro-ph.HE].
-  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. D 75, 043005 (2007) [astro-ph/0611201].
-  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. D 75, 043006 (2007) [astro-ph/0611418].
-  P. Achard et al. [L3 Collaboration], “Measurement of the atmospheric muon spectrum from 20-GeV to 3000-GeV,” Phys. Lett. B 598, 15 (2004) [hep-ex/0408114].
-  P. Adamson et al. [MINOS Collaboration], “Measurement of the atmospheric muon charge ratio at TeV energies with MINOS,” Phys. Rev. D 76, 052003 (2007) [arXiv:0705.3815 [hep-ex]].
-  OPERA Collaboration, N. Agafonova et al., “Measurement of the TeV atmospheric muon charge ratio with the complete OPERA data set,” Eur. Phys. J. C 74, 2933 (2014).
-  T. K. Gaisser, “Spectrum of cosmic-ray nucleons, kaon production, and the atmospheric muon charge ratio,” Astropart. Phys. 35, 801 (2012).
-  G. D. Barr, T. K. Gaisser, P. Lipari, S. Robbins and T. Stanev, “A Three - dimensional calculation of atmospheric neutrinos,” Phys. Rev. D 70, 023006 (2004)
-  Anatoli Fedynitch, presentation at ISVHECRI 2014.
-  V.A. Naumov & T.S. Sinegovskaya, “Simple Method for Solving Transport Equations Describing the Propagation of Cosmic-Ray Nucleons in the Atmosphere,” Physics of Atomic Nuclei 63, 1927-1935 (2000).
-  T. S. Sinegovskaya, A. D. Morozova and S. I. Sinegovsky, “High-energy neutrino fluxes and flavor ratio in the Earth atmosphere,” arXiv:1407.3591 [astro-ph.HE].
-  A. Fedynitch, J. Becker Tjus and P. Desiati, Phys. Rev. D 86, 114024 (2012) [arXiv:1206.6710 [astro-ph.HE]].
-  T. K. Gaisser, “Atmospheric leptons, the search for a prompt component,” EPJ Web of Conferences 52, 09004 (2013) arXiv:1303.1431 [hep-ph].
-  T. K. Gaisser, “Semianalytic approximations for production of atmospheric muons and neutrinos,” Astropart. Phys. 16, 285 (2002) [astro-ph/0104327].
-  B. Peters, ”Primary Cosmic Radiation and Extensive Air Showers,” Il Nuovo Cimento XXII, 800-819 (1961).
-  J. R. Hoerandel, “On the knee in the energy spectrum of cosmic rays,” Astropart. Phys. 19, 193 (2003) [astro-ph/0210453].
-  T. K. Gaisser, T. Stanev and S. Tilav, “Cosmic Ray Energy Spectrum from Measurements of Air Showers,” Front. Phys. China 8, 748 (2013) [arXiv:1303.3565 [astro-ph.HE]].
-  P. Lipari, “Establishing the astrophysical origin of a signal in a neutrino telescope,” arXiv:1308.2086 [astro-ph.HE].
-  T. K. Gaisser, T. Stanev, M. Honda and P. Lipari, “Primary spectrum to 1-TeV and beyond,” Proc. 27th ICRC (Hamburg, 2001).
-  S. Schoenert, T. K. Gaisser, E. Resconi and O. Schulz, “Vetoing atmospheric neutrinos in a high energy neutrino telescope,” Phys. Rev. D 79, 043009 (2009) [arXiv:0812.4308 [astro-ph]].
-  T. K. Gaisser, K. Jero, A. Karle and J. van Santen, “A generalized self-veto probability for atmospheric neutrinos,” Phys. Rev. D 90, 023009 (2014) [arXiv:1405.0525 [astro-ph.HE]].
-  J. W. Elbert “Multiple Muons Produced by Cosmic Ray Interactions,” in Proceedings of the DUMAND Summer Workshop (Scripps Institution of Oceanography, La Jolla), 1979, pp. 101-121.
-  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) [arXiv:0906.4113 [hep-ph]].
-  M. G. Aartsen et al. [IceCube Collaboration], “Atmospheric and Astrophysical Neutrinos above 1 TeV Interacting in IceCube,” arXiv:1410.1749 [astro-ph.HE].
-  R. Enberg, M. H. Reno and I. Sarcevic, “Prompt neutrino fluxes from atmospheric charm,” Phys. Rev. D 78, 043005 (2008) [arXiv:0806.0418 [hep-ph]].
-  Felix Riehn, “Charm production in SIBYLL,” Proc. ISVHECRI 2014 (this conference).
-  G. A. Alves et al. [E769 Collaboration], “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 (1996) [Erratum-ibid. 81, 1537 (1998)].
-  E. M. Aitala et al. [E791 Collaboration], “Asymmetries between the production of and mesons from 500-GeV/c - nucleon interactions as a function of and ,” Phys. Lett. B 371, 157 (1996) [hep-ex/9601001].
-  B. Abelev et al. [ALICE Collaboration], “Measurement of charm production at central rapidity in proton-proton collisions at TeV,” JHEP 1201, 128 (2012) [arXiv:1111.1553 [hep-ex]].
-  R. Aaij et al. [LHCb Collaboration], “Prompt charm production in pp collisions at sqrt(s)=7 TeV,” Nucl. Phys. B 871, 1 (2013) [arXiv:1302.2864 [hep-ex]].