Improving the prediction of the Atmospheric neutrino flux using the atmospheric muon flux

It is well known that the correlation of atmospheric neutrinos and muons are simply correlated in the energy region of 1–10 GeV, and used for the test bench of the hadronic interaction model used for the calculation of the atmospheric neutrino flux. However, the correlation becomes unclear for neutrinos in the energy range below 1 GeV, which is important for the study of mass ordering of neutrino and CP phase of the neutrino mass. We extend the study of the correlation to the lower neutrino energies and find that the atmospheric muon flux observed at high altitude shows a good correlation to the atmospheric neutrino flux, and could be used to calibrate the hadronic interaction model.

Neutrino oscillation physics has now entered into the precision era, after it has been well established that θ 13 has a finite value and clear evidence for neutrino oscillations of one flavor to another have been seen. Now experimentalists are working hard to determine very precisely the neutrino mixing parameters, mass hierarchy and CP violation in the lepton sector, and tomography of the Earth. For those purposes, new generation detectors, such as INO [1], JUNO [2], South Pole [3], Hyper [4], DUNE [5] etc, are under the construction.
In those detectors, the atmospheric neutrino is still one of the main targets, and high precision measurements are expected, especially at low energies ( 1GeV). However, it is difficult to calculate the atmospheric neutrino flux below 1 GeV accurately, only with the results of high energy experiments, or with the observed atmospheric muon flux as we have done in Ref [6]. In this paper we study the correlation of atmospheric muons and neutrinos at low energies. Introducing a pseudo-analytic expression for the calculation we note that the simulation and the pseudoanalytic expression for the calculation of the atmospheric muon and neutrino are similar, but there are large differences between them. The important region in the phase space of hadronic interactions for atmospheric muons are different from that of the neutrino, and varies with the observation altitude. Therefore the atmospheric muon response to the variation of the hadronic interaction model is different from that of the atmospheric neutrino depending on the observation site, especially at low energies.
In the following sections, we study the correlation of atmospheric muons and neutrinos, and study the possibility to reduce the calculation error of the atmospheric neutrino flux.
Let us consider the projection of the atmospheric lepton flux calculation process to the hadronic interaction * e-mail: mhonda@icrr.u-tokyo.ac.jp phase space spanned by the projectile particle momentum and the momentum of the produced meson which creates the lepton in the decay cascade. Writing the whole process in a pseudo-analytic expression as eq. 1, we may consider the integrand is the projection for each combination of projectile and meson, where M2L(M, p born M , x born , L, p L , x obs ) represents the probability that a meson M born with momentum p born M at x born decays and results in the lepton L obs with momentum p obs L at x obs L without a hadronic interaction with air nuclei, H int (N pro j , p pro j N , M born , p born M ) is the probability that the M born with momentum p born M is produced in the hadronic interaction of a hadronic particle N pro j with momentum p pro j N with air nuclei, σ prod (N pro j , E pro j ) is the production cross section of N pro j particle and air nuclei, ρ air (x int is the nucleus density of the air at x int , and Φ pro j (N pro j , p pro j N , x int ) the flux of N pro j -particle primary or secondary cosmic ray at x int with momentum p pro j N , which could be the projectile particle in the hadronic interaction.
We note we can make the pseudo analytic expression eq.1 to a realistic one, and carry out an analytic calculation in a 1-dimensional calculation scheme with some approximations. However, it is not practical to carry out the analytic calculation in a 3-dimensional scheme, and normally the Monte Carlo simulation is used. We use eq. 1 just for illustrating the study of the uncertainty of the hadronic interaction model. In the following study, it is convenient to change the integral variable and rewrite eq.1 as, and The Monte Carlo simulation we use here is the same one which we have used for our calculation of atmospheric neutrinos. Tagging all the particles appearing in the simulation, we record two momenta (p born M , p pro j N ), one is that of the meson which created the lepton in the decay and the other is that of the projectile particle which created the meson in the hadronic interaction with the air nuclei. The set of two momenta (p born M , p pro j N ) is considered as a point in the phase space of the hadronic interaction. We show the projection of such points in the phase space in Fig. 1 for muon flux calculations and in Fig. 2 for neutrino flux calculations, summing all the different kinds of projectiles but for each kind of meson.
The difference of the two hadronic interaction models may be understood by the difference of the probability density in the phase space. As such a variation, we consider an expression; where ∆ int (N pro j , M born , p pro j N , p born M ) is a "variation function" defined for each combination of projectile particle and meson defined in the phase space of the hadronic interaction.
If we assume that the Φ pro j (N pro j , p pro j N , x int ) is not largely affected by the variation of the hadronic interaction, we can calculate the lepton flux with a small variation of the interaction model by substituting the corresponding probability functionH int (N pro j , p pro j N , x int , M born , p born M ) in the expression eq.1. It means we can study a small variation of the interaction model and resulting atmospheric lepton flux without repeating the Monte Carlo simulation.
The lepton flux calculated by the variation of the interaction model may be written as; In the following, we study the effect of the the variation of hadronic interaction model, introducing an artificial variation in the interaction model, with a set of random EPJ Web of Conferences 208, 07001 (2019) https://doi.org/10.1051/epjconf/201920807001 ISVHECRI 2018 numbers and 3rd order b-spline. The 3rd order b-spline function with constant knots separation is represented as where ∆ is the knots separation, and x 0 is for a shift less than ∆. The linear combination of 3rd order b-spline function (eq.7) is contiguous up to 2nd order derivative, and is often used to connect the discrete data points or to fit them. We can extend it to a 2-dimensional spline function as Here, we took the same knots separation for x and y. We may construct an artificial variation function in eq.4 with those 2-dimensional b-spline functions as where {R i j N } is the set of random numbers with a normal distribution. Therefore, for each combination of the set of random numbers {R i j N } and δ, we have a variation of the interaction model and a flux spectrum as, Note, we assume δ has a small value ( 1) so that we consider a small variation from the hadronic interaction we used in our calculation of atmospheric lepton calculations. With this expression (eq.11), we can write the variation of the lepton flux as Once we have calculated the atmospheric lepton flux with a Monte Carlo simulation, it is easy and quick to calculate the variation due to a variation of the hadronic interaction model with the above formalism. Applying this Let us study the correlation of atmospheric neutrinos and muons for the variation of hadronic interaction made as in eq.11 using the correlation coefficient; In Fig.3, we show the correlation coefficient as a function of vertical down going muon momentum observed at Kamioka, and for the vertical down going 1 GeV neutrino also observed at Kamioka. We consider here all the possible combination of the kind of neutrinos (ν µ ,ν µ , ν e ,ν e ) and of muons (µ + , µ − ) for the correlation here. In the left panel of Fig.4 we plotted the muon momenta at the maximum correlation and 90% of it for the combinations of neutrinos and muons, but when the correlation coefficient is larger than 0.3. We find that the muon momenta at the maximum correlation is about 3 times the neutrino energy above a few GeV, but it quickly decreases as the neutrino energy decreases.
Next we take δ = 1 and construct 3000000 set of random numbers {R i j N } to study the variation of muons and neutrino fluxes. In the right panel of Fig.4, we plot the variation of the neutrino flux calculating the ratio ∆Φ ν µ /Φ ν µ at 1 GeV, and the same quantity when the variation of interaction model satisfies ∆Φ µ /Φ µ < 0.1 for the calculated muon flux in the momentum range having a correlation coefficient larger than 90% of the maximum. We find that with the selection of hadronic interaction model variation with the muon flux variation, the variation of the neutrino flux is largely suppressed and the standard deviation is around 0.5.
With this observation, we may consider that when we reconstruct the accurately observed atmospheric muon flux accurately in the muon momentum region where the correlation coefficient of muon flux and the neutrino flux is larger than 90% of the maximum, the expected error of the neutrino calculation is also largely suppressed. and muons (µ + , µ − ) both for vertical down going direction, as a function of neutrino energy. When the maximum correlation coefficient is smaller than than 0.3, the lines are not plotted. Right panel: The outside curve is the distribution of ∆Φ νµ /Φ νµ at 1GeV with the variation of hadronic interaction model by eq.12 for 3000000 set of random numbers and δ = 1. The inside curve is the same as above, but selected for the case when the hadronic interaction model which gives ∆Φ µ /Φ µ < 0.1 is satisfied in the muon momentum region with the correlation coefficient larger than 90% of the maximum. However in the practical sense, it is difficult to measure the atmospheric muon flux accurately in the low momentum region. Therefore, we consider the cases that the available minimum muon is 0.1, 0.3, 0.5, and 1.0 GeV/c for the accurate muon flux data. Then we plotted the variation of the atmospheric neutrino variation standard deviation in the ratio to flux value ∆Φ ν µ /Φ ν µ for the variation of muon flux satisfying ∆Φ µ /Φ µ < 0.1 in Fig.5, both for muon fluxes observed at Kamioka (sea level) and Mt. Norikura (∼2700m A.S.L.). We find the suppression of the neutrino flux error works better at the higher observation altitude.

Summary
We have developed a quick method to study the effect of the variation of hadronic interaction model in the calculation of the atmospheric lepton flux. Using this method, we find a procedure to reduce the error in the calculation of the atmospheric neutrino flux at low energies with the accurately measured atmospheric muon flux. However, the observation site for muon flux seems important. The muon flux observed at Mt. Norikura (∼ 2700m A.S.L.) gives a more stringent suppression of neutrino flux error than that at Kamioka (sea level). It is expected that the observations of muon flux on a much higher mountain will give a more stringent suppression of atmospheric neutrino calculation errors.