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 $\nu_\mu$ and $\nu_e$ are discussed, along with the implications of the muon charge ratio for the $\nu_\mu/\bar{\nu}_\mu$ 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.


Introduction
Atmospheric leptons are of current interest in two contexts, as a beam for the study of neutrino oscillations and the mass hierarchy (energy range 100 MeV→TeV) and as the background in the search for high energy neutrinos of astrophysical origin (energy range 100 GeV→PeV). 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.

History
It is well known that the idea of using Cherenkov light in a large volume of water to detect interactions of highenergy neutrinos was suggested in 1960 by Markov [3], by Greisen [4] and by Reines [5]. 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 a e-mail: gaisser@bartol.udel.edu 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 3 of water.
Perhaps less well known is the early (1961) paper [6] 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 sec θ 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 [7], which remains a standard reference for fluxes of atmospheric neutrinos.

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 N i (E i , X)dE i is the flux of particles of type i at slant depth X in the atmosphere with energies in the interval E to E + dE [8,9]. The probability for a particle of type j to interact in dX is dX/λ j (E j ), where λ j is the corresponding interaction length (in g/cm 2 ). Similarly, dX/d j (E j ) is the probability that a particle of type j decays in dX. The where dn i is the number of particles of type i produced on average in the energy bin dE i around E i per collision of an incident particle of type j. Depending on the boundary condition, the solution of Eq. 1 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 N N (X = 0) = φ N (E N ) and N i (X = 0) = 0 for i N. Here φ N (E N ) is the spectrum of nucleons at the top of the atmosphere. Correspondingly, the solutions N i (E i , X, θ) give the inclusive rates of particles of type i per unit area at depth X independently of whether or not there are particles nearby. In contrast, for the air shower boundary condition (N A (E A , X = 0) = δ(E = E 0 )) the solution in principle contains all the correlations among the particles in the shower.

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. 1 for the flux of ν µ +ν µ is where the sum is over the contributions of charged pions, charged kaons and charmed hadrons. At low energy (< 100 GeV/ cos * θ) 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 θ > 60 • [7]. The quantity i / cos * θ is the energy for meson i above which re-interaction in the atmosphere becomes more likely than decay. At energies below the critical energy ( i ) 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 E > i 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 charm ∼ 10 7 GeV, so the small fraction of prompt neutrinos will be isotropic below this energy. The numerator of Eq. 3 contains the spectrumweighted moments for production of secondary particles in particle collisions in the atmosphere. For example, for nucleon→ K + → ν µ where and Z NN has a similar definition in terms of the cross section for a nucleon of energy E to produce a secondary nucleon of energy E. Z K + ν is the spectrum weighted moment of the decay distribution. Eq. 5 is a general definition given in Ref. [10] 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 1 + γ, constant interaction cross sections and particle production that depends only on the ratio x = E/E , Eq. 5 simplifies to the more familiar, energy-independent form,

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 ν e and ν µ , which means that the fraction of prompt ν e is significantly larger than the fraction of prompt ν µ . The lower normalization of ν e from kaons also means that the contribution from muon decay is relatively more important to higher energy and lower zenith angle than for ν µ . A  Figure 2. Ratio of ν µ /ν µ inferred from the the µ + /µ − ratio of OPERA [16].
detailed calculation of electron neutrinos at high energy, including the contribution from the rare three-body decay of K S , is given in Ref. [11]. Because of its short lifetime, the critical energy for K S is 120 TeV. As a consequence, the spectrum of neutrinos from this channel is harder by one power of energy than that from K ± and K L . Above 100 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.

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 [12] are used to tune the model for hadron production for one of the standard neutrino flux calculations [13]. 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 K . The signature of the kaon channel nevertheless shows up in an important way in the charge ratio of muons, which increases from ≈ 1.28 around 100 GeV [14] to ≈ 1.4 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 ≈ 115 GeV compared to ≈ 850 GeV for kaons. There is an increase in the muon charge ratio in this energy range because of the importance of associated production, p → K + Λ, 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 [17]. 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. [17] to fit the proton excess in the primary spectrum and the Z-factor, the OPERA group found δ 0 = (p 0 − n 0 )/(p 0 + n 0 ) = 0.61 ± 0.02 and Z pK + = 0.0086±0.0004 at a mean muon energy of 2 TeV [16]. This result has implications for the ratio ν µ /ν µ in the same energy region. With the value of Z pK + from OPERA, the expected ratio for muon neutrinos increases from ν µ /ν µ ≈ 1.5 at low energy to ≈ 2.3 above a TeV, as shown in Fig. 2.

Analytic and numerical calculations of atmospheric neutrinos
A formal expression for the flux of neutrinos is where φ A is the primary flux of nuclei with total energy E A and Y ν (A, E ν , E A , θ) 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. [18], 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 log(E A ). 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 E ν = 10 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 [19] and the iterative scheme [20] of Sinegovskaya et al. [21]. 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.

Energy-dependent Z-factors
In their paper on neutrinos and muons from charm decay, Gondolo, Ingelman & Thunman [10] showed that the simple Eq. 3 valid for a power-law primary spectrum and scaling can be adapted to include energy dependence. With the EPJ Web of Conferences TG ! Honda ! TG e Honda e Figure 3. Left: Energy-dependent Z-factors from Ref. [24] (solid lines) compared to those of Ref. [13] (broken lines). Right: Spectrum of ν µ +ν µ and ν e +ν e for two models of particle production calculated with the primary spectrum of Ref. [17], which includes the knee. The crosses are the tabulated neutrino fluxes averaged over the sky from Ref. [13] (see text for discussion).
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. [22] and Ref. [23] 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 energydependent 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. [24] is sufficient, provided the parameters are tabulated at a sufficiently fine grid of energies. The production spectrum for with x = E i /E j and the two energy-dependent parameters c i, j (E j ) and p i, j (E j ). 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 and Having determined the parameters for particle production, The next step is to evaluate the spectrum-weighted mo-ments at each energy for an arbitrary spectrum using (For simplicity, the factor λ(E)/λ(E/x) in Eq. 5 has been approximated as unity here because the interaction cross section varies slowly over the range of x near 1 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. [24], which correspond to the Z-factors of Ref. [8] for γ = 1.7. In this case, the only energy dependence in the Z i, j (E) is from the knee of the spectrum, as in Ref. [23]. For comparison, the more realistic case Z i, j (E) corresponding to the hadronic interaction model of Ref. [13] are used. The energydependent Z-factors (Eq. 12) are shown in Fig. 3 (left) evaluated with the H3A primary spectrum of Ref. [17]. The corresponding fluxes of ν µ +ν µ and ν e +ν e 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. [13] 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.

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 nucleonnucleon collisions. We are interested in an energy range International Symposium on Very High Energy Cosmic Ray Interactions 2014  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 [25]. 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 [26] 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 E A = ZeR cut , where R cut , the cutoff rigidity, is tuned to fit the knee of the cosmic-ray spectrum. In Ref. [17], 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 R 1 = 4 PV (for the knee) and R 2 = 30 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. [27]. 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 [28] 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. [13] modified to include the knee as in H3a. This parameterization is the high helium version of a pa-SUPPL. FIG. 6. Comparison of zenith distributions for atmospheric neutrino flux with charm saturating previous limits [9] before (dashed purple line) and after (solid purple line) removal of events accompanied into the detector by muons from the neutrinos' parent air shower [27,28].  rameterization and extrapolation of direct measurements from Ref. [29]. The ∼30% range in the knee region in Fig. 4 translates into a comparable uncertainty for the atmospheric neutrino flux above ∼30 TeV.

Neutrino yields and neutrino self-veto
The High Energy Starting Event (HESE) analysis in Ice-Cube [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 [30]. 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 [31] consists of evaluating an expression similar to Eq. 7 but with an extra factor in the integrand: (13) EPJ Web of Conferences The exponential factor is the Poisson probability that a shower generated by a primary nucleus with energy E A has no muons with E µ > E min (θ), 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 [32] was obtained from simulations with Sibyll 2.1 [33] 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, p 3 = 1. In the expression for prompt neutrinos the factor 1/E i cos * θ is omitted. The differential yields are Then the passing rate for atmospheric neutrinos of energy E ν 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. [34]. 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 [2] 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 E ν > 60 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 cos θ = − sin δ. The veto suppresses the downward atmospheric background by 50% or more for zenith angles θ < 70 • , a shape completely different from the data.

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. [1]. Expected rates for flavor i 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 60 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 1 : 1 : 1 at Earth and for an astrophysical spectrum the expected numbers of events in 988 days are ν µ : ν τ : ν e ≈ 5 : 7 : 11. For comparison, the total number of events, including backgrounds in the 3-year analysis [2] is 37. The astrophysical spectrum in Eq. 17 is a fit with a differential spectral index of −2.3 from Ref. [2].
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 100 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. [31]. 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. [35] assuming a primary spectrum with the knee as in Ref. [17]. 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 ν e 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 ν e and most ν τ will be classified as cascades. At present the observations are consistent with a 1 : 1 : 1 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.

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 [36], 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 nonperturbative 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 [39] and LHCb [40].
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.