Spectrum of cosmic-ray nucleons and the atmospheric muon charge ratio
aa r X i v : . [ a s t r o - ph . H E ] M a r Spectrum of cosmic-ray nucleons, kaon production, andthe atmospheric muon charge ratio
Thomas K. Gaisser
Bartol Research Institute and Dept. of Physics and AstronomyUniversity of Delaware, Newark, DE, USA
Abstract
Interpretation of measurements of the muon charge ratio in the TeV rangedepends on the spectra of protons and neutrons in the primary cosmic ra-diation and on the inclusive cross sections for production of π ± and K ± inthe atmosphere. Recent measurements of the spectra of cosmic-ray nucleiare used here to estimate separately the energy spectra of protons and neu-trons and hence to calculate the charge separated hadronic cascade in theatmosphere. From the corresponding production spectra of µ + and µ − the µ + /µ − ratio is calculated and compared to recent measurements. The com-parison leads to a determination of the relative contribution of kaons andpions. Implications for the spectra of ν µ and ¯ ν µ are discussed.
1. Introduction
The muon charge ratio in the TeV range has been measured by MI-NOS [1, 2] and more recently by OPERA [3]. Both analyses use an analyticapproximation as a framework for making an inference about the separatecontributions of the pion and kaon channels to the charge asymmetry. Inthis paper a more detailed derivation of the muon charge ratio is used forthe analysis. The muon charge ratio is expressed in terms of the spectrum-weighted moments for production of π ± and K ± by protons and neutronsin the primary cosmic radiation, following the analysis of Lipari [4]. Theanalysis here accounts for the special contribution of associated productionof charged, positive kaons. Email address: [email protected] (Thomas K. Gaisser)
Preprint submitted to Elsevier June 4, 2018 his analysis also accounts for the effect of the energy dependence of thecomposition of the primary cosmic-ray nuclei. Measurements from ATIC [5]and CREAM [6, 7] indicate that the spectra of helium and heavier nucleibecome somewhat harder than the spectrum of protons above several hundredGeV. This feature for helium was recently confirmed by PAMELA [8].Because muon neutrinos are produced together with muons in the pro-cesses π ± → µ ± + ν µ (¯ ν µ ) and K ± → µ ± + ν µ (¯ ν µ ) , (1)these results also apply to ν µ and ¯ ν µ . In the TeV range and above thecontribution of muon decay to the intensity of muon neutrinos is negligible.For reasons of kinematics, kaons are relatively more important for neutrinosat high energy than for muons. An additional goal of this paper is to drawattention to the implications of the muon results for atmospheric neutrinosin the TeV energy range and beyond.
2. Muon charge ratio
The excess of µ + in atmospheric muons can be traced to the excess ofprotons over neutrons in the primary cosmic-ray beam coupled with thesteepness of the cosmic-ray spectrum, which emphasizes the forward frag-mentation region in interactions of the incident cosmic-ray nucleons. Theclassic derivation of the muon charge ratio [9] considers muon productionprimarily through the channel p → π ± + anything. The atmospheric cascadeequation for the intensity of nucleons as a function of slant depth X in theatmosphere is solved separately for N = n + p and ∆ = p − n subject tothe appropriate boundary conditions. For the total intensity of nucleons asa function of slant depth X (g/cm ) φ N ( E ) = φ N (0) × exp( − X Λ N ) (2)where the nucleon attenuation length is Λ N = λ/ (1 − Z NN ) and λ is theinteraction length of nucleons in the atmosphere. The corresponding resultfor ∆( X ) = p ( X ) − n ( X ) is∆( X ) = δ φ N (0) × exp( − X Λ − ) , (3)where δ = p (0) − n (0) p (0) + n (0) and 1Λ − = 1 − Z pp + Z pn Z pp + Z pn N . (4)2he Z -factors (like Z NN = Z pp + Z pn ) are spectrum-weighted momentsof the inclusive cross sections for the corresponding hadronic process. Forexample, a particularly important moment for this paper is Z pK + = 1 σ Z x γ d σ ( x )d x d x (5)for the process p + air → K + + Λ + anything . (6)The normalized inclusive cross section is weighted by x γ where γ is the in-tegral spectral index for a power-law spectrum and x = E K /E p . Feynmanscaling is assumed in these approximate formulas, so the parameters mayvary slowly with energy, especially near threshold. However, the scalingapproximation is relatively good because the moment weights the forwardfragmentation region. The next step is to solve the coupled equations for the production ofcharged pions by nucleons separately for Π + ( X ) + Π − ( X ) and for ∆ π =Π + ( X ) − Π − ( X ). The solutions are then convolved with the probability perg/cm for decay to obtain the corresponding production spectra of muonsand neutrinos. The decay kinematic factors are1 − r γ +1 π ( γ + 1)(1 − r π ) and ǫ π cos θE µ − r γ +2 π ( γ + 2)(1 − r π ) (7)for muons and (1 − r π ) γ ( γ + 1) and ǫ π cos θE µ (1 − r π ) ( γ +1) ( γ + 2) (8)for neutrinos. In each of Eqs. 7 and 8 the first expression is a low-energylimit and the second a high energy limit, where low and high are with respectto the critical energy ǫ π . The ratio r π = m µ /m π = 0 . r K = 0 . ν µ + ¯ ν µ the expression is φ ν ( E ν ) = φ N ( E ν ) × ( A πν B πν cos( θ ) E ν /ǫ π + A Kν B Kν cos( θ ) E ν /ǫ K + A charm ν B charm ν cos( θ ) E ν /ǫ charm ) . (9)Here φ N ( E ν ) = dN/d ln( E ν ) is the primary spectrum of nucleons ( N ) evalu-ated at the energy of the neutrino. The three terms in brackets correspondto production from leptonic and semi-leptonic decays of pions, kaons andcharmed hadrons respectively. The term for prompt neutrinos from decay ofcharm has been included in Eq. 9 (see Ref. [10]) but will not be discussedfurther here.The numerator of each term of Eq. 9 has the form A iν = Z Ni × BR iν × Z iν − Z NN (10)with i = π ± , K, charm and BR iν is the branching ratio for i → ν . The first Z -factor in the numerator is the spectrum weighted moment of the crosssection for a nucleon (N) to produce a secondary hadron i from a targetnucleus in the atmosphere, defined as in Eq. 5. The second Z -factor is thecorresponding moment of the decay distribution for i → ν + X , which iswritten explicitly in Eq. 8. The second term in each denominator is the ratioof the low-energy to the high-energy form of the decay distribution [11]. Theforms for muons are the same, but the kinematic factors differ in a significantway (Eq. 7 instead of Eq. 8). Explicitly, for neutrinos B πν = γ + 2 γ + 1 ! (cid:18) − r π (cid:19) Λ π − Λ N Λ π ln(Λ π / Λ N ) ! (11)and for muons B πµ = γ + 2 γ + 1 ! − ( r π ) γ +1 − ( r π ) γ +2 ! Λ π − Λ N Λ π ln(Λ π / Λ N ) ! . (12)The forms for kaons are the same as functions of r K and Λ K .The separate solutions for π + → µ + + ν µ and π − → µ − + ¯ ν µ have theform φ π ( E µ ) ± = φ N ( E µ ) A πµ × . ± α π βδ )1 + B ± πµ cos( θ ) E µ /ǫ π , (13)4here B ± πµ = B πµ ± α π βδ ± c π α π βδ . Here β = 1 − Z pp − Z pn − Z pp + Z pn ≈ . β π = 1 − Z π + π + − Z π + π − − Z π + π + + Z π + π − ≈ . α π = Z pπ + − Z pπ − Z pπ + + Z pπ − ≈ . c π = 1 − Λ N / Λ π − β Λ N / ( β π Λ π ) " ln ( β π /β ) ln (Λ π / Λ N ) ≈ . . The numerical values are based on fixed target data in the energy range ofhundreds of GeV [11]. The factors B ± πµ differ by less than one per cent. Tothis accuracy, the charge ratio of muons can therefore be written in the form µ + µ − ≈ βδ α π − βδ α π = f π + − f π + , (14)where f π + = (1 + βδ α π ) / The situation becomes more complex when the contribution from kaonsis considered. In the first place, because the critical energies are significantlydifferent for pions and kaons, the two contributions have to be followed sepa-rately. In addition the charge ratio of muons from decay of charged kaons islarger than that from pion decay because the process of associated produc-tion in Eq. 6 has no analog for forward production of K − . Instead, associatedproduction by neutrons leads to Λ ¯ K .For the charge separated analysis of kaons it is useful to divide kaonproduction by nucleons into a part in which K + and K − are produced equallyby neutrons and by protons and another for associated production, which istreated separately. Then in the approximation that kaon production by pionsin the cascade is neglected, the spectrum of negative muons from decay of K − is φ K ( E µ ) − = Z NK − Z NK φ N ( E µ ) A NK B Kµ cos( θ ) E µ /ǫ K . (15)5here is an equal contribution of central production to positive kaons,but in addition there is the contribution from associated production. Thetotal contribution of the kaon channel to positive muons is φ K ( E µ ) + = φ N ( E µ ) A NK × (1 + α K βδ )1 + B + Kµ cos( θ ) E µ /ǫ K . (16)Here α K = Z pK + − Z pK − Z pK + + Z pK − and B + Kµ = B Kµ × βδ α K βδ α K (1 − ln ( β ) /ln (Λ K / Λ N )) . Combining the expressions for µ + and µ − from pions (Eq. 13) and fromkaons (Eqs. 15 and 16), the muon charge ratio is µ + µ − = " f π + B πµ cos( θ ) E µ /ǫ π + (1 + α K βδ ) A Kµ /A πµ B + Kµ cos( θ ) E µ /ǫ K × " (1 − f π + )1 + B πµ cos( θ ) E µ /ǫ π + ( Z NK − /Z NK ) A Kµ /A πµ B Kµ cos( θ ) E µ /ǫ K − . (17)For the pion contribution, isospin symmetry allows the pion terms in thenumerator and denominator to be expressed in terms of f + π as defined afterEq. 14 above. The kaon contribution does not have the same symmetry.Numerically, however, the differences are at the level of a few per cent, asdiscussed in the results section.
3. Primary spectrum of nucleons
What is relevant for calculating the inclusive spectrum of leptons in theatmosphere is the spectrum of nucleons per GeV/nucleon. This is because, toa good approximation, the production of pions and kaons occurs at the levelof collisions between individual nucleons in the colliding nuclei. To obtainthe composition from which the spectrum of nucleons can be derived we usethe measurements of CREAM [6, 7], grouping their measurements into theconventional five groups of nuclei, H, He, CNO, Mg-Si and Mn-Fe.Direct measurements of primary nuclei extend only to ∼
100 TeV totalenergy. Because we want to calculate spectra of muons and neutrinos up to6 E . d N / d E ( m - s r - s - G e V . ) E total (GeV)p HeCNOMgSi FeGrigorovAkenoMSUKASCADEHEGRACasaMiaTibet-SIBYLLKASCADE-GrandeAGASAHiRes1&2Auger2009Allparticle fit E . d N / d E ( m - s r - s - G e V . ) E N (GeV/nucleon)All nucleon -2.7Polygonato Figure 1: Left: three-population model of the cosmic-ray spectrum from Eq. 21 comparedto data [12–22]. The extra-galactic population in this model has a mixed composition.Right: Corresponding fluxes of nucleons compared to an E − . differential spectrum ofnucleons and to the all nucleon flux implied by the Polygonato model (galactic componentonly) [25]. a PeV, we need to extrapolate the direct measurements to high energy ina manner that is consistent with measurements of the all-particle spectrumby air shower experiments in the knee region (several PeV) and beyond, asillustrated in the left panel of Fig. 1. To do this we adopt the proposal ofHillas [23] to assume three populations of cosmic rays. The first populationcan be associated with acceleration by supernova remnants, with the kneesignaling the cutoff of this population. The second population is a higher-energy galactic component of unknown origin (“Component B”), while thehighest energy population is assumed to be of extra-galactic origin.Following Peters [24] we assume throughout that the knee and other fea-tures of the primary spectrum depend on magnetic rigidity, R = pcZe , (18)where Ze is the charge of a nucleus of total energy E tot = pc . The motivationis that both acceleration and propagation in models that involve collision-less diffusion in magnetized plasmas depend only on rigidity. The rigidity7 c γ p He CNO Mg-Si Fe γ for Pop. 1 —- 1.66 1.58 1.63 1.67 1.63Population 1: 4 PV see line 1 7860 3550 2200 1430 2120Pop. 2: 30 PV 1.4 20 20 13.4 13.4 13.4Pop. 3 (mixed): 2 EV 1.4 1.7 1.7 1.14 1.14 1.14” (proton only): 60 EV 1.6 200. 0 0 0 0 Table 1: Cutoffs, integral spectral indices and normalizations constants a i,j for Eq. 21. determines the gyroradius of a particle in a given magnetic field B accordingto r L = R / B . (19)Peters pointed out that if there is a characteristic rigidity, R c above whicha particular acceleration process reaches a limit (for example because thegyroradius is larger that the accelerator), then the feature will show up intotal energy first for protons, then for helium and so forth for heavier nucleiaccording to E ctot = A × E N,c = Ze × R c . (20)Here E N is energy per nucleon, A is atomic mass and Ze the nuclear charge.The first evidence for such a Peters cycle associated with the knee of thecosmic-ray spectrum comes from the unfolding analysis of measurements ofthe ratio of low-energy muons to electrons at the sea level with the KAS-CADE detector [15].In what follows we assume that each of the three components ( j ) containsall five groups of nuclei and cuts off exponentially at a characteristic rigidity R c,j . Thus the all-particle spectrum is given by φ i ( E ) = Σ j =1 a i,j E − γ i,j × exp " − EZ i R c,j . (21)The spectral indices for each group and the normalizations are given explicitlyin Table 1. The parameters for Population 1 are from Refs. [6, 7], which weassume can be extrapolated to a rigidity of 4 PV to describe the knee. InEq. 21 φ i is d N/ dln E and γ i is the integral spectral index. The subscript i = 1 , F r a c t i on o f p r i m a r y H y d r ogen E N (GeV/nucleon) Figure 2: Solid line: charge ratio parameter δ for the model with parameters of Table 1.Dashed line: same for Polygonato model [25]. The composite spectrum corresponding to Eq. 21 and Table 1 is super-imposed on a collection of data in the left panel of Fig. 1. No effects ofpropagation in the galaxy or through the microwave background have beenincluded in this phenomenological model. For the two galactic components,however, a consistent interpretation could be obtained with source spectra γ ∗ ∼ . γ ∗ ∼ .
07 for population 2 together with anenergy dependent diffusion coefficient D ∼ E δ with δ = 0 .
33 for both com-ponents to give local spectra of γ = γ ∗ + δ of ∼ .
63 and ∼ . φ i,N ( E N ) = A × φ i ( A E N ) (22)for each component and then summing over all five components. The nucleonspectrum is shown in the right panel of Fig. 1.The energy-dependent charge ratio δ ( E N ) needed to calculate the muoncharge ratio follows from Eq. 22 and Table 1. To a good approximation, it isgiven by the fraction of free hydrogen in the spectrum of nucleons, as shownin Fig. 2. The fraction decreases slowly from its low energy value of 0 .
76 at10 GeV/nucleon [26] to a minimum of 0 .
63 at 300 TeV and then increases9omewhat at the knee. Note that, because of the relation among E tot , E N and R c in Eq. 20, the steepening at the knee occurs for nuclei at Z/A ≈ theenergy per nucleon as compared to protons. Hence the free proton fractionrises again at the knee.Also shown for comparison in Fig. 2 by the broken line is the δ param-eter for the rigidity-dependent version of the Polygonato model, which hasa common change of slope ∆ γ = 1 . E N > GeV/nucleon doesnot affect the charge ratio, which is measured only for E µ < GeV. It istherefore possible to consider the difference between the two versions of δ inFig. 2 as a systematic effect of the primary composition.
4. Comparison with data
We now wish to compare the calculation of Eq. 17 to various sets of datausing the energy-dependent primary spectrum of nucleons (Eq. 22) with pa-rameters from Table 1. There are two problems in doing so. First, expressionsfor the intensity of protons and neutrons from Eqs. 2 and 3 and the subse-quent equations are valid under the assumption of a power-law spectrumwith an energy independent value of δ . The assumption of a power lawwith integral spectral index of − . δ ( E N ) = δ (10 × E µ ).The other problem is that the data are obtained over a large range ofzenith angles, and the charge ratio also depends on angle. The first MINOSpublication [1] gives µ + /µ − as a function of the energy of the muon at thesurface. These data are shown in Fig. 3 along with older high energy datafrom the Park City Mine in Utah [27] and data at lower energy from L3 [28]and CMS [29]. The figure shows three calculations of the muon charge ratioin the vertical direction that follow from Eq. 17. The highest curve assumesa constant composition fixed at its low energy value, δ = 0 .
76 [26]. The10 µ + / µ - r a t i o E µ (TeV)CMSL3-CMINOSenergy independentcomposition(E)plus reduced K/ π Figure 3: Muon charge ratio compared to data of CMS [29] and L3-C [28] below 1 TeVand to MINOS [1] at higher energy. The L3-C data plotted here are averaged over0 . ≤ cos( θ ) ≤ . middle curve is the result assuming the energy-dependent composition pa-rameter δ ( E N ) that corresponds to the parameterization of Table 1 (solidline in Fig. 2), which is still higher than the data. Both the higher lines as-sume the nominal values of the spectrum weighted moments from Ref. [11].The lowest curve is obtained by reducing the level of associated production,by changing Z pK + from its nominal value of 0 . . θ ) ≈ . θ ) ≈ . Z pK + , which reflects p → K + . µ + / µ - r a t i o E µ cos θ (TeV)CMSOPERAMINOSUtahenergy independentreduced K/ π low-Ehigh-E Figure 4: Muon charge ratio compared to data of CMS [29], OPERA [3] and MINOS [2].Measurements from the near detector of MINOS [30] and from Park City [27].
The dependence of µ + /µ − on zenith angle enters Eq. 17 in the form E µ cos( θ ). For this reason the muon charge ratio is often presented as afunction of this combination. The data of OPERA are presented only interms of the product E µ cos( θ ). Because of the complex overburden at GranSasso, there is no simple relation between zenith angle and energy. TheMINOS data are also presented in this form in Ref. [2], but the mean energyfor each value of E µ cos( θ ) is not given. For a primary cosmic-ray compositionthat has no energy dependence, Eq. 17 depends only on E µ cos( θ ). Theeffect of the energy-dependence of composition can threrfore be assessed bycomparing the calculation for various fixed values of δ to the data, which isdone in Fig. 4. 12he upper curve in Fig. 4 is the same as the corresponding curve in Fig. 3,plotted for a constant composition with δ = 0 .
76, its value at 10 GeV. Theparameter δ decreases from 0.71 at 100 GeV/nucleon to 0.68 at a TeV, andfrom 0.64 at 10 TeV/nucleon to less than 0.62 at 100 TeV. The full curvethrough the data in Fig. 4 is evaluated for δ = 0 .
5. Summary
The muon charge ratio is sensitive both to the proton excess in the spec-trum of primary cosmic-ray nucleons and to the value of Z pK + . Using recentdata on primary composition, we find a proton excess that decreases steadilyfrom 10 GeV/nucleon to 500 TeV. This portion of the cosmic-ray spectrumproduces muons from a few GeV to well over 10 TeV. Assuming associatedproduction (Eq. 6) to be the major uncertainty, a level of associated produc-tion in the range Z pK + = 0 . ± . Z pK + = 0 . R K/π = Z pK + + Z pK − Z pπ + + Z pπ − = 0 . . .
046 + 0 .
033 = 0 . . (23)It is interesting that analyses of seasonal variations of TeV muons byMINOS [31] and IceCube [32] also suggest a somewhat lower value of R Kπ than its nominal value of 0.149. On the other hand, the value in Eq. 23still represents a significant contribution from the K + decay channel. If theenergy-dependent composition of the Polygonato model is used instead, agood fit is obtained with Z pK + = 0 . δ in the relevant energy range (0 .
68 as compared to 0 . R K/π = 0 . µ + µ − = " f π + B πµ cos( θ ) E µ /ǫ π + f K + A Kµ /A πµ B Kµ cos( θ ) E µ /ǫ K × " (1 − f π + )1 + B πµ cos( θ ) E µ /ǫ π + (1 − f K + ) A Kµ /A πµ B Kµ cos( θ ) E µ /ǫ K − (24)with A Kµ /A πµ = 0 . f + π = 0 .
55 and f + K = 0 .
67. For E µ ∼ TeV, δ ≈ .
64 forprimary energy per nucleon of 10 TeV. Thus f + π = (1 + α π βδ ) = 0 .
55 inagreement with the MINOS analysis. A precise comparison with the MINOSvalue for f + K is not possible for the reasons explained after Eq. 17. However,numerical differences are at the level of a few per cent. For example, B + Kµ ≈ . × B Kµ . From Eqs. 15 and 16, the value of f K + = φ K ( E µ ) + / ( φ K ( E µ ) + + φ K ( E µ ) − ) ≈ .
69 in the TeV region for δ ≈ .
64. However, if the expressionfor the kaon contribution in Eq. 17 is expressed in terms of f + K there is anadditional multiplicative factor less than unity. Thus, although the forms aredifferent, the fits are much the same.The role of kaons is relatively more important for neutrinos than formuons. Because the muon mass is close to that of the pion, the muon car-ries most of the energy of the decaying pion. Kaons split the energy almostequally on average between the µ and the ν µ . The steep spectrum enhancesthe effect so that kaons are the dominant source of muon neutrinos above afew hundred GeV. Forward production of K + is therefore particularly im-portant. The effect is illustrated in Fig. 5, in which the ratio ν µ / ¯ ν µ is plottedfor the same sequence of assumptions as in the plot of the muon charge ratio(Fig. 3). The implications of the muon charge ratio for neutrinos will be thesubject of a separate paper. Acknowledgments : I am grateful to Anne Schukraft and Teresa Montarulifor comments on an early version of this paper. I am grateful for helpfulcomments from an external reviewer and to Jeffrey de Jong and Simone Biagifor information about MINOS and OPERA. This research is supported inpart by the U.S. Department of Energy under DE-FG02-91ER40626.14 ν µ / an t i - ν µ r a t i o E ν (TeV)energy independentcomposition(E)plus reduced K/ π Polygonato K/ π Figure 5: ratio of ν µ / ¯ ν µ calculated with the same parameters as for the charge ration ofmuons. The dashed line shows the results with parameters of the Polygonato model [25]. References [1] P. Adamson et al. , (MINOS Collaboration) Phys. Rev. D , 052003(2007).[2] P.A. Schreiner, J. Reichenbacher, M.C. Goodman, Astropart. Phys. ,61 (2009).[3] N. Agafonova et al., (OPERA Collaboration), Eur. Phys. J. C , 25(2010).[4] Paolo Lipari, Astropart. Phys. 1 (1993) 195-227.[5] A.D. Panov, et al. , Bulletin of the Russian Academy of Sciences:Physics, 73 (2009)564-567.[6] H.S. Ahn et al. , Ap.J. 714 (2010) L89-L93.157] H.S. Ahn, et al. , Ap.J. 707 (2009) 593-603.[8] O. Adriani et al. (PAMELA) Science 332 (2011) 69-72.[9] W.R. Frazer et al., Phys. Rev. D 5 (1972) 1653.[10] P. Desiati & T.K. Gaisser, Phys. Rev. Letters 105 (2010) 121102.[11] Cosmic Rays and Particle Physics , T.K. Gaisser, (Cambridge UniversityPress, 1990).[12] N.L. Grigorov et al., Yad. Fiz. 11 (1970) 1058.[13] M. Nagano et al. (Akeno), J. Phys. G 10 (1984) 1295-1310 and 18 (1992)423-442.[14] Yu. A. Fomin, et al. (MSU), J. Phys. G: Nucl. Part. Phys. 22 (1996)1839.[15] T. Antoni et al. , Astropart. Phys. 24 (2005) 1.[16] F. Arqueros, et al. (HEGRA), Astronomy & Astrophysics 359 (2000)682.[17] M.A.K. Glasmacher et al. (CASA MIA) Astropart. Phys. 10 (1999) 291.[18] M. Amenomori, et al. (Tibet), Astrophys. J. 678 (2008) 1165[19] W.D. Apel et al. (KASCADE-Grande), Phys. Rev. Letters 107 (2011)171104.[20] M. Takeda et al. (AGASA) Astropart. Phys. 19 (2003) 447.[21] R.U. Abbasi et al. (HiRes) Astropart. Phys. 32 (2009) 53.[22] Pierre Auger Collaboration, Phys. Rev. Letters 104 (2010) 091101.[23] A.M. Hillas, arXiv:astro-ph/0607109.[24] B. Peters, Il Nuovo Cimento XXII (1961) 800-819.[25] J¨org R. H¨orandel, Astropart. Phys. 19 (2003) 193.1626] T.K. Gaisser & Todor Stanev, “Cosmic Rays” in
Reviews of ParticlePhysics , K. Nakamura et al., J. Phys. G 37 (2010) 075021.[27] G.K. Ashley, J.W. Keuffel & M.O. Larson, Phys. Rev. D , 20 (1975).[28] P. Achard et al. (L3 Collaboration) Phys. Lett. B 598 (2004) 15.[29] CMS Collaboration, Physics Letters B ,83 (2010).[30] P. Adamson et al. , (MINOS Collaboration) Phys. Rev. D83