Calculation of atmospheric neutrino flux using the interaction model calibrated with atmospheric muon data M. Honda∗ and T. Kajita† Institute for Cosmic Ray Research, the University of Tokyo, 5-1-5 Kashiwa-no-ha, Kashiwa, Chiba 277-8582, Japan 7 K. Kasahara‡ 0 0 Shibaura Institute of Technology, 307 Fukasaku, 2 n Minuma-ku, Saitama 330-8570, Japan a J 0 S. Midorikawa§ 1 Faculty of Software and Information Technology, 3 v Aomori University, 2-3-1 Kobata, Aomori 030-0943, Japan. 8 1 4 1 T. Sanuki¶ 1 6 International Center for Elementary Particle Physics, 0 / h the University of Tokyo, 7-3-1 Hongo, p - Bunkyo-ku, Tokyo 113-0033, Japan o r t (Dated: February 5, 2008) s a : v i X r a 1 Abstract Using the “modified DPMJET-III” model explained in the previous paper [1], we calculate the atmospheric neutrino flux. The calculation scheme is almost the same as HKKM04 [2], but the usage of the “virtual detector” is improved to reduce the error due to it. Then we study the un- certainty of the calculated atmospheric neutrino flux summarizing the uncertainties of individual components of the simulation. The uncertainty of K-production in the interaction model is esti- mated using other interaction models: FLUKA’97 and Fritiof 7.02, and modifying them so that they also reproduce the atmospheric muon flux data correctly. The uncertainties of the flux ratio and zenith angle dependence of the atmospheric neutrino flux are also studied. PACS numbers: 95.85.Ry, 13.85.Tp, 14.60.Pq ∗[email protected]; http://icrr.u-tokyo.ac.jp/∼mhonda †[email protected] ‡[email protected] §[email protected] ¶[email protected] 2 I. INTRODUCTION In the previous paper [1] (hereafter Paper I), we have studied the interaction model (DPMJET-III[3]) employed intheHKKM04atmospheric neutrino fluxcalculation[2], using the atmospheric muon flux data observed by precision measurements [4, 5, 6]. In the study, we foundthat thecalculatedandobserved muonfluxes didnot agree, especially for momenta above 30 GeV/c. We modified the interaction model to improve the agreement between the calculated and observed atmospheric muon fluxes. We call this modified interaction model the modified DPMJET-III in this paper. Note, the modification is actually applied to the “inclusive DPMJET-III” in a phenomenological way based on the quark-parton model. In this paper, we calculate the atmospheric neutrino flux with the modified DPMJET-III (Sec. II). The calculation scheme and the physical input data are basically the same as the HKKM04. However, we update the geomagnetic model from IGRF2000 to IGRF2005 [7], and improve the use of “virtual detector” in the 3-dimensional calculation to reduce the error due to it [8]. There are other hadronic interaction models which are used in the detector simulations of high energy experiments, such as FLUKA’97 [9] and Fritiof 7.02 [10]. We calculate the atmospheric neutrino fluxes with these interaction models applying the modification, so that they also reproduce the atmospheric muon flux observed by the precision measurements (Sec. III). Note, to reproduce the observed muon flux, modifying the primary flux model might be alternative solution. We calculated the atmospheric neutrino flux, changing the spectralindexofprimarycosmic rayprotonsfrom−2.71to−2.66above100GeV,whichalso reproduces the observed muon flux in µ++µ− sum correctly with the original DPMJET-III. Those calculations give almost the same atmospheric neutrino flux in the energy region below 100 GeV, where π’s are the main source of atmospheric neutrinos. With the modifica- tions based on the atmospheric muon data, the π production is almost the same in all three calculations. However the K’s are not related to the atmospheric muons below 1 TeV/c, and above this momenta, almost no muon flux is available from the precision measurements. There remain sizable differences in the K production, resulting in differences in atmospheric neutrino fluxes at higher energies. In Sec. IV, we estimate the uncertainties in our calculations. As the uncertainties of the predicted atmospheric neutrino flux is crucial for the study of neutrino oscillations, the 3 study of them is important [11]. Since our calculation reproduces the observed atmospheric muon flux data, the uncertainty due to that in the π production could be estimated from the experimental error and the residual of the reconstruction of atmospheric muon data. The uncertainty of K production is estimated from the variation of the atmospheric neutrino flux at higher energies calculated in the modified calculation schemes. The total uncertainty is estimated by summarizing individual uncertainties. Note, the stabilities of the (ν + µ ν¯ )/(ν + ν¯ ) ratio and the zenith angle dependence of the atmospheric neutrino flux are µ e e especially important in the study of neutrino oscillations. They are also studied with the uncertainty of the flux value. II. CALCULATION OF ATMOSPHERIC NEUTRINO FLUX Inthissection, we describe thecalculation oftheatmospheric neutrino flux with themod- ified DPMJET-III in detail. The calculation scheme is basically the same as HKKM04 [2]. For the primary flux model, we take the same primary flux model as HKKM04, based on AMS [12, 13] and BESS [4, 14] data, with a spectral index of −2.71 above 100 GeV (see also Refs. [15, 16]). For the model of the atmosphere, we used the US-standard ’76 [17], as the error due to the atmospheric density model is sufficiently small for the calculation of atmospheric neutrino flux [1]. Note, however, we update the geomagnetic field model from IGRF2000 to IGRF2005 [7], and improve of the usage of the “virtual detector” as explained in this section. Note, there are a considerable number of 3-dimensional calculations of the atmospheric neutrino flux [8, 18, 19, 20, 21, 22, 23]. However, all those calculations suffer from the small statistics athigher neutrino energies (&10 GeV),due tothe inefficiency ofthe 3-dimensional calculation scheme. In this paper, we calculate the atmospheric neutrino flux averaged over all azimuthal angles, combining 3-dimensional and 1-dimensional calculations. In HKKM04, it is shown that the atmospheric neutrino flux calculated with a 1-dimensional scheme agrees with that calculated with the 3-dimensional scheme above a few GeV, averaged over all azimuthal angles, although more than a few % azimuth angle dependence remains in the atmospheric neutrino flux at 10 GeV due to muon curvature in the geomagnetic field. For the 3-dimensional calculation, we assume the surface of the Earth is a sphere with R = 6378.180 km. We also assume three more spheres; the injection, simulation, and e 4 escape spheres. The radius of the injection sphere is taken as R = R + 100 km, the inj e simulation sphere as R = R +3000 km, and the escape sphere as R = 10×R . Note, a sim e esc e spheroid with an eccentricity of ∼1/298 is a better approximation for the Earth. However, the differece from the sphare approximation is small, and we estimate the errors of the atmospheric neutrino flux due to the sphere approximation is also small (. 1%). Cosmic rays are sampled on the injection sphere uniformly toward the inward direction, following the given primary cosmic ray spectra. Before they are fed to the simulation code for the propagation in air, they are tested to determine whether they pass the rigidity cutoff, i.e., the geomagnetic barrier. For a sampled cosmic ray, the ‘history’ is examined by solving the equation of motion in the negative time direction. When the cosmic ray reaches the escape sphere without touching the injection sphere again in the inverse direction of time, the cosmic ray can pass through the magnetic barrier following its trajectory in the normal direction of time. The propagation of cosmic rays is simulated in the space between the surface of Earth and the simulation sphere. When a particle enters the Earth, it loses its energy very quickly, and results in neutrinos with energy less than 100 MeV. Therefore, we discard such particles as soon as they enter the Earth, as most neutrino detectors which observe atmospheric neutrinos do not have sensitivity below 100 MeV. For secondary particles produced in the interaction of a cosmic ray and air-nucleus, there is the possibility that they go out from and re-enter in the atmosphere and create low energy neutrinos. Therefore, a simulation sphere which is toosmallmay miss such secondary particles. Ontheother hand, it isvery timeconsuming to followallparticles outto distances far from the Earth. In HKKM04, the simulation sphere with radius R = R +3000 km sim e was found to be large enough to suppress the error to well below the 1 % level. Note, neutrino detectors are very small compared with the size of the Earth, and are considered as the infinitesimal points on the surface of the Earth. We introduce a finite size “virtual detector” for each targert detector in the 3-dimensional calculation scheme. In HKKM04, the surface of the Earth within a circle around the target detector of radius θ = 10◦ (∼ 1000 km) is used as the virtual detector. When neutrinos pass throuth the d surface of the Earth inside of the circle (upward or downward), they are registered. We do not need the virtual detector in 1-dimensional calculation scheme, since it treats the propagation of the cosmic rays on a line which go through the neutrino detector. This is a 5 far less time consuming computation scheme than the 3-dimensional calculation scheme for the atmospheric neutrino flux, The finite size of the virtual detector introduces an error, since it averages the neutrino flux over positions where the geomagnetic conditions are different from the position of target detector [8]. To study the relation between the size and the error, we calculate the atmo- spheric neutrino flux with different size virtual detectors, θ = 10◦ (Φ (10◦)) and θ = 5◦ d ν d (Φ (5◦)). The fluxes Φ (10◦) and Φ (5◦) are compared in Fig. 1 for Kamioka from the ν ν ν HKKM04 calculation averaging over all azimuthal angles. We find a difference is seen for downward directions, and is almost constant for cosθ > 0 in ratio, where θ is the zenith z z angle of the arrival direction of the neutrinos. Therefore we expect the maximum error at E = 0.1 GeV is ∼ 5 % for downward directions averaging over azimuthal angles. The error ν due to the finite size virtual detector are smaller than those due to uncertainty of hadronic interaction model in HKKM04. All Direction Average At E ν = 0.1 GeV 0 )o 1.01 0 )o 1.02 1 1 5 )/φ (oν 1 νµ 5 )/φ (oν 1 φ (ν νµ φ (ν ν e 0.99 0.98 νµ ν e ν ν µ e 0.98 ν 0.96 e 0.97 0.94 0.1 1 2 −1. −.8 −.6 −.4 −.2 0 .2 .4 .6 .8 1. E ν (GeV) cosθ z FIG.1: Leftpanel: Ratioofalldirectionaveraged fluxwithasmallervirtualdetector(Φ (θ = 5◦)) ν d to that with the larger virtual detector (Φ (θ = 10◦)) used in HKKM04. Right panel: Zenith ν d angle (θ ) dependence of the ratio at E = 0.1 GeV in azimuthal average. z ν We can reduce the error due to the finite size virtual detector, with a little more com- putation. Let us assume the “true” atmospheric neutrino flux is expressed as an analytic function of the position. Note, we drop the arguments for arrival direction in the following expressions. The discussion here should apply to each arrival direction independently. We consider the power expansion the analytic function as, φ (θ ,θ ) = φ(0,0) +φ(1,0) ·θ +φ(0,1) ·θ +φ(2,0) ·θ2 +φ(1,1) ·θ θ +φ(0,2) ·θ2 +··· (1) ν x y ν ν x ν y ν x ν x y ν y 6 and 1 ∂m+nφ φ(m,n) ≡ · ν , (2) ν m!·n! ∂θxm∂θyn(cid:12)(θx,θy)=(0,0) (cid:12) where, θ ,θ are the distances from the target dete(cid:12)ctor in center angle to any directions x y perpendicular to each other, say, to South and East respectively, i.e., (θ ,θ ) constitute a x y local coordinate system. In the Monte Carlo calculation of the atmospheric neutrino flux, the calculated flux with a finite size virtual detector is the average flux over the virtual detector. With the increase of statistics, the flux Φ(θ ) calculated in Monte Carlo calculation should approach d 1 φ (θ ,θ )dθ dθ (θ ≪ 1), (3) S(θ ) Z ν x y x y d d θr<θd where S(θ ) ≃ πθ2 isthe“area”ofthevirtual detector, andθ ≃ θ2 +θ2. The integrations d d r x y p of terms proportional to θ or θ in Eq. 1 vanish, and non-vanishing terms start from the x y integrations of second order terms, φ(2,0)θ2 + φ(1,1)θ θ + φ(0,2)θ2, resulting in the terms ν x ν x y ν y proportional to θ4. For a sufficiently small θ , Φ (θ ) is expressed as, d d ν d Φ(2)θ4 Φ (θ ) ≃ Φ(0,0) + d = Φ(0,0) +Φ(2′)θ2 , (4) ν d S(θ ) d d where Φ(2′) ≡ Φ(2)θ2/S(θ ) ≃ Φ(2)/π. When we have the neutrino fluxes calculated with d d two virtual detectors with small enough radii θ and θ /2 for the same target, we expect d d Φ (θ )−Φ (θ /2) ≃ Φ(2′)·[θ2−(θ /2)2], thenΦ (0),thetruefluxvalueatthetargetdetector, ν d ν d d d ν is given as 4 Φ (0) = Φ(0,0) ≃ Φ (θ )− ·[Φ (θ )−Φ (θ /2)] . (5) ν ν d ν d ν d 3 As is seen in Fig. 1, the difference of the Φ (10◦) and Φ (5◦) is almost constant for ν ν cosθ > 0, so it should be sufficient to examine the assumption and procedure for vertical z down going directions. In the left panel of Fig. 2, we plotted the total neutrino flux for the vertical down going directions (cosθ > 0.9) calculated with different size of virtual z detectors, θ = 10◦,5◦, and 2.5◦ for Kamioka, Sudbury, and Gran Sasso with the HKKM04 d calculation. In the right panel of Fig. 2, we depicted the difference to the estimated true value with Eq. 5 in the ratio. We may say the convergence of the calculated fluxes to the “true value” agrees with the expectation of Eq. 4, and we apply the Eq. 5 with θ = 10◦ d and 5◦ to the atmospheric neutrino flux calculated in the 3-dimensional scheme. Note, the error due to the finite size of virtual detector does not exist in the 1-dimensional calculation scheme. 7 Vertical,E ν =100 MeV Vertical, E ν =100 MeV −1V ) 3 SNO 0) 1.1 −1sr Ge φ (ν −2−1m sec 2 Gansasso φ (θ )νd Kamioka 40 1 (x 1.0 Gran Sasso φ (θ )νd 1 Kamioka Sudbury 0 2.0 4.0 6.0 8.0 10. 0 2.0 4.0 6.0 8.0 10. θ (degree) θ (degree) d d FIG. 2: Left: Atmospheric neutrino fluxes for vertical directions calculated with the virtual detec- tors with different radii. Right: Ratio of fluxes calculated with virtual detectors with θ = 10,5, d and 2.5◦ to the flux estimated with Eq. 5. Thuscalculatedatmospheric neutrino fluxes inthe3-dimensionalscheme areshown inthe Appendix A up to 10 GeV for Kamioka, Sudbury, and Gran Sasso separately. The neutrino flux calculated for the Soudan2 site is almost identical to that calculated for Sudbury. The neutrino flux calculated for Frejus is ∼ 10 % larger than that for Gran Sasso at 0.1 GeV, and the difference is smaller at higher energies. Above 10 GeV, the atmospheric neutrino flux is calculated using the 1-dimensional scheme. They are tabulated in the Appendix B up to 10 TeV. We compared the atmospheric neutrino fluxes calculated with modified and original DPMJET-III in the ratio of flux values in Fig 3, averaging over all directions for Kamioka up to 1 TeV. The atmospheric neutrino flux calculated with the modified DPMJET-III shows an increase above 10 GeV from that with the original DPMJET-III, but the increase rate is different for different kinds of neutrinos. This is because the modification of the interaction model in Paper I enhances the productions of π±’s, K+’s, and K0, with no change for K− production. Therefore, the increase of ν and ν is larger than that of ν¯ and ν¯ . µ e µ e In Fig 4, we compared the atmospheric neutrin fluxes calculated with the modified and original DPMJET-III with those from other calculations based on the 3-dimensional calcu- lation scheme, Bartol [8, 24] and Fluka [19, 20]. The atmospheric neutrino fluxes calculated for Kamioka and averaged over all the directions are depicted in panel (a) and the ratios are compared in panel (b) up to 1 TeV. Although there are sizable differences in the flux values 8 −III T E J M P ν D 1.2 µ al n gi Ori −III/ νe T JE νµ M P D d e difi 1.0 ν +−5% Mo e −1 0 1 2 3 10 10 10 10 10 E ν (GeV) FIG. 3: The comparison of the atmospheric neutrino fluxes calculated with modified and original DPMJET-III in ratio. The denominator is the original DPMJET-III. among different calculations, the ratio (ν +ν¯ )/(ν +ν¯ ) is almost identical to each other µ µ e e below 100 GeV among them. 3 10 2V ) (a) os 5 −2−1−1m sec sr Ge 2 xνµ ννeµ utrino Flux Rati (b) THBFhlKauirKskto aMWl0o4rk ννµe ++ννµe ( 2 e 3ν10 N ( x 0.2) E φx ν 0.5 xνe 2 ν ν ( x1.5) µ µ 1 This Work 10 HKKM04 ν ν e e Bartol Fluka −1 0 1 2 3 1 −1 0 1 2 3 10 10 10 10 10 10 10 10 10 10 Eν(GeV) E ν (GeV) FIG. 4: The comparison of all direction average of the atmospheric neutrino fluxes with other calculations [8, 19, 20, 24]; (a) the absolute values of each kind of neutrinos and (b) the ratio of them. 9 III. MODIFICATION OF FLUKA’97 AND FRITIOF 7.02 Inthissection, we modifythe FLUKA’97andFritiof7.02interaction models, so thatthey reproduce the observed atmospheric muon data, following the procedure we used to modify DPMJET-III in Paper I. Then we calculate the resulting atmospheric neutrino flux to study the robustness of the modification procedure. The calculations in this section are carried out using the 1-dimensional calculation scheme for computation speed. The modification is applied to the hadronic interactions above 30 GeV, to study the muon flux above 10 GeV/c and neutrino flux above 3 GeV/c. > 1.7x Z K /Zπ N< 10 −1 Zπ 7 1. Z 3x ZK DPMJET−III 10 −2 FLUKA’97 Fritiof 7.02 1 2 3 4 5 6 10 10 10 10 10 10 E (GeV) proj FIG. 5: Z-factors for π++π− (Z ) and K++K− (Z ), and their ratio in each interaction model π K as the function of projectile energy. Solid lines are for DPMJET-III, dashed lines for FLUKA’97, and dotted lines for Fritiof 7.02. The secondary spectra of π and K productions differ between interaction models. The difference is seen in the Z-factors defined as p Z ≡ N < x1.7 > andx ≡ i , (6) i i i p proj where N is the multiplicity and p is the momentum of the i secondary particle, and p i i proj is the projectile momentum. The Z-factors are compared in the sums, Zπ+ + Zπ− and ZK++ZK−, for DPMJET-III, FLUKA’97, and Fritiof 7.02. in Fig. 5. Note, the contribution of neutral K’s to the neutrino flux is smaller than that of the charged K’s and π0 does not contribute to either muon flux or neutrino flux. The neutral π and K are not compared in the figure. 10