Observation of Atmospheric Neutrinos

In 1998, the Super-Kamiokande discovered neutrino oscillation using atmospheric neutrino anomalies. It was the first direct evidence of neutrino mass and the first phenomenon to be discovered beyond the standard model of particle physics. Recently, more precise measurements of neutrino oscillation parameters using atmospheric neutrinos have been achieved by several detectors, such as Super-Kamiokande, IceCube, and ANTARES. In addition, precise predictions and measurements of atmospheric neutrino flux have been performed. This paper presents the history, current status, and future prospects of the atmospheric neutrino observation.


History
Atmospheric neutrinos are generated when primary cosmic rays strike nuclei in the atmosphere and the hadron that results from these collision decays [1]. The principles and expected event rates of neutrino detection were first proposed in the 1960s [2,3]. The first discovery was achieved by two independently working groups in 1965. These groups placed detectors deep underground in the Kolar Gold Mines of South India [4] and a South African gold mine [5], respectively, to avoid cosmic-ray muons. They looked for upward-going muon events, because it must be generated through the interaction of atmospheric neutrino from the other side of the Earth with the surrounding rock, and found the signal clearly.
The next generation of neutrino detectors, Kamiokande and IMB (Irvine-Michigan-Brookhaven detector), appeared in the 1980s. These were kiloton-scale water Cherenkov detectors, whose detection principles will be mentioned later. The original purpose of the experiments using these detectors was to search for nucleon decay predicted by the Grand Unified Theory, and atmospheric neutrinos were one of the serious backgrounds of the nucleon decay search. However, nucleon decay was not discovered, thus, neutrinos became the main target of research. The successful of the neutrino observation was, at first, neutrino signals from Supernova 1987A [6]. In addition, Kamiokande detected clear neutrino signal from the Sun [7], called "solar neutrinos". The first discovery of the solar neutrinos were reported by Homestake experiment in 1968 [8]. The strength of the observed signal was about one third of the predicted amount. It was long standing problem called "solar neutrino puzzle" [9]. At the beginning of the 2000s, this was found to be attributable to neutrino oscillation.
The anomalous atmospheric neutrino measurements were also reported by several experiments in this period. The ratio of ν µ and ν e should be roughly 2:1 since atmospheric neutrinos are generated by the decay of pions (ν µ ) and muons (ν µ and ν e ); however, the ratio observed was not in agreement with the predicted ratio. The average discrepancy between the data and Monte Carlo simulation (MC) was (µ/e) data /(µ/e) MC = 0.57 +0.08 −0.07 ± 0.07 in Kamiokande, though the ratio had incident angle dependence [10,11]. The value reported in IMB was 0.54 ± 0.05 ± 0.12 [12,13] which was consistent with

Atmospheric Neutrino Flux Prediction
Atmospheric neutrinos are generated from the decay of π and K, which are secondary particles resulting from the interaction of cosmic rays with air molecules in the Earth's atmosphere. They come from all directions into detectors although their flux depends on the zenith and azimuth angles of the direction of their arrival. The flux in the horizontal direction is generally higher than that in the vertical direction, because of the longer path taken by parent particles through the atmosphere. The energy from atmospheric neutrinos exists in a wide range of 100 MeV to PeV scale measurements along a power-law spectrum, as shown in the Figure 1. The reason fewer neutrinos are produced at higher energies is mainly due to the falling primary cosmic-ray spectrum. The effect that the π decay lengths are longer than the paths in atmosphere and the parent particles reach the ground before decaying, plays a role. It is also relevant for the energy dependence of ν µ and ν e ratio and causes the neutrino flux to peak at the horizon as shown in the right Figure 1. The energy distribution is suppressed below the GeV energy region due to the rigidity cutoff effect of the primary cosmic rays by Earth's magnetic field. The trajectories of charged particles in primary cosmic rays are affected by geomagnetic fields. Only particles that interact in the atmosphere before curving back into the space produce atmospheric neutrinos. The trajectory depends on the particles' momentum and total charge; the ratio of them is called "rigidity". Geomagnetic fields affect particles with lower energy more strongly; therefore, low-energy atmospheric neutrino flux is suppressed, a process that is called "cutoff". Since Kamioka is located at a rather low geomagnetic latitude, it has a high local vertical rigidity cutoff. In addition, due to its strong azimuthal asymmetry, the cutoff is higher for particles arriving from the east than from the west. This east-west asymmetry arises from the fact that the primary cosmic rays are positively charged. On the other hand, in high-energy regions, above a few tens TeV, neutrinos from charm mesons' decay are considered, instead of π and K decay, due to the much shorter lifetimes of these charm mesons which are on the order of 10 −12 s. These are called "prompt" neutrinos [18,19], which are uniformly generated in the atmosphere, with equal fluxes of ν µ and ν e .
Precise predictions of the atmospheric neutrino have been performed by several research groups, including HKKM [20], Bartol [21], and FLUKA group [22]. The differences among these models are choices of hadronic models and measurements of the primary cosmic-ray spectrum. Figure 1 shows a comparison of atmospheric neutrino fluxes and neutrino flux ratios calculated by different groups. The flux calculations differ by about 10%. Another feature is that the flavor ratio between ν µ +ν µ and ν e +ν e below the GeV scale is approximately two; however, it increases with the energy levels, as shown in Figure 1, because the muons produced by π decay reach the ground before decaying. Time variation in the atmospheric neutrino flux is also expected. Long-term variation is due to solar activity with an average period of 11 years. The primary cosmic-ray flux at Earth is anticorrelated with the solar activity because the plasma from the Sun scatters the cosmic rays, and the cosmic-ray flux is reduced during periods of high solar activity. Consequently, the atmospheric neutrino flux is also predicted to be anticorrelated with solar activity. There is also yearly variation due to seasonal temperature variations that affect atmospheric density which increases in summer, when relatively more neutrinos are produced.
calculation is smoothly connected to the 1D calculation carried out in the previous work. As the modified JAM is used below 32 GeV, any difference above that is due to the difference of the calculation scheme between 3D and 1D. However, the difference between present and previous works is very small in the figure above 1 GeV. Note, we magnify the difference between present and previous works in the ratio in Fig. 8. The difference is less than a few percent above 1 GeV. On the other hand, the atmospheric neutrino fluxes calculated with the modified-JAM show an increase from the previous one below 1 GeV, as is expected from the increase of atmospheric muon spectra below 1 GeV=c.
In the right panel of Fig. 7, we show the ratios of the atmospheric neutrino flux averaged over all directions. It is seen that the differences in the flux ratios are very small between present and previous works above 1 GeV as the absolute atmospheric neutrino fluxes. It is noticeable that the ð þ Þ =ð e þ e Þratios are very close to each other including the ones calculated by different authors. We find there is a visible difference in e = e between the present and previous works. As is seen in the left panel of Fig. 5 in Sec. II, the original JAM interaction model has a smaller þ = À ratio, and is responsible for the smaller e = e ratio below 1 GeV. Note, we do not modify the interaction model more than the muon flux data require. It is difficult to examine the þ = À ratio at the momenta corresponding to the e = e below 1 GeV, due to the small statistics in the observation at balloon altitude.
The muons at sea level or mountain altitude are not useful to examine the atmospheric neutrino of these energies, since the muons result from higher energy pions at higher altitude.
In Fig. 9, we show the atmospheric neutrino fluxes as a function of the zenith angle averaging over all the azimuthal angles at 3 neutrino energies; 0.32, 1.0, and 3.2 GeV for Kamioka. In Fig. 10, we show the comparison of the present and previous works in the ratio as the function of zenith angle. There is a difference due to the increase of the  [20,21], dotted lines for the FLUKA group [22], and dash dot for our previous calculation (HKKM06).
IMPROVEMENT OF LOW ENERGY ATMOSPHERIC . . . PHYSICAL REVIEW D 83, 123001 (2011) the variation of atmospheric neutrino flux as a function the geomagnetic field [23]. Note, the variation of upward    [20], Bartol (dashed) [21], FLUKA (dotted) [22], and a previous HKKM06 model [23]. Some plots are applied by several factors for easy viewing. Lower figures show zenith angle dependence of atmospheric neutrino flux averaged over all azimuthal angles for Kamioka site calculated by HKKM11. These figures are taken from [20].

Neutrino Interaction
Interactions with nuclei in water (or ice) and rocks surrounding the detector are the norm in atmospheric neutrino observations, and that with electrons is negligible as the cross-section is three orders of magnitude smaller than that with nuclei. Interactions can be classified into charged-current (CC) or neutral-current (NC) interactions according to the type of bosons that are exchanged, W ± or Z 0 , Universe 2020, 6, 80 4 of 24 respectively. A CC interaction produces a charged lepton, electron or muon, whose flavor corresponds to that of a neutrino, ν e or ν µ . Therefore, the original neutrino flavor is identified by distinguishing the flavor of the related charged lepton. However, an NC interaction does not indicate the neutrino flavor since the outgoing lepton is a neutrino. The following neutrino interactions are dominant in the atmospheric neutrino energy region,

Neutrino Oscillation
Neutrino oscillation occurs if flavor eigenstates are mixed with mass eigenstates, and a difference in mass exists. The mixing between flavor eigenstates (ν α ) and mass eigenstates (ν i ) can be written as Universe 2020, 6, 80

of 24
In the three-flavor neutrino framework, U is Pontecorvo-Maki-Nakagawa-Sakata (PMNS) mixing matrix [44][45][46], and usually parametrized by three mixing angles (θ 12 , θ 23 , θ 13 ), and one CP-violating Dirac phase (δ CP ) in neutrino oscillations as follows: where c ij and s ij represent cos θ ij and sin θ ij , respectively. Neutrino oscillation frequencies are determined by the neutrino mass differences, ∆m 2 21 = m 2 2 − m 2 1 and ∆m 2 32 = m 2 3 − m 2 2 , where m 1 , m 2 , and m 3 are the three mass eigenvalues. Among these oscillation parameters, θ 12 and ∆m 2 21 have been measured by solar and reactor neutrino experiments, θ 23 and ∆m 2 32 have been measured by atmospheric and accelerator neutrino experiments, and θ 13 has been determined by reactor neutrino experiments. The order of absolute mass (m 3 m 2 > m 1 in a normal ordering or m 2 > m 1 m 3 in an inverted ordering) and the CP phase (δ CP ) are still unknown parameters. The dominant neutrino oscillation in atmospheric neutrinos is a channel between ν µ and ν τ caused by the parameters of the mass eigenstate between ν 2 and ν 3 . In addition, the appearance of oscillation from ν µ to ν e is considered to be a sub-leading effect. The dominant neutrino oscillation probabilities in a vacuum for atmospheric neutrinos are expressed as follows: (4) P(ν µ ↔ ν e ) sin 2 θ 23 sin 2 2θ 13 sin 2 ∆m 2 31 L 4E (5) P(ν µ ↔ ν τ ) sin 2 2θ 23 cos 4 θ 13 sin 2 ∆m 2 31 L 4E (6) P(ν e ↔ ν τ ) cos 2 θ 23 sin 2 2θ 13 sin 2 ∆m 2 31 L 4E (7) where L is the flight length of neutrinos and E is the neutrino energy. When neutrinos traverse the Earth, their matter potential due to the difference in the forward-scattering amplitudes of ν e and ν µ,τ , which induces a matter dependent effect on of neutrino oscillations [47], must be taken into account. In this scenario, the oscillation value in Equation (7) where V e = ± √ 2G F N e is the effective matter potential, and the sign is positive for neutrinos or negative for anti-neutrinos; N e is the electron density, which is assumed to be constant; and G F denotes the Fermi constant. The potential is derived from the difference that ν e which has both CC with electrons and NC with electrons and nucleons, while ν µ and ν τ have only NC. In this form, when 2EV e /∆m 2 31 = cos 2θ 13 , the effective mixing angle is resonantly enhanced. Since cos 2θ 13 is positive, the enhancement only occurs for neutrinos if ∆m 2 31 is positive which is as normal mass ordering, while it occurs for anti-neutrinos only in the case of the inverted mass ordering. Figure 3 shows the ν µ survival and Universe 2020, 6, 80 6 of 24 ν µ → ν e transition probabilities for neutrinos and anti-neutrinos assuming a normal mass ordering. Earth matter effects suppress the disappearance of ν µ and enhance the appearance of ν e especially in upward-going neutrinos, where the cosine zenith angle is negative, with energies in the range of 2-10 GeV. The appearance of ν e in neutrino (b) in Figure 3 is enhanced in this energy region, while no clear enhancement appears in the anti-neutrino plot (d). If the mass ordering is inverted, this feature is switched and appears as an anti-neutrino case. Therefore, the difference in this level of enhancement between neutrinos and anti-neutrinos can help to determine the mass ordering in neutrino oscillation. The ν µ → ν e transition probability is also affected by the δ CP parameter, which results in a change of about 2% in the maximum total ν e flux observed in the energy region less than 1 GeV. quantity of underground water available to fill the detector and maintain its temperature. These changes impact the water transparency and subsequent performance of the detector and therefore must be corrected through calibrations. Since neutrino oscillations are a function of the neutrino energy, a thorough understanding of the detector energy scale is important for precision measurements.
At the same time the range of energies of interest to atmospheric neutrino analysis spans from tens of MeV to tens of TeV, eliminating the possibility of calibration through radioactive isotopes. Accordingly, the energy scale is calibrated using natural sidebands covering a variety of energies. Neutral pions reconstructed from atmospheric neutrino interactions provide a calibration point via the π 0 momentum and stopping cosmic ray muons of various momenta are used to measure photoelectron production as a function of muon track length (Cherenkov angle) for multi-GeV (sub-GeV) energies. Here the muon track length is estimated using the distance between the entering vertex and the position of the electron produced in its subsequent decay. The energy spectrum of these Michel electrons additionally serves as a low energy calibration point. Figure 3 shows the absolute energy scale measurement using each of these samples.
In the oscillation analysis the absolute energy scale uncertainty is conservatively taken to be the value of the FIG. 2. Oscillation probabilities for neutrinos (upper panels) and antineutrinos (lower panels) as a function of energy and zenith angle assuming a normal mass hierarchy. Matter effects in the Earth produce the distortions in the neutrino figures between two and ten GeV, which are not present in the antineutrino figures. Distortions in the ν μ survival probability and enhancements in the ν e appearance probability occur primarily in angular regions corresponding to neutrino propagation across both the outer core and mantle regions (cosine zenith < −0.9) and propagation through the mantle and crust (−0.9 < cosine zenith < −0.45). For an inverted hierarchy the matter effects appear in the antineutrino figures instead. Here the oscillation parameters are taken to be Δm 2 32 ¼ 2.5 × 10 −3 eV 2 , sin 2 θ 23 ¼ 0.5, sin 2 θ 13 ¼ 0.0219, and δ CP ¼ 0.

Super-Kamiokande
The Super-Kamiokande (SK) detector is a cylindrical tank of 39.3 m diameter and 41.4 m height, filled with 50 kilotons of pure water as shown in Figure 4. It is located 1000 m underground (2700 m water equivalent) in the Kamioka mine in Gifu Prefecture, Japan. The detector is divided into two regions called "inner" and "outer", and lined with 11,129 twenty-inch PMTs in the inner detector and 1885 eight-inch PMTs in the outer detector. The signal detection method of SK is that Cherenkov light generated in water from the charged particles is observed by PMTs [49]. The fiducial volume used for the data analysis is 22.5 kilotons, which is within 2 m of the inner wall to maintain the detector performance. SK was launched in April 1996, and until now, there have been five experimental phases. The first phase lasted five years until July 2001. There were several successes in the detection of atmospheric neutrinos and solar neutrinos. However, a serious accident occurred in November 2001, in which most of the PMTs were broken. The experiment was resumed in October 2002, with about half the number of PMTs. Fortunately, even with half the number of PMTs, the atmospheric neutrino analysis was not affected greatly. After a three-year operation, full reconstruction was completed in 2006. The third phase, with almost the same number of PMTs as in the first phase, started in July 2006, and ended in August 2008. In the fourth phase, new electronics (QBEE [50]) were installed, which improved the detection efficiency of the decay electron from stopping muon higher. It is also possible to tag a 2.2 MeV gamma-ray, which is produced by the neutron capture by protons in water, even though it is about 20% detection efficiency. Neutron signals are useful for the improvement of atmospheric neutrino analysis. Although it is difficult to separate the neutrino and anti-neutrino events, information about the number of neutrons can be used to statistically differentiate between them because the number of neutrons in anti-neutrinos is larger than that in neutrinos. The fourth phase continued for 10 years until May 2018. SK was refurbished during the rest of 2018 to allow the loading of gadolinium to pure water. If gadolinium is loaded, the efficiency of neutron detection significantly improves by up to 90% depending on the gadolinium concentration. The main purpose of this gadolinium loading is the discovery of Diffuse Supernova Neutrino Background, but it also improves the analysis of atmospheric neutrinos due to the high neutron detection efficiency. After a successful refurbishment, the fifth phase started in January 2019. In the atmospheric neutrinos observation in SK [51], the data sample is categorized by three distinct topologies: fully contained (FC), partially contained (PC), and upward-going muon (UPMU). FC events produce reconstructed interaction vertices inside the fiducial volume, while PC events also produce interaction vertices within the fiducial volume of the inner detector but are accompanied by considerable light in the outer detector. FC can be further sub-divided into single-and multi-ring and electron-and muon-like events. UPMU events are caused by muon-neutrino interactions in the surrounding rock, which produce penetrating muons. These muons either stop in the inner detector volume (stopping events) or continue through the inner detector (through-going events). The leptons, muons or electrons, generated by neutrino interactions preserve information from the original neutrinos, including the neutrinos' flavor, energy, and direction. To identify muons or electrons, the Cherenkov ring pattern is used. An electron produces diffused ring patterns because of the electromagnetic shower and multiple scattering, while a muon produces a ring with a sharp edge. The SK atmospheric neutrino analysis uses event categories, reconstructed energy and direction, and particle types (muon-like or electron-like).
Discovery of the Neutrino Oscillation Figure 5 shows several plots of zenith angle dependence on the first results obtained from SK in 1998 [15]. The electron-like events were consistent with the predicted results, while the number of upward-going muon-like events were significantly smaller than that of the downward-going events. This was evidence of the neutrino oscillation that muon-neutrinos converted to other flavor of neutrinos through the flight inside the Earth.
As a cross-check of the above analyses, we have reconstructed the best estimate of the ratio LyE n for each event. The neutrino energy is estimated by applying a correction to the final state lepton momentum. Typi-events show a significant deficit at large LyE n . At large LyE n , the n m have presumably undergone numerous oscillations and have averaged out to roughly half the initial rate.
The asymmetry A of the e-like events in the present data is consistent with expectations without neutrino oscillations and two-flavor n e $ n m oscillations are not favored. This is in agreement with recent results from the CHOOZ experiment [22]. The LSND experiment has reported the appearance of n e in a beam of n m produced by stopped pions [23]. The LSND results do not contradict the present results if they are observing small mixing angles. With the best-fit parameters for n m $ n t oscillations, we expect a total of only 15-20 events from n t chargedcurrent interactions in the data sample. Using the current sample, oscillations between n m and n t are indistinguishable from oscillations between n m and a noninteracting sterile neutrino. Figure 2 shows the Super-Kamiokande results overlaid with the allowed region obtained by the Kamiokande 1566 Figure 5. Zenith angle distributions of muon-and electron-like events at sub-GeV and multi-GeV scales in the first 535 days of SK data produced [15]. Upward-going events are at less than zero and downward-going events are at more than zero on the horizontal axis. The vertical axis shows the number of events in each bin. The hatched regions show the expected results without neutrino oscillation. The lines show the best-fit expectations with neutrino oscillation.

IceCube
The IceCube is a cubic-kilometer neutrino detector buried in the Antarctic ice, as shown in Figure 6. It comprises 5160 digital optical modules (DOMs) along 86 vertical strings, with 60 DOMs per string. Each DOM houses a downward-facing 10-in PMT with electronics in a glass pressure sphere. Of the 86 strings, 78 are deployed with an inter-string distance of ∼125 m and a space of ∼17 m between each DOM at depths between 1450 and 2450 m below the surface. This part of the detector is optimized for the neutrino energy range of 100 GeV to 100 PeV. The remaining eight strings, located at the bottom center of the detector, are set more densely with DOMs. This detector, called DeepCore, comprises 647 DOMs with high-quantum-efficiency PMTs placed 2100 m under the clearest ice. DeepCore has a volume of ∼10 7 m 3 , and is optimized for the detection of lower-energy neutrinos down to 5.6 GeV.
IceCube observes Cherenkov light generated in ice by charged particles resulting from neutrino interactions. As the target neutrino energy for IceCube is high, deep inelastic scattering is a dominant feature in neutrino interaction. The observed events are categorized into two types of patterns: long, straight tracks produced by muons (track-like) and spherical cascades produced by electromagnetic and/or hadronic showers (cascade-like).
IceCube observation was fully commissioned in 2011. One important result was the discovery of two ultrahigh-energy (∼PeV) neutrinos in 2013 [52]. Compared to the expected number of atmospheric neutrinos, 0.082 ± 0.004 +0.041 −0.057 , these two events were possibly of astrophysical origin.
lactic nuclei sources assuma-ray flux observed by the atory [5] is produced by the ertae objects that emit TeV ted to be optically thin to The model calculated by t charged cosmic rays are ough the decay of escaping of neutrinos from the pregamma-ray bursts is calcuthe gamma-ray emission to energy cosmic rays. n the search for diffuse asheric muons and neutrinos ced extensive air showers. ng atmospheric muon backenergy range from primary a GeV to the highest mea-100 EeV [8]. These events arth as a filter in order to os traversing through the eric neutrinos were conside decay of pions and kaons neutrino flux) and neutrinos rm-containing mesons (the o flux). Detailed threethe energy spectrum and

ANTARES
The ANTARES neutrino telescope is in the Mediterranean Sea at a depth of 2475 m [54]. A schematic view of the detector is shown in Figure 7. It comprises 12 detection lines: 11 equipped with 25 storeys of three optical modules and one line with 20 storeys of optical modules, giving a total of 885 optical modules. Each optical module has a 10-in PMT, whose axis points 45 • downward. The detector was completed in 2008, and a total of 2830 days of data had been analyzed for atmospheric neutrino data by 2016.  Table 1 shows the summary of three detectors for atmospheric neutrino research.

Super-Kamiokande
The main contribution to the understanding of neutrino oscillation from atmospheric neutrinos is the determination of θ 23 and ∆m 2 32 . Recently, all sub-leading effects of neutrino oscillation can be investigated due to precise observation and large statistics of data. Here, the earth matter effect plays an important role in neutrino oscillation schemes, which resolves mass ordering, two possible θ 23 regions, and δ CP . It has different behavior of appearance between ν e andν e for the mass ordering. However, SK is insensitive to the charge sign of particles; therefore, CC neutrino and anti-neutrino interactions cannot be distinguished on an event-by-event basis. Instead, they are statistically separated based on the number of decay electrons, number of Cherenkov rings, and transverse momentum. In the most sensitive energy region between 2 and 10 GeV, not only CC quasi-elastic interactions but also single pion production via ∆ resonance excitation and the deep inelastic scattering process should be considered. In single pion production, π − generated in an anti-neutrino reaction, such as ν e + n → e + nπ − , will be captured on an 16 O nucleus, leaving the positron as the only detected particle and no delayed electron signal. In neutrino reactions, on the other hand, π + is generated via ν e + n → e − nπ + . It is not captured in this manner and produces a delayed electron signal through its decay chain. Thus, an anti-neutrino tends to produce a single-ring event without any delayed electron, while a neutrino event has the opposite effect. Due to the V-A structure of the weak interaction, the angular distribution of the leading lepton fromν is more forward than those from ν reaction, which means that the transverse momentum inν is expected to be smaller than that in ν. The statistical separation of ν e andν e is performed by a likelihood method using the above variables.
For a constraint on the neutrino oscillation parameter, the data obtained in SK are fitted to expectations by MC simulation using a binned χ 2 method. There are 520 analysis bins in total (energy and zenith angles in each event category) for each SK phase and 155 systematic uncertainties. Figure 8 shows the chi-square differences from the minimum as a function of ∆m 2 32 (or ∆m 2 13 ), sin 2 θ 23 and δ CP , for the normal and inverted mass ordering cases using 5326 days of data from SK measurements [48]. From the results, the normal mass ordering showed better agreement than the inverted mass ordering with ∆χ 2 ≡ χ 2 NH,min − χ 2 I H,min = −3.48. The best-fit neutrino oscillation parameters in the normal mass ordering were ∆m 2 32 = 2.50 +0.13 −0.31 × 10 −3 , sin 2 θ 23 = 0.587 +0.036 −0.069 , δ CP = 4.18 radians. In addition, a neutrino oscillation analysis with constraints from other experiments was provided. One concerned reactor short-baseline neutrino experiments, Daya Bay, RENO, and Double Chooz, for θ 13 , and the central value of these experiments was sin 2 θ 13 = 0.0219 ± 0.0012 [55]. The other concerned Tokai-to-Kamioka (T2K) long-baseline neutrino experiment [56]. Since SK is the far detector of T2K, many experimental aspects, for example, the detector simulation, the neutrino interaction generator (NEUT), and the event reconstruction tools, are common between them. Thus, it is possible add published binned T2K data to the SK atmospheric neutrino fit for determination of the neutrino oscillation parameters. When the constraints on θ 13 from the reactor experiments were applied to atmospheric neutrinos, the preference for the normal mass ordering was ∆χ 2 = −4.33. When the T2K constraints were added, it became stronger at ∆χ 2 = −5.27.
A. Results and discussion Figure 9 shows one-dimensional allowed regions for jΔm 2 32;31 j, sin 2 θ 23 , θ 13 and δ CP . In each plot the curve is drawn such that the χ 2 for each point on the horizontal axis is the smallest value among all parameter sets including that point. When the atmospheric neutrino data are fit by themselves with no constraint on θ 13 , the normal hierarchy hypothesis yields better data-MC agreement than the inverted hierarchy hypothesis with χ 2 NH;min − χ 2 IH;min ¼ −3.48. The preferred value of sin 2 θ 13 is 0.018(0.008) assuming the former (latter). Though both differ from the The data's preference for both nonzero sin 2 θ 13 and the normal mass hierarchy suggest the presence of upwardgoing electron neutrino appearance at multi-GeV energies driven by matter effects in the Earth (cf. Fig. 2). Figure 10 shows the up-down asymmetry of the multi-GeV singleand multiring electronlike analysis samples. Here the asymmetry is defined as N U − N D =N U þ N D , where N U ðN D Þ are the number of events whose zenith angle satisfy cos θ z < −0.4 (cos θ z > 0.4). Small excesses seen between a few and ten GeV in the multi-GeV e-like ν e and the multiring e-like ν e andν e samples drive these preferences. . Constraints on neutrino oscillation parameters from atmospheric neutrino data in SK [48]. The difference in the χ 2 distributions as a function of ∆m 2 32 (left), sin 2 θ 23 (center), and δ CP (right) assuming normal (blue) and inverted (orange) hierarchies. The minimum χ 2 value in the normal mass ordering is set to zero.

IceCube
Measurements of the atmospheric neutrino oscillation parameters of θ 23 and ∆m 2 32 in IceCube were reported in 2018 [57]. The unique characteristics of IceCube detector for the atmospheric neutrino search are, at first, the high statistics due to the huge volume. It is also possible to detect downward-going muon events induced by atmospheric neutrino interaction using the surrounding IceCube detector to distinguish with cosmic-ray muons. To provide a constraint for the neutrino oscillation parameters, the reconstructed neutrino energy (8 bins) and zenith angle (8 bins), both track-like and cascade-like, were fit to the expectations using a binned χ 2 method. Here, the reconstruction was performed by calculating the likelihood of the observation of photoelectrons by DOMs as a function of the neutrino interaction position, direction and energy. The best-fit neutrino oscillation parameters were ∆m 2 32 = 2.31 +0.11 −0.13 × 10 −3 , sin 2 θ 23 = 0.51 +0.07 −0.09 assuming a normal mass ordering. Figure 9 shows the L/E distribution along with the corresponding predicted counts given the best-fit neutrino oscillation parameters broken down by event type for both track-like and cascade-like events. The two peaks corresponded to down-going and up-going neutrino trajectories. As the track-like sample was enriched in ν µ CC events, up-going events were strongly suppressed in track-like events due to neutrino oscillation. While the cascade-like sample was evenly divided between ν µ CC events and interactions without a muon in the final state, some suppression could also be observed.
ting the relative optical ident photon angle. The n is modeled by two the shape of the DOM e lateral angular accephotons traveling roughly ontal) and is fairly well data. Five MC data sets to þ1σ uncertainty from arametrized in the same ncy described above. A data is used. s sensitivity to photons triking the DOMs head by string-to-string LED d using a dimensionless responding to a bubble DOM face for vertically tion). Zero corresponds to idence from 0°to 30°from range from −5 to 2 were

ANTARES
The measurement of the atmospheric neutrino oscillation parameters of θ 23 and ∆m 2 32 in ANTARES was reported in 2019 [58]. The energy for the atmospheric neutrino analysis in ANTARES ranged from a few tens of GeV up to 100 GeV. Track-like events originating from the penetrating muons produced via CC interactions of ν µ were used for the analysis. On the other hand, shower-produced events, for example, electromagnetic showers from ν e CC interactions or hadronic showers from NC interactions, were regarded as the background events. Here, muon-track reconstruction was essential, and two different algorithms were used: one in which PMT hits were selected to find the best muon track and another involving fitting to a chain at each step to improve the track estimation. Once the muon track was reconstructed, the muon energy was estimated from its track length, given a constant energy loss of muons in the sea at 0.24 GeV/m in the energy range of 10-100 GeV.
To obtain a constraint for the neutrino oscillation parameters, a logarithmic base-10 scale of reconstructed neutrino energy in GeV was divided into eight bins, seven from 1.2 to 2.0 plus an additional underflow bin for log 10 (E/GeV) < 1.2. The cosine zenith angle was divided into 17 bins from 0.15 to 1.0. The fit was performed using a log-likelihood approach. The best-fit neutrino oscillation parameters were ∆m 2 32 = 2.0 +0.4 −0.3 × 10 −3 , θ 23 = 45 +12 −11 degree. Figure 10 shows the reconstructed energy divided by the cosine of the reconstructed zenith angle. The lowest bin of this plot is expected to be affected by neutrino oscillation. The data showed good agreement with neutrino oscillation assumptions.  Figure 11 shows the allowed region of neutrino oscillation parameters (sin 2 θ 23 and ∆m 2 32 ) at 90% C.L. overlaying several experiments both for atmospheric and accelerator neutrinos. Even though there are different sources of neutrinos, and the sensitive energy range is different in atmospheric neutrino experiments, all results showed good agreement with each other.
Here, the MINOS far detector, which is 5.4 kTon mass of iron-scintillator tracking calorimeter, is used for the study of atmospheric neutrinos as well as neutrinos originating from the Fermilab NuMI accelerator beam. The far detector is located at underground (2070 m water equivalent) in Soudan mine, MN, USA. The detector is magnetized and it enables the separation of ν µ andν µ on an event-by-event basis using the curvature of the produced charged muon. The contour of MINOS in the figure was combined results of accelerator neutrinos and 60.75 kt-year data of atmospheric neutrinos. Figure 11. Contours of the 2D allowed region of ∆m 2 32 and sin 2 θ 23 at 90%C.L. from several experiments. Atmospheric neutrinos in SK (green) [48], IceCube/DeepCore (red) [57] and ANTARES (black) [58] are described in this paper. The accelerator neutrino results from T2K (blue) [59], NOvA (purple) [60] and MINOS (light blue) [61] are also shown. This figure is taken from [58].

Tau-Neutrino Appearance
The dominant channel of the atmospheric neutrino oscillation is transition between ν µ and ν τ in the standard scenario. Actually, the long-baseline neutrino beam experiment from CERN to GranSasso, OPERA, which is sensitive to the same neutrino oscillation parameters as atmospheric neutrinos, found the appearance from ν µ to ν τ neutrino oscillation [62]. Finally, OPERA observed 10 ν τ events with a background expectation of 2.0 ± 0.4, which was equivalent to 6.1σ level of discovery significance. The ν τ appearance should also be seen in the atmospheric neutrino data.

Super-Kamiokande
In the atmospheric neutrino observation by the SK detector, the deficit of upward-going ν µ events due to the neutrino oscillation passing through the Earth was observed. The original ν µ is considered to change to ν τ . The appearance of ν τ was also searched by SK [63], but the detection was challenging.
Since the atmospheric neutrino flux falls as 1/E 3 and ν τ charged-current interactions only occur above the τ lepton production threshold, 3.5 GeV, the expected rate at SK is only one event per kiloton per year. Furthermore, τ events are difficult to identify individually as they tend to produce multiple visible particles in the SK detector, as shown in Figure 12.
An analysis was performed that employed a neural network technique to discriminate between "tau-like" and "non-tau-like" events from the hadronic decays of τ detected by SK from atmospheric ν e and ν µ background events. The following seven variables were used as inputs to the neutral network: (1) total visible energy: τ signal events are expected to have higher visible energy compared to the background events; (2) shower-like events: hadronic decay of τ events tend to make a shower in the ring pattern; (3) number of decay electron candidates: pions produced by the hadronic decay of τ produce more decay electrons than the background events; (4) distance between the primary interaction point and decay electron vertex: pions are expected to have smaller momentum compared to the background; (5) sphericity, which is the evaluation whether isotropic or not: hadronic decay of τ is more isotropic than the background; (6) number of Cherenkov ring fragments: τ events are expected to have more ring candidate; (7) ratio of the observed photons and the most-energetic ring in an event: τ events are expected to be small because energy is carried by multiple particles in the hadronic decay of τ. When "tau-like" events are selected from neural network output in this analysis, 76% of signal events and only 28% of background events remain which is estimated by the simulation.  To evaluate whether a τ signal was being observed, the maximum likelihood method was used. The output of a neural network and reconstructed zenith angles were used to construct the probability density function (PDF) for both signals and the background using simulation. As ν τ is generated via neutrino oscillation through the Earth, signals tend to be an upward-going event, which is why the zenith angle was used for the PDF. In this analysis, the 5326 days of atmospheric neutrino data in SK was fitted to the following function: where α is a normalization factor of tau signal, i is a i-th systematic error. According to the normalization factor (α), the failure of tau to appear is zero and the appearance of a tau signal is one. It was found to be 1.47 ± 0.32, assuming the normal mass ordering of neutrino mass splitting, which is equivalent to a 4.6σ level of ν τ appearance. Assuming an inverted mass ordering, it was 1.57 ± 0.31, and the significance was 5.0σ. Figure 13 shows the zenith angle distribution for both tauand non-tau-like events, along with the expectation fitted to tau and the background. These plots show good agreement between the data and MC simulations.
Based on the significant discovery of the τ event, the charged-current tau-neutrino cross-section was also calculated. The measured cross-section was expressed as (1.47 ± 0.32)× < σ theory >, where < σ theory > is the flux-averaged theoretical cross-section. It was calculated by the integral of the differential charged-current ν τ cross-section weighted with the energy spectrum of atmospheric ν τ from neutrino oscillations, and was found to be 0.64 × 10 −38 cm 2 . Finally, the measured cross-section integrated from 3.5 to 70 GeV was calculated as (0.94 ± 0.20) × 10 −38 cm 2 .

IceCube
A search for ν τ appearance by neutrino oscillation using atmospheric neutrino samples by IceCube/DeepCore was performed [64]. The energy region ranges from 5.6 to 56 GeV, which was the same as that in the neutrino oscillation analysis. The identification of individual ν τ events by DeepCore is difficult because of the tau lepton produced via CC interaction decays with a small track length of ∼1 mm, compared to the position resolution of the detector. To search for the appearance of ν τ by neutrino oscillation, the distortion of neutrino energy and direction in cascade-like events compared to non-neutrino oscillation assumption were investigated. In the data analysis, two independent methods were applied: one targeted a high acceptance of all neutrino events, whose background estimation was simulation-driven, and the other was optimized for a higher rejection of non-neutrino events, with data-driven background estimation. Both methods used a boosted decision tree for event selection and background rejection. The neutrino energy and direction were reconstructed by the maximum likelihood method using the charge and time observed by DOMs.
To calculate the flux-averaged theoretical cross section, the differential CC ν τ cross section as a function of neutrino energy is weighted with the energy spectrum of atmospheric tau neutrinos from neutrino oscillations. Because CC ν τ interactions are not distinguishable from CCν τ interactions in Super-K, the theoretical cross section is a flux average of ν τ andν τ cross sections. The flux-averaged theoretical cross section, hσ theory i, is calculated as where dΦð E ν Þ dE ν is the differential flux of tau neutrinos as a function of neutrino energy as shown in Fig. 2, and σð E ν Þ is the differential charged-current tau neutrino cross sections used in NEUT code as seen in Fig. 3. The range of the integral is determined to be between 3.5 and 70 GeV from the tau neutrino energies in the simulation. As shown in Fig. 15, the neutrinos have energies more than 3.5 GeV in the CC ν τ interactions because of the energy threshold, and the expectation of CC ν τ interactions with more than 70 GeV is less than one in the entire run period. The measured cross section is shown togeth theoretical cross sections and the MC sim Fig. 15. The measured and theoretical cross se are consistent at the 1.5σ level. Zenith angle distribution along with the expectation, assuming normal mass ordering, for tau-like (upper) and non-tau-like (lower) events [63]. The fitted tau signal is shown in gray. Figure 14 shows the L/E distribution along with the corresponding predicted counts, given the best-fit neutrino and cosmic-ray muon broken down for the first method of analysis. The plot showed good agreement between the data and the model. The normalization factor of τ appearance in DeepCore was found to be 0.73 +0. 30 −0.24 which is equivalent to a 3.2σ level of ν τ appearance. The result was consistent with those obtained for SK measurements.
The corresponding values for the nuisance parameters can be found in Table II. Figure 17 shows the expected and observed Δχ 2 values for a ν τ normalization ranging from 0 to 2.0. The band of expected values assumes standard oscillations with a ν τ normalization of 1.0. Our main result for the CC þ NC measurement has a best fit value of 0.73 with the 68% confidence interval (C.I.) covering the range (0.49,1.07) and the 90% C.I. covering (0.34,1.30). For the CC-only normalization, we observe the best fit at 0.57 with the 68% C.I. (0.30,0.98) and the 90% C.I. (0.11,1.25).
These measured values are compatible with corresponding values obtained from analysis B within less than 1σ standard deviation. These confirmatory results of analysis B are 0.59 þ0. 31 −0.25 (0.43 þ0.36 −0.31 ) for the CC þ NC (CC-only) measurement, also see Fig. 18.
All values are also compatible within the 90% confidence interval with expectations assuming the three-flavor neutrino oscillation paradigm (i.e., ν τ normalization ¼ 1.0) and the assumed ν τ CC cross sections. The significance at which we can reject the null hypothesis of no ν τ appearance is 3.2σ and 2.1σ for the CC þ NC and the CC-only case for analysis A, respectively. The confirmatory analysis B yields slightly weaker limits of 2.5σ (1.4σ).
The confidence intervals for the measurements presented here, shown in Fig. 18, are calculated using the approach of Feldman and Cousins [69] to ensure proper coverage.
The presented results are of a comparable precision to those of SK and OPERA (see Fig. 18), and complementary to those measurements in terms of energy scale, L=E range, systematic uncertainties, and statistics. Specifically, the SK measurement is based on lower-energy events where   Figure 14. L/E distribution along with the expectation of the best-fit neutrino and cosmic-ray muon [64]. The bottom portion shows the ratio of the data and the expectation.

Sterile Neutrino Analysis
The existence of neutrino oscillation has been established by a wide range of experiments using different sources of neutrinos: the atmosphere, the Sun, nuclear reactors, and accelerators. However, not all neutrino experiments show the three standard flavors of neutrinos. For example, an excess of electron neutrinos in a muon-neutrino beam with ∆m 2 ∼ 1eV 2 was found in the LSND [65] and MineBooNE [66] experiments. Additional anomalies appeared in theν e and ν e rates in several reactor experiments [67] and gallium-based experiments [68]. These results made a hint of neutrino oscillations driven by a ∆m 2 > 1eV 2 . On the other hand, the number of neutrinos was measured as 2.980 ± 0.0082 light neutrino flavors by large electron-positron(LEP) collider using the width of the Z 0 mass peak [69]. This result required that its mass of an additional neutrino is either heavier than half the Z 0 boson; however, it is difficult to be a player of neutrino oscillations, or not interact via weak interactions, which is called sterile neutrinos.
To introduce N additional sterile neutrinos, U in Equation (1) should be (3 + N) × (3 + N) matrix. The matter effect in neutrino oscillation for sterile neutrinos should also be considered. Since sterile neutrinos have no CC nor NC interactions, the effective matter potential is expressed by V n = ±G F / √ 2N n . It is derived from the difference from ν µ and ν τ which have only NC. Here, NC depend only on neutron density (N n ) because the interactions with electrons and protons are equal and opposite, and their densities are identical in neutral matter.
Atmospheric neutrino samples can also make a constraint on sterile neutrinos due to a wide range in both energy and travel distance. The results of single additional sterile neutrino scenario, called the "3+1" model, have been reported by SK [70], IceCube [71], and ANTARES [58]. Actually, it hard for this model to explain all the anomalies consistently; however, it can be extended to models with more than one sterile neutrino although it needs more parameters. The "3+1" model, denoted as U, is that a neutrino flavor eigenstate ν s with mass eigenstate ν 4 should be added in the mixing matrix defined in Equation (2). In this model, six new parameters are added: three mixing angles θ 14 , θ 24 , θ 34 , two CP-violating phases, and one mass differences ∆m 2 41 To search for sterile neutrinos using atmospheric neutrinos, several assumptions were made. First, since the assumed mass difference between the mass eigenstate ν 4 and the other is large, sin 2 (∆m 2 41 L/4E) is averaged as 1/2. Second, the other parameters, except for |U µ4 | 2 = sin 2 θ 24 , |U τ4 | 2 = cos 2 θ 24 sin 2 θ 34 , (δ 24 only for ANTARES), are ignored since they have negligible impact on atmospheric neutrino oscillation. The effect of these neutrino oscillation parameters results in a different ν µ disappearance pattern compared to the standard neutrino oscillation scheme.
The analysis method is similar to the standard three-neutrino hypothesis, i.e., fitting the reconstructed energy and cosine zenith angle from the experimental atmospheric neutrino data to the prediction but using only the ν µ disappearance channel. No clear evidence of sterile neutrinos was found in these experiments. Figure 15 shows the exclusion region at 90% and 99% C.L. in the |U µ4 | 2 and |U τ4 | 2 parameters.  Figure 15. Contours in 2D exclusion regions of |U µ4 | 2 = sin 2 θ 24 and |U τ4 | 2 = cos 2 θ 24 sin 2 θ 34 at 90% (left) and 99% (right) C.L. obtained by the atmospheric neutrino data in SK (blue) [70], IceCube/DeepCore (red) [71] and ANTARES (black) [58]. The dashed lines show the normal mass ordering, while the solid lines show the inverted mass ordering (IceCube/DeepCore) and unconstraint δ 24 (ANTARES). The markers show the best-fit values for each experiment. This figure is taken from [58].

Atmospheric Neutrino Flux Measurements
Predictions and observations of the atmospheric neutrino flux have been reported by several groups. Its energy distribution is derived from charged particles, mainly muons or electrons, induced by neutrino interaction. Figure 16 shows several predictions and measurements of the energy distribution in atmospheric neutrinos. In this figure, SK applies the so-called "unfolding" method to all experimental phases [72]. The event rate of the observed charged particles induced by atmospheric neutrinos is expressed as the convolution of neutrino flux, neutrino oscillation probability, neutrino cross-section, and detector efficiency. The obtained atmospheric neutrino spectrum is derived from the deconvolution of the above observed values, which is the unfolding method. The observed flux of atmospheric ν µ and ν e is indicated by the red square and blue circle in Figure 16, and are consistent with the predictions that include an assumption of neutrino oscillation. The χ 2 values (p-value), including ν e and ν µ together, are 22.2 (0.51) for HKKM, 30.7 (0.13) for Bartol, and 25.6 (0.32) for FLUKA, with a degree of freedom of 23.
IceCube reported the results of atmospheric ν µ and ν e flux. In the ν µ analysis, the observable value is the deposit energy per track length of muons in the detector (dE/dX). It is a convolution of the atmospheric neutrino flux and the response matrix, which accounts for the effect of propagation through the Earth, neutrino interaction, detector response, and event selection. Unfolding of the atmospheric muon-neutrino flux from 100 GeV to 400 TeV was performed [73], and the results are indicated by the pink triangle in Figure 16. In addition, the so-called forward-folding analysis was applied, in which the dE/dX distribution was tested against the hypotheses of muons arising from atmospheric ν µ by π or K decay, prompt ν µ , and astrophysical ν µ [53]. At first, no evidence was found for an astrophysical or prompt ν µ . In a comparison with HKKM, the normalization factor of the absolute atmospheric neutrino flux was found to be 0.96 ± 0.16 and the spectra index was found to be steeper by E −0.032±0.014 . The allowed regions of these parameters from 332 GeV to 84 TeV are indicated by the pink band in Figure 16. As for the atmospheric ν e measurements by IceCube, two results were reported: one was a low-energy region (80 GeV to 6 TeV) using DeepCore [74] and the other was a high-energy region (0.1 to 100 TeV) using the full IceCube detector [75]. After several steps of background reduction, related to ν µ in the main sample and cosmic-ray muons from the cascade sample, the data were fitted to the prediction to obtain the atmospheric ν e flux. The observed flux was consistent with the prediction, as shown in Figure 16. In addition, no prompt neutrino signal was found for ν e .
ANTARES reported the atmospheric ν µ energy spectrum in the energy range between 0.1 and 200 TeV [76]. It was derived from the measured muon energy distribution through a response matrix, determined from simulations, and an unfolding method. The results are shown in the figure, although it is admittedly busy. The overall normalization factor is 25% higher than the prediction in [21], and the flux is compatible with the IceCube results. The prompt neutrino was not observed in the ANTARES data.
An east-west flux asymmetry of atmospheric neutrinos due to the rigidity cutoff is predicted. SK reported the measurement of this effect [72] using an FC sample, selecting both electron-and muon-like events with a single reconstructed Cherenkov ring. Figure 17 shows the azimuthal distribution of the subsample event selected to optimize the significance of the east-west dipole asymmetry. The effect is observed at a significance level of 6.0 (8.0) sigma for the muon-like (electron-like) samples. The dependence of the asymmetry on the zenith angle was also investigated and was observed at the 2.2 sigma level. These effects were consistent with the prediction within the given uncertainties.
SK also reported long-term and seasonal atmospheric neutrino flux variations [72]. An anticorrelation with solar activity was predicted for the long-term variation and a weak preference for a correlation was observed at the 1.1 sigma level using 20 years' observation. The seasonal variation is considered to be occurred by the change in the atmospheric density profile over the year; however, no such correlation was observed in SK. data sets.
Although we thoroughly validated the accuracy of the unfolding procedure using the HKKM-like pseudodata and included the estimated regularization bias as a systematic uncertainty, there has still been some concern about the ability of the unfolding procedure to accurately reproduce more complicated spectral shapes, such as the wavy shape that was eventually obtained, and, in particular, whether or not such a shape would be more strongly affected by the imately within the estimated regularization uncertainties. Therefore, we conclude that the shapes of our unfolded spectra are not due to an unexpected additional bias from the unfolding procedure.
In Fig. 14 this flux measurement is shown with the results from other experiments. Our measured data provide significantly improved precision below 100 GeV. The minimum of the observed energy range is extended below 1 GeV, and at higher energies overlaps with ν μ FIG. 14. The measured energy spectra of the atmospheric ν e and ν μ fluxes by SK, shown with measurements by other experiments, i.e., Frejus [39], AMANDA-II [40,41], IceCube [42][43][44][45], and ANTARES [75]. The phrase "forward folding" used by IceCube and AMANDA-II is synonymous with forward fitting. The HKKM11 flux model predictions for the Kamioka site are also shown in solid (with oscillation) and dashed (without oscillation) lines. The error bars on the SK measurement include all statistical and systematic uncertainties. Figure 16. Energy spectra of the atmospheric ν e and ν µ fluxes according to several experiments with an overlay of HKKM predictions for the Kamioka site in solid (with neutrino oscillation) and dashed (without oscillation) lines. SK used the unfolding method [72]. Both forward-folding and unfolding methods were reported by IceCube in the energy range of 100 GeV to 400 TeV for ν µ and 100 GeV to 100 TeV for ν e [73][74][75]. ANTARES used the unfolding method in the energy range of 100 GeV to 200 TeV [76]. This figure is taken from [72]. a larger decrease in flux. However, if two NMs are located at different rigidity cutoffs, it is assumed that while the gradient of the correlation will be different, a linear . The subselection is optimized to obtain the highest significance of the final A parameters, by using only events with 0.4 < E rec < 3.0 (GeV) and j cos θ rec j < 0.6.

MEASUREMENTS OF THE ATMOSPHERIC NEUTRINO FLUX …
PHYSICAL REVIEW D Figure 17. Azimuthal distribution of a subsample of electron-like (top) and muon-like (bottom) events from the data with statistical error [72]. The boxes are prediction with systematic error. A subsample is selected as it is optimized to obtain the highest significance by using only events with 0.4 < E < 3.0 GeV and | cos θ zenith | < 0.6.

Future Prospects
In this section, the search for atmospheric neutrinos by next-generation neutrino detectors is briefly described. Table 2 shows the summary of the detectors. Succession to the current water Cherenkov detector experiments are now proposed. Here, an enlargement of the volume in underground detector and denser photo-sensor in string type detector are crucial to determine unknown neutrino oscillation parameters such as mass ordering and CP-violating phase.
Hyper-Kamiokande (HK) [77] is underground (1750 m water equivalent) water Cherenkov detector, located 8 km from the SK site. It is based on the well-established technology with a fiducial volume of 187 kilotons, which is 8.3 times larger than that obtained for SK, and will start operation in 2027. It consists of 40,000 high-quantum-efficiency PMTs for inner detector. The second HK detectors in future is also proposed [78].
The IceCube plans detector upgrade as "IceCube-Gen2" [79]. The first step of the upgrade is scheduled for deployment in the 2022/2023 polar season [80]. It is proposed that 7 additional strings with 125 optical modules spaced 2.4 m apart, which is denser than IceCube/DeepCore, will set in a small part of IceCube observatory (2 MegaTon of ice). A denser detector, PINGU [81], that consists of 26 strings with 192 optical modules spaced 1.5 m apart, is proposed as a goal of IceCube-Gen2. The denser photo-sensor enables the neutrino oscillation parameter determination with high precision [82].
The successor of the ANTARES is KM3NeT [83]; it is a deep-sea neutrino detector in the Mediterranean Sea. There are two installation sites for different purpose; one is "ARCA" at 3500 m depth in Italy which aims an observation of high-energy neutrino sources in the Universe, the other is "ORCA" at 2500 m depth in France which aims a neutrino oscillation parameter determination using atmospheric neutrinos. The multi-PMT is set in the DOM as an optical sensor, and 18 DOMs (spaced 9 m apart for ORCA) are hosted in one string that is called "Detection Units (DU)" The detector construction was started, and its performance has been checked [84]. The first DU for ORCA had already been installed, and 5 additional DUs will be installed soon, and finally 115 DUs will be installed.
Apart from water Cherenkov technique, several different types of detector available for atmospheric neutrinos are proposed. DUNE [85] is a long-baseline neutrino experiment using a massive liquid argon time-projection chamber (LAr TPC). The far detector is located at a depth of 4300 m water equivalent at the Sanford Underground Research Facility, in Lead, SD, USA. One of the advantages of LAr TPC is excellent particle identification and better angular and energy resolutions than water Cherenkov detectors. The well reconstructed momentum of each final state particles enables the determination of the energy and direction of the incident neutrinos. The fiducial volume of one LAr TPC module is 10 kilotons. The first module construction has started, and 4 modules are the goal. The DUNE operation together with the neutrino beam will start in 2026; however, the atmospheric neutrino observation will be possible to start as soon as the first module completes when it is expected before 2026.
Main topics of the atmospheric neutrino observation in neutrino oscillations is determination of mass ordering. Any future detectors described here will reach more than 3σ level to reject the incorrect mass ordering, assuming that either the normal or inverted mass ordering is true, after 5-10 years of operation, although this depends on other parameters of the neutrino oscillation. The precise tau appearance observation and studies of any exotic model are also possible.

Summary
The research into atmospheric neutrinos has had great success in the last two decades, notably the discovery of neutrino oscillation by the SK detector. This paper reviews the recent achievements made in atmospheric neutrino observations by SK, IceCube, and ANTARES. First, the standard three-neutrino oscillation scheme, expressed by the PMNS matrix, was established and the parameters were determined. The appearance of tau neutrinos in atmospheric neutrinos, predicted by this theory, was confirmed. An exotic scenario, such as the appearance of sterile neutrinos, was tested; however, no clear signal was observed. The observed atmospheric neutrino flux was consistent with the predictions. Among the unknown neutrino oscillation parameters, the normal mass ordering is preferred in the current measurements; however, it is not significantly determined. Atmospheric neutrino observations in the next generation of neutrino detectors are expected to provide final determination of the mass ordering.
Funding: This research received no external funding.