Neutrino oscillations in the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) paradigm are parametrized by three mixing angles, two squared-mass differences, and a CP-violating phase [1,2]. Experiments measuring neutrinos of different flavors, energies, and baselines have constrained many of the PMNS parameters with increasing levels of precision. However, the octant of the mixing angle θ 23 , the phase δ CP , and the sign of the larger of the two squared-mass differences, Δm 2 32 , which determines the neutrino mass ordering, are all presently unknown. To date, the long-baseline accelerator neutrino experiments T2K [3] and NOvA [4] have made the world's most precise measurements of θ 23 , Δm 2 32 , and δ CP , but they have yet to definitively resolve the remaining questions.
Atmospheric neutrinos are an independent and natural counterpart to accelerator neutrinos for studying neutrino oscillations. Neutrinos created in the Earth's atmosphere span a range of energies and baselines that make their oscillations sensitive to the θ 23 mixing angle and the magnitude of the Δm 2 32 squared-mass difference. Additionally, atmospheric neutrinos which pass near or through the dense core of the Earth experience matter effects which alter their oscillation probabilities. An observation of these modified oscillation probabilities in either atmospheric neutrino or antineutrino data would provide important information toward resolving the neutrino mass ordering.
In this work, we analyze 6511.3 live days of atmospheric neutrino data from the Super-Kamiokande (SK) detector. This analysis improves upon the previous work [5] in three major ways: We use the number of tagged neutrons to enhance the separation of neutrino events from antineutrino events, we enhance the efficiency of classifying multi-ring events using a boosted decision tree (BDT), and we add 48% exposure by analyzing events from an expanded fiducial volume and from 1186 additional live days, including data collected after a major detector refurbishment in 2018. In addition to the atmospheric-only analysis, we present an analyses of SK data with an external constraint on the mixing angle θ 13 from the average measurements of the reactor neutrino experiments Daya Bay [6], RENO [7], and Double-Chooz [8].
The paper is organized as follows: Section I presents an overview of neutrino oscillation phenomenology relevant to atmospheric neutrino oscillations. Section II provides a description of the Super-Kamiokande detector and its capabilities for reconstructing neutrino interactions. Section III describes the simulation used to model atmospheric neutrinos interactions at SK. Section IV describes the analysis methodology and presents the results of the analyses without external constraints and with constraints on sin 2 θ 13 . We provide an interpretation and summary of the results in Sec. V.
Neutrinos are produced as flavor eigenstates of the weak interaction, which may be treated as superpositions of mass eigenstates via the PMNS matrix:
where α is a label for each lepton flavor-one of e, μ, or τand U αi is an element of the PMNS matrix. The PMNS matrix is parametrized by three mixing angles and a phase, and it factorizes into three submatrices which describe rotations by each mixing angle from the neutrino mass basis into the neutrino flavor basis:
0 0 0 c 23 s 23 0 -s 23 c 23 1 C A 0 B @ c 13 0 s 13 e -iδ CP 0 1 0 -s 13 e iδ CP 0 c 13 1 C A × 0 B @ c 12 s 12 0 -s 12 c 12 0 0 0 1 1 C A:
In Eq. ( 2), the sines and cosines of the mixing angles are written as cos θ ij ≡ c ij and sin θ ij ≡ s ij , respectively, and the phase δ CP changes sign for the antineutrino case. The probability of a neutrino of one flavor jν α i oscillating to a different flavor jν β i after some time t-or, equivalently, along a baseline L-is found by computing the amplitude jhν β jν α ij 2 . The probability is nonzero for the case α ≠ β if the mass states have nonzero mass differences given by the signed quantity Δm 2 ij ¼ m 2 i -m 2 j . In the simplest case, neutrinos oscillate in a vacuum, and oscillation probabilities may be computed by propagating neutrino states according to their vacuum Hamiltonian, written here in the mass basis:
This leads to oscillation probabilities of the form
where Δ ij ¼ 1.27Δm 2 ij L=E. Here, Δm 2 ij is expressed in units of eV 2 , L is the oscillation baseline in kilometers, and E is the neutrino energy in GeV. Experiments have measured all mixing angles and squared-mass differences to be significantly different from zero, while the value of the phase δ CP is still unknown. 1 In addition, solar neutrino oscillation experiments observe evidence for matter effects in the Sun which imply that the Δm 2 21 squared-mass difference is positive, establishing an ordering for two of the neutrino masses, m 2 > m 1 [9][10][11][12]. However, current experiments are consistent with either the normal ordering, m 3 ≫ m 2 ; m 1 , or the inverted ordering, m 2 ; m 1 ≫ m 3 . Consequently, the sign of the squared-mass difference between m 3 and the next-most-massive neutrino, given by either Δm 2 32 or Δm 2 31 , is not known. We use the notation Δm 2 32;31 or simply Δm 2 for this squared-mass difference where the ordering is not explicitly specified.
Numerically, Δm 2 32;31 has been measured to be approximately 30 times larger than Δm 2 21 , such that Δm 2 32 ≈ Δm 2 31 . The difference in magnitude between Δm 2 21 and Δm 2 32;31 also implies that the terms in Eq. ( 4) containing one or the other squared-mass differences dominate for different ranges of L=E. For long-baseline beam and atmospheric neutrinos, where neutrino baselines range from tens of kilometers to several thousand kilometers, and typical neutrino energies range from MeV to several GeV, the Δm 2 21 terms are subdominant, leading to approximate flavor oscillation probabilities of the form
sin 2 θ 23 sin 2 2θ 13 sin 2 1.27 Δm 2 L E ; Pðν μ → ν μ Þ ≈ 1 -4 cos 2 θ 13 sin 2 θ 23 ð1cos 2 θ 13 sin 2 θ 23 Þ × sin 2 1.27 Δm 2 L E ; Pðν e → ν e Þ ≈ 1sin 2 2θ 13 sin 2 1.27 Δm 2 L E :
The approximate oscillation probabilities in Eq. ( 5) are primarily functions of the mixing angles and the absolute value of the squared-mass difference Δm 2 . The phase δ CP , and the neutrino mass ordering-i.e., the sign of the squared-mass difference-are subleading effects which make them challenging experimental signatures.
Neutrino oscillations in matter enhance the dependence of oscillation probabilities on the neutrino mass ordering. In matter, due to an increased forward scattering amplitude, electron-flavor neutrinos experience a larger potential relative to μ and τ flavors, which modifies the vacuum Hamiltonian via an additional term,
Here, G F is the Fermi constant, N e is the electron density, and U is the PMNS matrix. The sign of a is positive for neutrinos and negative for antineutrinos. Propagating the neutrino states according to the matter Hamiltonian leads to an effective squared-mass difference and mixing angle:
where Γ ≡ 2aE=Δm 2 . Equation (7) shows that the effective quantities depend on the sign of Δm 2 . In particular, for neutrinos in the normal ordering, Γ ≈ cos 2θ 13 maximizes the effective mixing angle sin 2 θ 13;M . A maximum also occurs for antineutrinos in the inverted ordering. This maximum effective mixing angle predicts a resonant enhancement of muon-to-electron flavor conversions for either neutrinos or antineutrinos according to the neutrino mass ordering.
Atmospheric neutrinos are produced when cosmic rays interact with nuclei in the Earth's atmosphere. These interactions result in hadronic showers of primarily pions and kaons, which decay into neutrinos. The atmospheric neutrino energy spectrum extends from a few MeV to several TeV and has an approximate flavor ratio in the few-GeV range of ðν μ þ νμ Þ=ðν e þ νe Þ ≈ 2∶1. While present, tau neutrinos intrinsic to the atmospheric neutrino flux are suppressed by many orders of magnitude relative to electron-and muon-flavor neutrinos due to kinematic restrictions on their production.
The zenith angle θ z describes atmospheric neutrino baselines. Neutrinos produced directly above a detector are downward-going, θ z ¼ 0, and are produced at an average distance of 15 km above Earth's surface. Neutrinos produced on the other side of the Earth from a detector are upward-going, θ z ¼ π, and travel an approximate distance of 13 000 km through the Earth. Oscillation signatures are most evident in upward-going atmospheric neutrinos due to the longer baselines.
A general atmospheric neutrino baseline begins at a production point in the atmosphere and passes through the Earth before ending at a detector near the surface. We model the matter effects induced by passage through the Earth assuming a simplified version of the preliminary reference Earth model (PREM) [13], where the Earth is treated as a sphere with radius R Earth ¼ 6371 km and contains concentric spherical shells of decreasing densities. Table I lists the Earth layers and corresponding densities assumed in this work.
To compute neutrino oscillation amplitudes through layers of different matter densities, amplitudes along steps through matter of fixed densities are multiplied together [14]. The general matrix form of the propagated mass eigenvectors X for neutrinos passing through a fixed matter density is
TABLE I. Neutrino propagation layers and corresponding densities used for calculating neutrino oscillation probabilities in this analysis, based on a simplified PREM [13]. Layer R Min (km) R Max (km) Density (g=cm 3 ) Atmosphere 6371 Á Á Á 0 Crust 5701 6371 3.3 Mantle 3480 5701 5.0 Outer core 1220 3480 11.3 Inner core 0 1220 13.0 where Mfoot_1 i =2E are the eigenvalues of H Matter . This definition allows the neutrino probability along a baseline of changing matter density to be written as
where L i and ρ i are the baseline and density of the ith step, respectively. In the case of a spherically symmetric Earth, as is assumed in this work, the neutrino oscillation baseline L only depends on the zenith angle and production height and does not depend on the azimuth. Figure 1 shows the calculated oscillation probabilities for atmospheric neutrinos as a function of the cosine of the zenith angle and the neutrino energy E ν . The two sets of figures show the cases for muon-to electron-flavor neutrino and antineutrino oscillation probabilities in each massordering scenario. The resonance in ν μ → ν e or νμ → νe oscillations due to matter effects is exclusively visible in the normal and inverted scenarios, respectively, and occurs at baselines of several thousand kilometers and neutrino energies around a few GeV. This resonance in atmospheric neutrinos is the experimental signature of the unknown neutrino mass ordering.
In addition to mass-ordering sensitivity via Earth matter effects, atmospheric neutrino oscillations also provide sensitivity to other oscillation parameters. Muon-to-tau flavor conversions provide sensitivity to jΔm 2 32;31 j and sin 2 θ 23 . While the tau neutrinos are often too low-energy to produce CC interactions, the Pðν μ → ν μ Þ survival probability manifests as a disappearance of upward-going muon neutrinos. Atmospheric neutrinos also provide modest sensitivity to the combined effects of δ CP and sin 2 θ 13 through electron neutrino or antineutrino appearance for neutrinos of all energies.
Super-Kamiokande (SK) is a 50 kiloton cylindrical water-Cherenkov detector located within the Kamioka mine in Gifu, Japan [15,16]. The detector consists of two optically separated regions: an inner detector (ID), which contains 32 kilotons of water and is viewed by over 11 000 inward-facing 20-inch photomultiplier tubes (PMTs), and a 2 m-thick outer detector (OD) with over 1800 outward-facing 8-inch PMTs for vetoing cosmic backgrounds. To increase the light collected in the OD, the OD walls are covered with reflective Tyvek, and wavelength-shifting plates are mounted to the OD PMTs. The ID PMTs have typical quantum efficiencies of 20% for wavelengths near 400 nm and timing resolutions of 2 ns for 1 photoelectron [17].
SK has been operational since its construction in 1996 and, until July 2020, has operated with pure water. During this period, there were five distinct data-taking phases, SK I-V, which had similar operating conditions with a few notable exceptions: At the end of the SK I phase (1996)(1997)(1998)(1999)(2000)(2001), an accident 2 resulted in the loss of roughly half of the experiment's PMTs. During the SK II phase (2002)(2003)(2004)(2005), the remaining PMTs were rearranged to provide uniform but reduced (19%) photocoverage. New PMTs were installed to restore photocoverage to original conditions starting with the SK III phase (2006)(2007)(2008). In 2008, the experiment's electronics were upgraded [18], marking the start of the SK IV phase (2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018). The SK IV electronics upgrade extended the window of hit times recorded following neutrino-like events, which improved the efficiency for detecting decay electrons and enabled the detection of the 2.2 MeV gamma emission following neutron captures on hydrogen [19].
During 2018, the SK detector was drained to conduct work in preparation of loading gadolinium sulfate, a compound with a high neutron capture cross section, into the detector's water. The work comprised of installing a new water circulation system capable of continuously purifying gadolinium sulfate, sealing the welding joints of the tank walls to repair and prevent leaks, and replacing several hundred PMTs that had failed during the SK III and SK IV phases. The subsequent SK V (2019-2020) phase resumed data-taking with pure water and reaffirmed the detector's stability and performance after the refurbishment work. In July 2020, gadolinium was dissolved into the detector's water for the first time [20], marking the start of the SK Gd phase. This work includes data from the purewater SK I-V phases only. Future analyses using data from the SK Gd phase will feature enhanced neutron-tagging efficiency due to the high-neutron-capture cross section of the gadolinium and its subsequent 8 MeV γ cascade. The operating conditions of the SK phases are summarized in Table II.
The SK detector observes Cherenkov light from charged particles with sufficient momentum produced following neutrino interactions. The light projected onto the PMTlined walls of the detector forms ring patterns of hit PMTs. The ring patterns are reconstructed using the APFIT [21] algorithm: The timing information of hits establishes an event vertex, and a fit considering the spatial distribution and observed charges of hit PMTs determines the particle's direction, momentum, and particle type. APFIT separates rings into e-like and μ-like: Electrons and photons tend to scatter and produce electromagnetic showers, resulting in many overlapping rings which appear as a single ring with blurred edges. Heavier particles such as muons and charged pions do not create showers, and therefore have sharp ring edges. Higher-momentum particles produce more light, so the charge contained within the ring provides an estimate of the particle's momentum. Figure 2 shows an example of the Cherenkov light patterns observed in SK following a neutrino candidate interaction, and their fitted properties. The event contains multiple ring patterns, each corresponding to a different particle.
In addition to Cherenkov rings, SK identifies electrons from muon decays. Decay electrons are found by scanning for time-clustered hits following a primary neutrino interaction trigger. A hit-time-based fitter estimates the decay electron vertex for each candidate hit cluster, and candidates are accepted if there are 50 or more hits within a 50 ns time window. The overall decay electron tagging efficiency is estimated to be 96% for μ þ and 80% for μ -in the SK IV and SK V periods. The reduced efficiency for μ -is due to μ -capture in the water, in which no decay electron is produced.
Neutrons in the SK detector are captured on hydrogen, producing deuterium in an excited state. The decay of the excited deuterium produces a 2.2 MeV γ, which results in a few time-coincident and spatially clustered PMT hits. These γ emissions from neutron captures are identified using a two-step process. In the first step, a sliding 10 ns hit-time window finds candidate neutron captures from clusters of 7-50 hits with fewer than 200 hits in a surrounding 200 ns window. The lower bound of the hit range suppresses spurious coincidences from noise, while the upper bound avoids tagging decay electrons. In the second step, variables which quantify the isotropy, likelihood of single-vertex origin, and the time spread of the hits are calculated for each candidate cluster. A neural network classifies candidates as either signal or background based on these variables. When applied to SK IV-V atmospheric neutrino MC events, the neural network has an average neutron-tagging efficiency for neutron capture on hydrogen of 26% with a background rate of 0.016 false neutrons per event. The uncertainty on the neutron-tagging efficiency is evaluated using an americium-beryllium (AmBe) source embedded in a scintillating box placed at various locations throughout the detector, and is estimated to be 9%. A detailed description of the neutrontagging algorithm and its development may be found in Ref. [22].
In the SK detector, the charges of particles-and therefore neutrino and antineutrino interactions-cannot be differentiated on an event-by-event basis. However, statistical separation is possible. For example, in the process
in which an antineutrino interacts with a proton, the outgoing negatively charged pion is more likely to be captured by an 16 O nucleus before decaying than is a positively charged pion produced in the equivalent ν μ interaction. Captured pions do not produce decay electrons, so requiring one or more decay electrons preferentially selects more neutrino than antineutrino events for this process. The statistical separation can be further improved by also considering the number of neutrons, which will be described in Sec. II B 2.
Calibration ensures an accurate and consistent response of the detector to particle interactions. Calibration studies
TABLE II. Super-Kamiokande data-taking phases. An electronics upgrade at the start of SK IV enabled neutron tagging on hydrogen (H), utilized in the SK IV and SK V phases. During 2020, gadolinium (Gd) was added to the detector's water to increase the neutron-tagging efficiency. At the time of this writing, SK Gd is ongoing, and data from the SK Gd phase are not included in this analysis. Phase Dates Live time (Days) Photo-coverage (%) Neutron tagging SK I 1996-2001 1489.2 40 Á Á Á SK II 2002-2005 798.6 19 Á Á Á SK III 2006-2008 518.1 40 Á Á Á SK IV 2008-2018 3244.4 40 H SK V 2019-2020 461.0 40 H SK Gd 2020-Present Á Á Á 40 H þ Gd μ-like e-like FIG. 2. Event display of a multi-ring atmospheric neutrino cadidate event in SK V data. Hit ID PMTs are displayed on an unrolled view of the cylindrical detector, with the color and radius of each hit PMT corresponding to the detected charge. The reconstruction algorithm APFIT identifies three Cherenkov rings, indicated by dashed outlines: one bright μ-like ring, with p μ ≈ 1010 MeV=c, and two fainter e-like rings, each with p e ≈ 320 MeV=c.
based on the detector geometry, PMT responses, and properties of the SK water are documented in Ref. [16]. We assess APFIT's energy determination using calibration sources at multiple energies, which span as much of the atmospheric neutrino energy spectrum as possible. Cosmic ray muons which stop within the detector and produce a decay electron provide a way to measure photoproduction as a function of the muon's track length. The track length of these muons is determined by their ID entrance point and decay electron vertex. The expected momentum of each muon, assuming minimum ionization along the track length, may be compared to the fitted Cherenkov ring momentum. The energy spectrum of decay electrons from these muons also provides a low-energy calibration source. We use the momenta of the two e-like rings from neutral pion decay, π 0 → γγ, which add to form an invariant mass distribution of neutral pion events, as a third calibration source.
Figure 3 shows the difference between data and Monte Carlo (MC) for the different calibration sources during each SK phase. In this analysis, the decay electron and π 0 mass measurements were also performed separately using events with vertices greater than 200 cm (conventional fiducial volume) and between 100 cm and 200 cm (additional fiducial volume; see Sec. II B 4) from the detector walls. The largest difference, 4%, is observed for the decay electron calibration events in the additional fiducial volume region during the SK II phase.
Neutrino events at SK are broadly categorized as fully contained (FC), partially contained (PC), or upward-going muons (Up-μ). FC and PC events have a reconstructed event vertex within the ID and are differentiated based on the number of hits detected in the OD: FC events have minimal OD activity, while PC events have OD activity following the primary event trigger. Up-μ events are neutrino interactions within the rock below the SK tank or in the OD water which produce muons traveling upward. Figure 4 shows the average number of FC, PC, and Up-μ events observed per day during the SK I-V phases. The average event rates for each category are consistent across all phases within statistical uncertainties.
We separate neutrino candidate interactions from each category into analysis samples to enhance the different oscillation signals present in the atmospheric neutrino data. The data taken during the SK IV and SK V periods, 3705 days, or 57% of the total SK I-V exposure, uses the observed number of neutron captures on hydrogen as an additional classification handle to enhance the purity of neutrino and antineutrino samples.
Fully contained events span the energy range of 100 MeV to 100 GeV and are a mixture of charged-current (CC) and neutral-current (NC) neutrino interactions of all flavors. Because the normal and inverted neutrino mass
6 4 2 0 2 4 6 SK I SK II 2 10 3 10 4 10 SK III 6 4 2 0 2 4 6 2 10 3 10 4 10 SK IV 2 10 3 10 4 10 SK V Conventional FV Michel-e Expanded FV Michel-e Mass 0 S Conventional FV Mass 0 S Expanded FV Sub-GeV Stopping Muon Multi-GeV Stopping Muon Momentum (MeV/c) Momentum (MeV/c) (MC -Data) / Data (%)
FIG. 3. Performance of energy reconstruction for various sources across all SK periods. Solid points indicate measurements made within the conventional fiducial volume (vertex is greater than 200 cm from the detector walls), while dashed points correspond to measurements in the additional fiducial volume (vertex is between 100 cm and 200 cm from the detector walls). Michel-e refers to the stopping muon decay electron spectrum calibration.
ordering scenarios predict an enhancement to the number of either ν μ → ν e or νμ → νe events, and SK cannot distinguish the sign of neutrino interactions on an event-by-event basis, the FC sample definitions are designed to increase the statistical purity of ν e and νe events.
Fully contained events with a single Cherenkov ring are first separated by the ring's particle identification (PID) score, either e-like or μ-like. Events are further divided based on the visible energy E vis into sub-GeV, E vis < 1330 MeV, and multi-GeV, E vis > 1330 MeV. Next, events are separated by the number of decay electrons. For sub-GeVevents, we use the number of decay electrons to separate events by likely interaction processes. Samples enhanced with quasielastic interactions are formed by requiring no decay electrons for e-like events and either zero or one decay electron for μ-like events. Sub-GeV e-like and μ-like events with one or more, or two or more decay electrons, respectively, are separated into additional samples enhanced in interaction processes which produce a pion below Cherenkov threshold. For multi-GeV events, the number of decay electrons, either zero, or one or more, is used to separate e-like events into antineutrino-and neutrino-enhanced samples. Fully contained single-ring events in the SK IVand SK V phases have been revised in this analysis to additionally incorporate the number of tagged neutrons and are discussed in Sec. II B 2.
Multi-ring events can contain mixtures of e-like and μ-like rings, making the neutrino flavor ambiguous. However, multi-ring events also provide extra information which is useful for separating ν e events from νe events. For multi-GeV multi-ring events, we use a boosted decision tree (BDT) to classify these events as ν e -like, νe -like, μ-like, or "other." The "other" sample primarily selects NC events. More details on the BDT are presented in Sec. II B 3.
Sub-GeV multi-ring events are not included in the analysis due to their poor direction resolution and consequently minimal sensitivity to oscillation effects, with two exceptions: First, multi-ring events with E vis > 600 MeV, where the most energetic ring is μ-like and has a reconstructed momentum of at least 600 MeV=c, are included in the multi-ring μ-like sample. Second, FC NC interactions which produce a π 0 are an oscillation-insensitive background to the other samples, but they are included in the analysis to constrain NC interactions. These neutral-current π 0 events are identified from sub-GeV events using a dedicated fitter which assumes there are two rings present, regardless of the number of reconstructed rings. Events are classified as π 0 -like based on the likelihood that the two fitted rings originate from a π 0 decay. Events which are classified as π 0 -like are separated into two samples based on the number of reconstructed rings without the two-ring assumption, either one or two.
Partially contained events have typical energies between 1 GeV and 1 TeV, and they are nearly all ν μ CC interactions, as the muon produced in the interaction often exits the ID. The muon momentum in PC events can only be estimated using the portion of the track within the detector. Partially contained events which exit the ID and stop within the OD, determined by comparing the amount of light in the OD to simulated PC events, are classified as stopping, while PC events that completely exit the detector are classified as through-going.
Upward-going muon events are the highest-energy events, up to ∼10 TeV, observed at SK and are classified into three samples: "stopping" if the muon stops within the ID, and otherwise "showering" or "nonshowering," depending on whether or not the exiting muon's charge deposition is consistent with radiative losses. The Earth shields the Up-μ sample from cosmic ray backgrounds from below the horizon. However, upward-scattered muons from downward-going cosmic rays are an irreducible background near the horizon. We estimate the background rate in each Up-μ sample at the horizon by comparing the cosmic muon rate at the azimuthal directions with the smallest and largest mountain overburdens. The estimated number of background events is subtracted from the number of Up-μ data events with reconstructed directions near the horizon. The estimated Up-μ background rate is approximately 2%-3%.
The event classification outlined in this section corresponds to 19 distinct analysis samples. There are eight FC CC-enhanced single-ring samples which are used for data collected during the SK I-III phases, while the remaining 11 FC multi-ring, NC π 0 , PC, and Up-μ samples are used for data taken during all phases. The data and MC counts for the 19 standard analysis samples are presented in Fig. 5 as a function of reconstructed lepton momentum or zenith angle.
This analysis modifies the FC single-ring event selection during the SK IV and SK V phases based on the number of 0 100 200 Sub-GeV e-like 1 d.e.* 0.5 1 1.5 0 500 1000 Sub-GeV e-like 0 d.e.* 0.5 1 1.5 0 100 200 300 400 -like 0 d.e.* P Sub-GeV 0.5 1 1.5 0 200 400 600 800 -like 1 d.e.* P Sub-GeV 0.5 1 1.5 0 200 400 Stopping P Up-0.5 1 1.5 0 200 400 -like 0 S 1-ring 0.5 1 1.5 0 50 100 * e Q Multi-GeV e-like 0.5 1 1.5 0 100 200 300 * e Q Multi-GeV e-like 0.5 1 1.5 0 50 100 150 PC Stopping 0.5 1 1.5 0 500 1000 1500 Non-showering P Up-0.5 1 1.5 0 50 100 150 -like 2 d.e.* P Sub-GeV 0.5 1 1.5 0 200 400 e Q Multi-Ring e-like 0.5 1 1.5 0 100 200 300 e Q Multi-Ring e-like 0.5 1 1.5 0 500 1000 PC Through-going 0.5 1 1.5 0 100 200 300 Showering P Up-0.5 1 1.5 0 500 1000 -like 0 S 2-ring 0.5 1 1.5 0 100 200 300 -like* P Multi-GeV 0.5 1 1.5 0 200 400 600 -like P Multi-Ring 0.5 1 1.5 0 100 200 Multi-Ring Other 0.5 1 1.5 V Super-Kamiokande I Data CC e Q CC e Q CC P Q CC IV, with best-fit systematic pulls applied. Different colors in the MC histograms correspond to the true neutrino flavors present in each sample. Panels beneath each distribution show the data-MC ratio, and all error bars are statistical. FC (non-PC and non-Up-μ) sub-GeV and multi-GeV single-ring samples, marked with an asterisk (*), contain events from only SK I-III, while the corresponding events from SK IV-V are separated into samples using tagged neutron information, shown in Fig. 6. observed neutron captures on hydrogen. The modification is motivated by the greater average neutron production in antineutrino interactions relative to neutrino interactions: Additional neutrons are expected in antineutrino events from proton-to-neutron (p
where l denotes lepton flavor, and also in both CC and NC deep inelastic scattering (DIS) processes due to the larger fraction of energy transferred to the recoiling hadronic system. We have incorporated neutron information into the analysis sample definitions to create ten additional samples, five for sub-GeV events and five for multi-GeV events, as follows:
TABLE III. Monte Carlo CC and NC purities by sample and data event counts used in this analysis. Purities and MC counts are shown with oscillation probabilities applied and without the effects of systematic pulls. The "ν τ CC" column shows the purity of both ν τ and ντ CC events. The Up-μ data counts are shown after background subtraction. MC purity Events Sample Energy bins cos θ z bins ν e CC νe CC ν μ CC νμ CC ν τ CC NC MC Data Fully contained (FC), single ring, sub-GeV SK I-III: e-like 0 decay-e 5 e AE momentum 10 in ½-1; 1 0.733 0.226 0.003 0.001 0.000 0.036 6409.1 6647 1 decay-e 5 e AE momentum single bin 0.796 0.016 0.086 0.020 0.001 0.081 612.0 682 μ-like 0 decay-e 5 μ AE momentum 10 in ½-1; 1 0.027 0.008 0.704 0.149 0.001 0.112 2153.9 2419 1 decay-e 5 μ AE momentum 10 in ½-1; 1 0.001 0.000 0.677 0.291 0.000 0.030 4241.4 4476 2 decay-e 5 μ AE momentum single bin 0.001 0.000 0.948 0.029 0.001 0.022 330.7 336 SK IV-V: ν e -like 5 e AE momentum 10 in ½-1; 1 0.794 0.016 0.090 0.024 0.001 0.074 943.7 1093 νe -like 0n 5 e AE momentum 10 in ½-1; 1 0.789 0.175 0.003 0.001 0.000 0.031 5961.8 6669 νe -like 1n 5 e AE momentum 10 in ½-1; 1 0.582 0.367 0.003 0.002 0.001 0.044 2266.6 1668 ν μ -like 5 μ AE momentum 10 in ½-1; 1 0.011 0.003 0.755 0.173 0.000 0.057 6596.0 7879 νμ -like 5 μ AE momentum 10 in ½-1; 1 0.001 0.000 0.533 0.417 0.001 0.049 2150.4 1793 Fully contained (FC), single ring, multi-GeV SK I-III: ν e -like 4 e AE momentum 10 in ½-1; 1 0.568 0.086 0.102 0.014 0.039 0.190 359.6 383 νe -like 4 e AE momentum 10 in ½-1; 1 0.556 0.341 0.013 0.002 0.012 0.075 1359.8 1339 ν μ -like 2 μ AE momentum 10 in ½-1; 1 0.002 0.001 0.621 0.371 0.003 0.002 1588.5 1564 SK IV-V: ν e -like 4 e AE momentum 10 in ½-1; 1 0.607 0.087 0.098 0.015 0.033 0.159 584.1 643 νe -like 0n 4 e AE momentum 10 in ½-1; 1 0.637 0.287 0.008 0.002 0.007 0.058 866.4 986 νe -like 1n 4 e AE momentum 10 in ½-1; 1 0.435 0.460 0.009 0.002 0.015 0.079 736.1 616 ν μ -like 2 μ AE momentum 10 in ½-1; 1 0.002 0.001 0.695 0.297 0.003 0.001 1464.0 1619 νμ -like 2 μ AE momentum 10 in ½-1; 1 0.001 0.000 0.446 0.549 0.004 0.002 593.1 503 SK I-V common samples: Fully contained (FC) Sub-GeV NC π 0 -like Single-ring 5 e AE momentum single bin 0.219 0.064 0.018 0.002 0.001 0.696 748.0 868 Two-ring 5 π 0 momentum single bin 0.096 0.028 0.015 0.001 0.000 0.860 2095.3 2494 Fully contained (FC) Multi-GeV, multi-ring ν e -like 3 visible energy 10 in ½-1; 1 0.495 0.066 0.175 0.014 0.035 0.215 2149.7 2411 νe -like 3 visible energy 10 in ½-1; 1 0.519 0.260 0.052 0.007 0.025 0.138 1210.0 1131 μ-like 4 visible energy 10 in ½-1; 1 0.028 0.003 0.713 0.200 0.006 0.050 3257.7 3427 Other 4 visible energy 10 in ½-1; 1 0.203 0.023 0.257 0.014 0.074 0.429 837.9 982 Partially contained (PC) Stopping 2 visible energy 10 in ½-1; 1 0.089 0.034 0.559 0.262 0.011 0.045 641.6 689 Through-going 4 visible energy 10 in ½-1; 1 0.006 0.002 0.638 0.341 0.007 0.006 3310.2 3397 Upward-going muons ðUp -μÞ Stopping 3 visible energy 10 in ½-1; 0 0.008 0.003 0.646 0.340 0.000 0.003 1574.3 1753.8 Nonshowering single bin 10 in ½-1; 0 0.002 0.001 0.662 0.334 0.000 0.001 5315.8 6423.9 Showering single bin 10 in ½-1; 0 0.001 0.000 0.671 0.327 0.000 0.001 1051.4 1110.6
FC single-ring e-like: Both sub-GeVand multi-GeVevents are divided into three samples: Events with one or more decay electron and any number of neutrons are classified as ν e -like. Events with no decay electrons are considered νe -like, and they are further separated into two samples based on whether or not there is at least one tagged neutron.
FC single-ring μ-like: Both sub-GeV and multi-GeV events are divided into two subsamples: Events with exactly one decay electron and one or more tagged neutrons are considered νμ -like; otherwise, they are considered ν μ -like.
More details about the neutron-tagging algorithm and the event selection modifications may be found in Refs. [22,23]. Figure 6 shows the data and MC comparisons for the SK IV-V neutron tagged samples. Compared to the equivalent antineutrino samples shown in Fig. 5, the SK IV-V νe samples requiring at least one tagged neutron have a higher purity of true νe MC events, summarized in Table III.
Multi-ring events in SK often result from the presence of energetic final-state pions. Charged pions are typically reconstructed as μ-like rings, while neutral pions generate e-like rings from their decays. The mixture of ring flavors in multi-ring events can cause ambiguity when attempting to classify the event by the interacting neutrino's flavor. On the other hand, the additional rings provide extra observable information not available in single-ring events. To classify multi-ring events based on a synthesis of the available information, we consider seven reconstructed quantities [24]. As in the single-ring classification, the primary quantity used to separate μ-like events from e-like events is the PID score of the most energetic ring in the event. To further improve the e -μ separation, we also consider the maximum distance of any decay electrons from the event vertex, L Decaye . L Decaye differs between μlike and e-like events because, while both electron-and muon-neutrino CC interactions may have decay electrons produced from pion decays, the muons produced in muon neutrino CC interactions typically have higher momenta than additional charged pions, and so the muons travel further before decaying. To separate ν e -like events from νelike events, we again consider the number of decay electrons. We also utilize the additional kinematic information present in multi-ring events by comparing the number of rings and two additional variables based on the distribution of ring momenta, the fraction of the momentum carried by the most energetic ring, p Max =p Tot , and the fraction of momentum transverse to the most energetic ring's direction, p T =p Tot These kinematic variables differ between ν e and νe events due to the different angular and momentum dependencies between neutrino and antineutrino cross sections. Finally, we use the visible energy to separate ν τ events which only occur at higher energies due to the energy threshold of their interactions. Based on these variables, we create four analysis samples for multi-ring events: CC ν e -like, CC νe -like, CC μ-like (both ν μ and νμ ), and "other." The "other" category is intended to separate NC and ν τ interactions from the CC-enhanced categories.
Previous versions of this analysis used a likelihoodbased classification algorithm with the seven reconstruction quantities as inputs [24]. This analysis replaces the likelihood-based selection with a boosted decision tree (BDT) utilizing the same set of variables as in the likelihood-based classification [25]. The BDT is trained on labeled true multi-ring MC events equivalent to a 430-year exposure and validated on MC events equivalent to a 70-year exposure for each SK phase. We optimize the BDT for sensitivity to the neutrino mass ordering by weighting the training events. Adjusting the weights of CC ν e and νe training events relative to CC ν μ and "other" events changes the expected signal purities obtained from the BDT selection, which in turn affects the expected mass-ordering sensitivity. Compared to the likelihood-based classification, we find that the optimized BDT increases the efficiencies of correctly classifying true CC events from 35% to 55% for ν e , from 58% to 65% for νe , and from 74% to 81% for ν μ , while decreasing the efficiency for correctly classifying NC events from 57% to 31%. The purity of true ν e CC events in the ν e -like sample decreases from 56% to 50%, while the purities of correctly classified CC events in the νe -like and ν μ -like samples (26% and 91%, respectively) are unchanged, and the purity of correctly classified NC events in the "other" sample increases from 30% to 43%. Accounting for differences in exposure, the increased efficiencies in the CC samples result in proportionally larger statistics in the CC ν e -like, νe -like, and ν μ -like multiring samples and decreased statistics in the "other" sample compared to the previous analysis.
Figure 7 shows the distributions of the BDT input variables from MC events. The largest differences between electron-and muon-flavor neutrino interactions are in the PID of the most energetic ring and the maximum distance travelled by decay electrons, L Decaye . Electron neutrino and antineutrino multi-ring events differ in the expected number of decay electrons, the number of rings, and in the fractional momentum variables p Max =p Tot and p T =p Tot Figure 8 demonstrates the ν e and νe separation performance of the BDT in the SK IV-V phases. Neutron information from the SK IV-V phases was also considered as an input variable to the BDT, but it was found to only result in a marginal improvement to ν e -ν e separation while introducing a large amount of systematic uncertainty.
Previous SK analyses used a 22.5 kiloton fiducial volume for FC events defined by requiring that the reconstructed event vertex be at least 200 cm away from the ID walls (conventional fiducial volume). This analysis includes events reconstructed within the 4.7 kiloton region 100-200 cm from the ID walls (additional fiducial volume) for all SK phases, representing a 20% increase in exposure. To include events in the additional fiducial volume region, both the reconstruction algorithms and background estimations required reevaluation [26]. In 2018, the pretabulated expected charge distributions, which are the basis of the ring counting and PID likelihoods used in APFIT, were recomputed as a function of a particle's distance to the nearest detector wall. These new charge tables encode the reduced number of hits and increased angular dependence expected from events near the detector walls. Adopting the updated charge tables reduced the e-μ misidentification of single-ring events in the additional region by up to 35%.
Cosmic muons with vertices misreconstructed inside the fiducial volume also become more prevalent closer to the detector walls. To estimate the size of these backgrounds, events in the additional fiducial volume region were eyescanned by experts. A slightly higher background rate of 0.5%, compared to 0.1% in the conventional fiducial volume, was found. The increased background rate in the additional fiducial volume is acceptable, although the background rate rises quickly for events with reconstructed vertices less than 100 cm from the detector walls, prohibiting further expansion.
We note that Ref. [26] included events from the additional fiducial volume region in the context of a proton decay search which was primarily concerned with sub-GeV atmospheric neutrino backgrounds. In contrast, this oscillation analysis includes both sub-GeV and multi-GeV atmospheric neutrino events. Additional studies to those presented in Ref. [26] were needed to confirm the performance of the reconstruction algorithms for multi-GeV events in the additional fiducial volume region. We studied the reconstructed momentum and direction biases and
1 2 3 4 5 6 Number of rings 0 0.1 0.2 0.3 0.4 0.5 0.6 Events (arb. units) 2000 1000 0 1000 2000 Most energetic ring PID 0 0.02 0.04 0.06 0.08 0.1 0.12 0 1 2 3 4 5 6 7 8 9 10 Number of decay electrons 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0 500 1000 1500 (cm) decay e.
L 5 10 4 10 3 10 2 10 1 10 1 0 5000 10000 (MeV) Vis. E 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Events (arb. units) 0 0.2 0.4 0.6 0.8 1 Tot. /p T p 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0 0 . 5 1 Tot. /p Max. p 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 Super-Kamiokande IV-V MC CC e Q CC e Q CC P Q & 7. Monte Carlo distributions of input variables for the FC multi-GeV, multi-ring BDT. Distributions are area-normalized. L Decay e is the maximum distance of any decay electron from the event's primary vertex. p T is the momentum transverse to the most energetic ring, p Max . is the momentum of the most energetic ring, and p Tot is the sum of the momenta from all reconstructed rings. 8. Data and MC distributions of the multi-ring BDT classifier scores for multi-ring ν e and νe events. The BDT makes a final classification decision for each event based on the class with the highest score: Negative values correspond to νe -like events, and positive values correspond to ν e -like events. Distributions are shown for the SK IV-V phases, including events in the expanded fiducial volume region. MC events are shown with oscillations applied, and without the effect of any fitted systematic uncertainty parameters. The "MC total" distribution includes contributions from true ν μ , ν τ , and NC events. resolutions using MC events, and we confirmed that these quantities were similar between events in the conventional and additional fiducial volume regions. We also studied the data versus MC agreement between the ring PID and ringcounting likelihoods for events in the additional fiducial volume region, as shown in Fig. 9. The figure compares the data and MC distributions for sub-GeV and multi-GeV events in both fiducial volume regions. The level of agreement in the likelihood distributions using the updated reconstruction algorithms is found to be equivalent between events in the conventional and expanded fiducial volume regions at both sub-GeV and multi-GeV energies.
This analysis is the first SK atmospheric neutrino oscillation analysis to include FC events from the additional fiducial volume region, resulting in a total fiducial volume of 27.2 kilotons for all SK phases. The total exposure for FC events in this analysis is 484.2 kilotonyears, a 48% increase over the previous published analysis.
We produce simulated Monte Carlo (MC) events to provide a prediction of atmospheric neutrino data. The SK MC consists of simulated neutrino interactions generated according to the flux model of Honda et al. [27] and the cross-section models of the NEUT simulation software [28,29]. NEUT also propagates intermediate particles within the nuclear medium and produces the final-state particles which exit the nucleus. These final-state particles are stepped through a GEANT3-based [30] simulation of the SK detector which models particle propagation in the detector water, Cherenkov radiation emission, and the detection of Cherenkov radiation by PMTs [16]. The SK detector simulation also implements a data-driven model of PMT electronics, dark noise, and disabled PMTs.
Because atmospheric neutrinos span several orders of magnitude in energy, multiple interaction processes are relevant across the different SK neutrino samples. Fully contained events, which are the lowest-energy atmospheric neutrinos in this analysis, have leading contributions from charged-current quasielastic (CCQE) and single-pion production processes. Partially contained and Up-μ events consist of higher-energy neutrinos where deep inelastic scattering (DIS) becomes dominant. The contribution to the event rate from each cross section process is shown in Fig. 10.
This analysis uses NEUT version 5.4.0 [29] for the nominal neutrino interaction models, which is an update from version 5.3.6 used in the previous analysis. Notable changes in this version include the replacement of the nominal one-particle-one-hole (1p1h) model, responsible for the CCQE and NCQE processes, from the Smith and Moniz relativistic Fermi gas (RFG) model [32] to the local Fermi gas (LFG) model by Nieves et al. [33]. The nominal value of the axial mass parameter M QE A was decreased from In these panels, the "CC QE" histograms include both 1p1h and 2p2h processes, and "CC Other" refers to single, nonpion hadron production. All MC distributions are calculated assuming neutrino oscillations and without the effect of any systematic uncertainties.
1.21 GeV=c to 1.05 GeV=c 2 . This change results in an overall decrease of CCQE events in the nominal MC by ∼20%, although the previous M QE A value is still within the 1σ uncertainty; see Sec. IV B. The LFG model also includes a correction to the cross section at low fourmomentum (q 2 ) transfers due to weak charge screening calculated via the random phase approximation (RPA) technique [34]. The BBBA05 [35] vector form factors are used in the nominal 1p1h model, which is unchanged from the previous analysis.
NEUT implements the Rein-Sehgal resonant single-pion production model [36], which is unchanged from the previous analysis. The coherent pion production crosssection calculation, which was previously based on the Rein-Sehgal model, has been updated to use the Berger-Sehgal model [37] for neutrino events with E ν < 10 GeV. This was done to improve the agreement in total cross section and angular distribution of the outgoing pion with recent scattering experiments [38].
Deep inelastic scattering cross sections in NEUT are calculated using the GRV98 parton distribution functions (PDFs) [39], with corrections to the low-q 2 regime from Bodek and Yang [40]. NEUT simulates the production of multiple hadrons using two separate models, selected based on the invariant mass of the hadronic system, W. For W < 2 GeV, a custom multi-pion generator is used, while for W > 2 GeV, NEUT uses PYTHIA v5.72 [41]. The DIS hadron production models are only used to produce final states with multiple hadrons to avoid overlap with the single-hadron models.
Final-state interactions (FSIs) are processes which modify particles exiting the nucleus due to intranuclear effects. Secondary interactions (SIs) are equivalent processes which occur in the detector medium instead of within the nucleus. NEUT implements four FSI þ SI processes: quasielastic scattering, charge exchange, pion absorption, and hadron production. Six parameters adjust the probabilities of these processes. NEUT 5.4.0 updated the default values of these parameters [42], which increased the pion absorption probability by 27%, and decreased the probability of charge exchange for pions with momentum < 400 MeV=c by 30% compared to version 5.3.6.
The cross section for the two-particle-two-hole (2p2h) process, in which no pions are produced but a pair of nucleons is ejected from the nucleus, is calculated using the Valencia model of Nieves et al. [43]. This is unchanged from the previous version of NEUT.
Atmospheric neutrino data are analyzed via the fit of binned MC counts to data. We bin atmospheric neutrino MC and data events using two-dimensional bins of reconstructed cosine zenith angle and momentum. The bin definitions for the 29 samples used in this analysis are listed in Table III. The zenith angle bin definitions were updated for this analysis and are discussed in Sec. IVA.
Systematic uncertainties are incorporated in the analysis as additional nuisance parameters which modify the nominal MC prediction in the fit. We estimate the effect of each systematic uncertainty by quantifying the change induced -μ). The simulated event rates are for an oscillated MC, with oscillation parameters set to the Particle Data Group best-fit values [31], except for δ CP , which is set to 0. "CC Other" modes refer to the production of single, nonpion hadrons.
in each analysis bin due to þ1σ and -1σ deviations from its nominal value. The systematic uncertainties considered in this analysis are described in Sec. IV B.
In the previous analysis, FC and PC events were binned into ten evenly spaced cosine zenith angle bins on the interval ½-1; 1. This separated the expected oscillation resonance region across two bins, reducing the signal-tobackground ratio in samples sensitive to the neutrino mass ordering. In this analysis, the zenith angle bins for FC and PC events have been updated to more precisely cover the resonance region, demonstrated in Fig. 11. The updated bins are defined by the edges -1, -0.839, -0.644, -0.448, -0.224, 0, 0.224, 0.448, 0.644, 0.839, and 1. With the present statistics, adopting the updated bins results in a negligible effect on the expected sensitivity of this analysis to the mass ordering. However, the updated bins reduce the ambiguity as to whether or not events in the signal region occur in the expected zenith angle range or slightly outside of it.
The analysis includes 193 independent systematic uncertainty sources: 48 describe atmospheric neutrino flux and cross section effects common to all phases, and the remaining 29 × 5 describe reconstruction efficiencies which are separately estimated for each data-taking phase. Additional details on the formulation of the systematic uncertainties may be found in Ref. [44].
The atmospheric neutrino flux uncertainties in this analysis are unchanged from the previous analysis [5]. Two energy-dependent uncertainties scale the normalization of the flux above and below 1 GeV. Additionally, there are three energy-dependent uncertainties which modify the ratios of muon-to-electron, electron-to-antielectron, and muon-to-antimuon flavor neutrinos present in the atmospheric neutrino flux. An uncertainty based on the ratio of the Honda flux calculation with a modified kaon-to-pion ratio over the nominal model is also included [45]. There are two additional uncertainties which allow modifications to the flux shape. The first modifies the ratio between upwardgoing (cos θ z < 0) and downward-going (cos θ z > 0) events, while the second modifies the ratio between vertical (cos θ z > 0.9) events and horizontal (0 < cos θ z < 0.1) events. The effect size of the shape uncertainties is separately estimated for each neutrino flavor and for sub-GeVand multi-GeV events.
The effects of uncertainties on parameters in the CCQE and single-pion cross-section models are computed by reweighting MC events by the ratio of the doubledifferential cross section after þ1σ and -1σ changes in each parameter. For CCQE events, M QE A is taken to be ð1.05 AE 0.16Þ GeV=c 2 . There are three parameters which control the single-pion production cross section in the Rein-Sehgal model: the axial mass, M Res A ¼ ð0.95AE 0.15Þ GeV=c 2 ; the axial form factor coefficient, C 5 A ¼ 1.01 AE 0.12; and the isospin-1 2 background contribution scaling parameter, I1
2 ¼ 1.30 AE 0.20.
A uncertainty, we place additional uncertainty on CCQE events due to differences between the LFG and RFG models. There are five uncertainties where the 1σ effect is computed by forming ratios between the predictions of the two models: the absolute event rate (both sub-GeV and multi-GeV), the shape of the E ν dependence, and the ratios of ν μ =ν e and νμ =ν e . These uncertainties are computed as a function of neutrino energy. The uncertainty on the 2p2h contribution to CCQE interactions is set to 100% due to a lack of direct measurements, and it is unchanged from the previous analysis.
An uncertainty on FSIs and SIs is implemented by recomputing the final states for MC events using multiple sets of the six NEUT FSI parameters. Each set consists of variations of one or more FSI parameters by their AE1σ uncertainties as determined by a fit to external pion scattering from Ref. [46]. The set which produces the largest change in classification outcome is taken as a conservative estimate of the FSI þ SI uncertainty.
Two neutron-related systematic uncertainties are included in this analysis to account for the dependence on neutron production models in the neutron-based event selection used for SK IV and SK V data. Variations in these systematic uncertainties move events between samples with no tagged neutrons and samples with one or more tagged neutrons. The largest of the two uncertainties originates from measurements of neutron production as a function of transverse momentum by T2K [47]. An additional neutronrelated systematic uncertainty accounts for differences in neutron multiplicity predicted by the different neutrino interaction generators NEUT and GENIE [48].
In addition to model parameter uncertainties, we assign a 20% uncertainty on the ratio of NC to CC events and a 25% uncertainty on the tau neutrino CC cross section [49]. The tau neutrino cross section uncertainty was estimated from a comparison between NEUT and the cross section model of Hagiwara et al. [50].
Ring reconstruction systematic uncertainties include efficiencies of ring-based particle identification (PID), ring counting, and identification of π 0 decays. The estimation of these systematic uncertainties is unchanged from the previous analysis and is performed using a "scale-andshift" procedure: First, MC events are labeled as either signal or background, where the particular label differs for each systematic uncertainty source. For example, the ringcounting systematic uncertainty estimation considers events with a single charged particle with true momentum sufficient to produce a Cherenkov ring as signal, and events with multiple particles capable of producing rings as background. Similarly, the PID uncertainty estimation considers true e-like events as signal and true μ-like events as background. Once labeled, the signal and background distributions of each reconstruction quantity x (e.g., the ring-counting or PID likelihoods) are scaled and shifted by a linear function, x 0 ¼ β 0 þ β 1 x, and fit to a subset of atmospheric neutrino data. The signal and background distributions have independent β parameters. The fitted β parameters are then fluctuated according to their fitted uncertainties to generate toy datasets with random amounts of signal and background events. The maximum fractional change in signal purity observed in the toy datasets is taken as a conservative estimate of the uncertainty.
The scale-and-shift procedure is used to separately evaluate ring-related systematic uncertainties for events in the conventional and additional fiducial volume regions during each data-taking phase. The ring-counting and PID systematic uncertainties are separately computed for each combination of e-like and μ-like, and sub-GeV and multi-GeV events. As in the previous analysis, the scale-and-shift procedure is used to place an additional uncertainty on the normalization of ν μ contamination present in the e-like samples, estimated from the fitted fraction of ν μ events classified as e-like. We also use the scale-and-shift procedure to estimate the uncertainty associated with the multi-ring BDT: Two sources of systematic uncertainty are considered for multi-ring events based on the data-MC agreement in the input distributions used to train the BDT, and in the output BDT scores themselves. The changes in the event rate in each multi-ring sample obtained through fluctuating the input and output distributions according to their fitted uncertainty are taken as the 1σ effects.
As in the previous analysis, this analysis takes the maximum data-MC disagreement between the decay electron, π 0 mass, and stopping muon calibration sources during each SK phase, shown in Fig. 3, as an estimate of the overall energy scale uncertainty. Variations in this uncertainty increase or decrease the reconstructed momenta of all MC events, causing these events to move between momentum bins. The energy scale uncertainty was newly estimated for the SK V phase and is 1.8% in the conventional fiducial volume and 2.0% in the additional fiducial volume.
This analysis introduces a systematic uncertainty for the neutron-tagging efficiency and false detection rate during the SK IV and SK V phases, which is factored into two components: detector conditions and neutron travel distance. The contribution to the uncertainty due to detector conditions was estimated from comparisons between data from an AmBe source [51] and MC events generated with variations in PMT gain, PMT dark rate, and water transparency within the ranges observed during the SK IV-V phases. The neutron travel distance contribution to the uncertainty was estimated from an atmospheric neutrino data-MC comparison of the distance between the primary neutrino event vertex and the reconstructed neutron capture vertex, convolved with the estimated capture efficiency as a function of the distance. The quadrature sum of these effects sets the overall neutron-tagging uncertainty. Similarly to the neutron production uncertainties, variations in the tagging efficiency uncertainty move events between samples requiring no tagged neutrons and samples requiring one or more tagged neutrons.
The systematic uncertainties with the largest effect on the analysis's sensitivity to the neutrino mass ordering are those which can change the ratio of upward-going (cos θ z < 0) versus downward-going (cos θ z > 0) events in the multi-GeV e-like samples. This is because the mass-ordering signal region occurs in the upward-going bins of these samples, and downward-going neutrino events act as an unoscillated dataset which constrains uncertainties without zenith angle dependence independent of oscillation effects. Figure 12 shows the estimated largest sources of uncertainty which affect the ratio of upward-versus-downward going events in the multi-GeV e-like samples. The largest source of uncertainty is statistics, indicating the mass-ordering sensitivity is statistics-limited. The ν τ cross-section uncertainty is the next largest: Tau neutrino interactions at SK are almost exclusively upward-going, since they are the result of ν μ → ν τ oscillations, so there are no downward-going tau neutrino events to constrain the ν τ cross-section uncertainty. Additionally, tau neutrino interactions tend to be reconstructed as multi-GeV, multi-ring e-like events, which places them in the mass-ordering signal region [52]. In contrast, uncertainties on the overall flux and cross-section normalizations are less important for the mass-ordering analysis, since these modify the predicted number of events independent of zenith angle.
This work presents two fits of atmospheric neutrino data: First, we perform an "atmospheric-only" fit which measures Δm 2 32;31 , sin 2 θ 23 , δ CP , sin 2 θ 13 , and the neutrino mass ordering using the SK atmospheric neutrino data with no external constraints. Next, we analyze the same data including a constraint on sin 2 θ 13 from reactor neutrino experiments.
In both fits, we evaluate data-MC agreement by computing a χ 2 statistic at each point on a fixed grid of neutrino oscillation parameters, defined for each fit in Table IV. Oscillation probabilities are applied to the MC events using the oscillation parameters at each point in the grid, and systematic uncertainty nuisance parameters are fit to determine the lowest χ 2 value. The best-fit neutrino oscillation parameters are taken to be the grid point with the lowest χ 2 value.
The χ 2 statistic is computed via a summation over n bins assuming Poisson statistics, and with systematic uncertainty pull terms, ϵ i , which have units of σ. Each ϵ i scales the effect of the ith systematic uncertainty and is included in the χ 2 calculation as an additional penalty term:
In Eq. (10), n indexes each bin, and O n are the observed counts in the nth bin. The expected counts E n are calculated from a nominal MC prediction E 0 n scaled by the effect of systematic uncertainties. The effect of the ith systematic uncertainty is calculated as ϵ i times the fractional change induced in the nth bin by a 1σ change, f i;n :
The constraint ∂χ 2 =∂ϵ i ¼ 0 at the minimum χ 2 value yields a system of equations which may be solved to find the configuration of ϵ i 's that produce the smallest χ 2 for a given MC prediction.
We studied the impact of using neutron captures on hydrogen with a tagging efficiency of 26% on the FIG. 12. Estimated fractional uncertainty on the ratio of upward-going, cos θ z < 0, versus downward-going, cos θ z > 0, events in the multi-GeV e-like samples (both single-ring and multi-ring) from statistical and systematic uncertainties. The fractional uncertainties are estimated by randomly fluctuating the nominal MC counts either by Poisson statistics for the "Statistics" entry, or by Gaussian fluctuations of one or more systematic uncertainties. Uncertainties on the up/down ratio in the multi-GeV e-like samples have the largest impact on the sensitivity to the neutrino mass ordering.
TABLE IV. The oscillation parameter grid definitions used in the three analyses presented in this work. The grids are evenly spaced and include the minimum and maximum values. The grids are the same for both the normal and inverted mass orderings in all analyses. Both analyses treat Δm 2 21 and sin 2 θ 12 as constrained parameters. Parameter Minimum Maximum Points SK, Atmospheric only sin 2 θ 23 0.3 0.775 20 Δm 2 32 (10 -3 eV 2 ) 1.2 3.6 25 δ CP (rad) -π π 21 sin 2 θ 13 0 0.075 16 SK, sin 2 θ 13 Constrained sin 2 θ 23 0.3 0.7 35 Δm 2 32 (10 -3 eV 2 ) 1.0 4.9 40 δ CP (rad) -π π 37 sin 2 θ 13 0.0220 AE 0.0007 1 S /2 S 0 / 2 S S CP G 0 5 10 2 F ' SK I-V expanded FV sensitivity With neutron tagging Inverted No neutron tagging Normal FIG. 13. 1D Δχ 2 sensitivity profile of δ CP for the normal and inverted mass-ordering scenarios, with (solid lines) and without (dashed lines) neutron-based event classification for SK IV and SK V data. The sensitivity is computed assuming the true neutrino oscillation parameters are the global best-fit parameters from [31] and the normal mass ordering. The parameters Δm 2 32 , sin 2 θ 23 , and δ CP are free parameters in the fit, while the other oscillation parameters are fixed.
sensitivity to the mass ordering. As discussed in Sec. II B 2, tagged neutrons provide an additional handle during event reconstruction which improves the statistical separation of neutrinos and antineutrinos in SK versus only considering the number of decay electrons. This corresponds to an increased purity of correct-sign ν e and νe events in samples where sensitivity to the mass ordering is expected. The increased purity can be seen by comparing the FC multi-GeV single-ring e-like samples listed in Table III: The SK I-III multi-GeV νe -like sample requiring no decay electron has a MC purity of 34% true νe CC interactions, while the SK IV-V multi-GeV νe -like sample requiring no decay electrons and one tagged neutron has a MC purity of 46% true νe CC interactions.
Figure 13 shows a comparison of the sensitivities to the mass ordering and δ CP obtained using the standard 19 samples discussed in Sec. II B for all SK phases versus the present configuration, which separates out the FC singlering samples from the SK IV and SK V phases (approximately 57% of the total exposure) into additional samples which use neutron information. The figure demonstrates that using the neutron-tagged samples improves the expected sensitivity to the mass ordering and δ CP by producing higher Δχ 2 values [see Eq. ( 10)] away from the true oscillation parameters. Using the neutron-tagged samples has a negligible effect on the sensitivities to sin 2 θ 23 and Δm 2 32 , which are primarily constrained by the ν μ -enriched samples.
Results in this section and the next are presented as Δχ 2 contours of one or two oscillation parameters taken with respect to the global best-fit point, with and without constraints on θ 13 . The contours are profiled: We draw the minimum Δχ 2 values among all other combinations of oscillation parameters scanned with one or two oscillation parameters fixed.
Figure 14 shows the one-dimensional (1D) Δχ 2 contours in the atmospheric-only analyses for the fitted neutrino oscillation parameters δ CP , sin 2 θ 13 , Δm
2 32;31 , and sin 2 θ 23 , with respect to the best-fit point across both S /2 S 0 / 2 S S CP G 0 2 4 6 8 10 12 14 16 2 F ' V expanded FV SK I Data fit Inverted MC expectation Normal 0 0.02 0.04 0.06 13 T 2 sin 0 2 4 6 8 10 12 14 16 2 F ' 1.5 2 2.5 3 3.5 ) 2 eV -3 (10 32,31 2 m ' 0 2 4 6 8 10 12 14 16 2 F ' 68% 90% 95% 99% 0.4 0.5 0.6 23 T 2 sin 0 2 4 6 8 10 12 14 16 2 F ' FIG. 14. 1D Δχ 2 profiles of oscillation parameters in the analysis, with sin 2 θ 13 treated as a free parameter. Solid lines correspond to the data fit result, while dashed lines correspond to the MC expectation at the data best-fit oscillation parameters, cf. Table IV. Dotted lines show critical values of the χ 2 distribution for 1 degree of freedom, with the lowest to highest corresponding to 68%, 90%, 95%, and 99% probabilities. mass orderings. The normal ordering is preferred: The difference between the minimum χ 2 in the inverted and normal orderings is Δχ 2 I:O:-N:O: ¼ 5.23. The MC expectation for the mass ordering at the best-fit oscillation parameters is Δχ 2 I:O:-N:O: ¼ 1.53, indicated by dashed lines in Fig. 14. The difference between the data result and the MC expectation is discussed in Sec. V.
Sensitivity to δ CP in the SK data originates from both sub-GeV and multi-GeV e-like samples, where values of δ CP near -π=2 indicate increased ν e appearance and decreased νe appearance relative to δ CP ¼ 0. The best-fit value of δ CP is -1.89 in both orderings, signaling increased ν e appearance. The constraints on δ CP and sin 2 θ 13 assuming the inverted ordering are weaker, although consistent with the best-fit values assuming the normal ordering. These weaker constraints are an expected consequence of the freedom to adjust sin 2 θ 13 and δ CP simultaneously, combined with the reduced antineutrino statistics relative to neutrino statistics in the atmospheric neutrino data.
The preferred values of sin 2 θ 13 are 0.020 and 0.010 in the normal ordering and inverted ordering, respectively. Both results are consistent with the reactor-preferred value of sin 2 θ 13 ¼ 0.0220 at the 90% confidence level, and sin 2 θ 13 ¼ 0 is disfavored in the normal ordering at the 99% level. The data favor the magnitude of the fitted squared-mass difference jΔm 2 32;31 j ¼ 2.40 × 10 -3 eV 2 in both orderings. The nonsmooth behavior of the constraints on this parameter, especially evident in the inverted ordering fit, is a consequence of rapidly varying oscillation probabilities in the sub-GeV samples. Finally, the atmospheric neutrino data place the best-fit value of sin 2 θ 23 in the lower octant, sin 2 θ 23 ¼ 0.45, although values in each octant are allowed at the 68% level.
As discussed in Sec. I A, the combination of nonzero sin 2 θ 13 and a normal neutrino mass ordering leads to electron neutrino appearance for upward-going multi-GeV events. We observe excess electron-flavor upward-going multi-GeV single-ring and multi-ring events in the SK data. Up-down asymmetry for multi-GeV e-like events. The y axis is the asymmetry parameter, the ratio between the difference and sum of upward-going (cos θ z < 0.6) and downward-going (cos θ z > 0.6) events. The x axis is the reconstructed neutrino energy: For single-ring events, the reconstructed energy is the visible energy of the ring assuming the reconstructed ring is an electron, while for mult-iring events, it is the total visible energy of the event. All error bars are statistical. MC lines for the normal and inverted orderings are drawn assuming the best-fit oscillation parameters of the analysis with sin 2 θ 13 constrained. SK IV-V multi-GeV single-ring events are selected using the number of tagged neutrons, and so they are separated from the SK I-III multi-GeV single-ring samples.
Figure 15 shows a projection of the multi-GeV e-like samples as an up-down asymmetry
where "Up" is the number of upward-going (cos θ z < -0.6) events, and "Down" is the number of downward-going (cos θ z > 0.6) events in each sample.
The figure plots the asymmetry for these data as a function of reconstructed energy and the expected asymmetry for the normal and inverted ordering scenarios, assuming the bestfit oscillation parameters from the fit to all atmospheric neutrino data. The ν e -enhanced samples, multi-GeV ν e -like and multi-ring ν e -like, have the largest excesses relative to either ordering, and drive the preference for the normal mass ordering in the analysis.
2. Results with reactor constraints on sin 2 θ 13
Figure 16 shows the 1D Δχ 2 profiles for the fitted neutrino oscillation parameters, assuming the constraint sin 2 θ 13 ¼ 0.0220 AE 0.0007 from reactor antineutrino disappearance experiments [31]. The constraint on sin 2 θ 13 is incorporated by introducing an additional systematic uncertainty for this fit, where the 1σ effect is defined as the change induced by varying sin 2 θ 13 by its measured 1σ uncertainty.
The best-fit value of δ CP in both the normal and inverted orderings for the fit with sin 2 θ 13 constrained is -1.75, which is consistent with the atmospheric-only analysis at the 1σ level. This fit also finds improved constraints on δ CP in the inverted ordering for values near π=2: The constraint on sin 2 θ 13 fixes the effect size of the mass ordering, such that the separate modifications to ν e appearance from δ CP are more readily resolved.
In this fit, the preference for the normal ordering increases to Δχ 2 I:O:-N:O: ¼ 5.69. This improvement is consistent with the observed preference for smaller values of sin 2 θ 13 in the inverted ordering fit with sin 2 θ 13 free: The χ 2 value in the inverted ordering increases with the added constraint, while the χ 2 value in the normal ordering remains similar to the result without the constraint.
Figure 17 shows the constraints on sin 2 θ 23 and Δm 2 32 from the θ 13 -constrained analysis of SK atmospheric neutrino data compared with the constraints from MINOS/MINOS+ [53], NOvA [4], T2K [3], and IceCube [54]. The SK atmospheric neutrino data are consistent with the other experiments at the 90% level. While the atmospheric neutrino data find a best-fit value of sin 2 θ 23 in the lower octant, we note that the previous publication found a best-fit value in the upper octant [5], and IV. Dotted lines show critical values of the χ 2 distribution for 1 degree of freedom corresponding to 68%, 90%, 95%, and 99% probabilities.
0.35 0.4 0.45 0.5 0.55 0.6 0.65 23 T 2 sin 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 ) 2 eV -3 (10 32 2 m ' Normal ordering, 90% C.L. MINOS/MINOS+ 2020 NOvA 2020 T2K 2023 IceCube 2023 Super-K
FIG. 17. 2D constant Δχ 2 contours of neutrino oscillation parameters Δm 2 32 and sin 2 θ 23 for the normal mass ordering. Contours are drawn for a 90% critical χ 2 value assuming 2 degrees of freedom, with the Δχ 2 computed for each experiment with respect to the best-fit point in the normal mass ordering. The Super-K contour shows the result of this analysis, and other contours are adapted from publications by MINOS+ [53], NOvA [4], T2K [3], and IceCube [54]. Best-fit points are indicated with markers for each experiment. that values of sin 2 θ 23 in both octants are allowed at the 1σ level.
Table V summarizes the fit results of the analyses presented in this work. In both analyses, the normal ordering is preferred, and the best-fit oscillation parameters predict weaker sensitivities to the neutrino mass ordering than the observed Δχ 2 I:O:-N:O: . To quantify the significance of the mass-ordering preference from the fit results, we generated ensembles of toy datasets to produce the distribution of the Δχ 2 I:O:-N:O: statistic. 3 Each toy dataset consists of fluctuated counts according to each bin's statistical uncertainty, which are scaled by randomly sampling the systematic uncertainty coefficients. Ensembles were generated assuming both the normal and inverted mass orderings and with oscillation parameters fixed at the best-fit points. We fit each toy dataset in each ordering with Δm 2 32 , sin 2 θ 23 , and δ CP as free parameters to compute Δχ 2 I:O:-N:O: . Figure 18 shows the distribution of Δχ 2 I:O:-N:O: compared with the data fit result for the atmospheric analysis with sin 2 θ 13 constrained. The probability of observing a more extreme result than the data (the p-value) is given by the area to the right of the data line in the normal-ordering scenario, and by the area to the left of the data line in the inverted-ordering scenario. While the figure shows the pvalue determined from simulated datasets for the inverted mass ordering is 0.0091, we note that with the present SK statistics, the expected sensitivity remains weak for rejecting either ordering. Indeed, the p-value for the data result within the normal ordering, 0.88, is not especially likely either. For the situation in which the data must select between two mutually exclusive hypotheses, the CL s method [56] provides an estimate of the p-value that considers simultaneous agreement from both hypotheses:
where p N:O: and p I:O: refer to the p-values in the normal or inverted ordering. This prescription decreases the significance for rejecting the inverted hypothesis proportional to the simultaneous significance of accepting the normalordering hypothesis. The CL s for the atmospheric analysis with sin 2 θ 13 constrained is 0.077, corresponding to a rejection of the inverted mass ordering at the 92.3% confidence level. This number is similar to the previous SK result, CL s ¼ 0.070 [5], despite originating from a larger Δχ 2 , 5.69 versus 4.33. While the mass-ordering sensitivity and data result both increased for this analysis, the probability of obtaining the data result simultaneously decreased in both orderings, resulting in a similar CL s value. The p-value obtained from toy datasets depends on the choice of oscillation parameters. While the atmospheric analysis places sin 2 θ 23 in the lower octant, values of sin 2 θ 23 spanning both octants are allowed at the 1σ level.
TABLE V. Best-fit neutrino oscillation parameters from the analyses presented in this work. The uncertainties on each oscillation parameter are the AE1σ allowed regions assuming a χ 2 distribution with 1 degree of freedom. The second-to-last column shows the total χ 2 . Both analyses have 930 bins. Fit result Ordering jΔm 2 32;31 j (10 -3 eV 2 ) sin 2 θ 23 sin 2 θ 13 δ CP ð-π; πÞ χ 2 Δχ 2 I:O:-N:O: SK, Atmospheric only Normal 2.40 þ0.07 -0.09 0.45 þ0.06 -0.03 0.020 þ0.016 -0.011 -1.89 þ0.87 -1.18 1022.06 5.23 Inverted 2.40 þ0.05 -0.33 0.48 þ0.07 -0.05 0.010 þ0.021 -0.008 -1.89 þ1.32 -1.97 1027.29 SK, sin 2 θ 13 Constrained Normal 2.40 þ0.07 -0.09 0.45 þ0.06 -0.03 Á Á Á -1.75 þ0.76 -1.25 1022.06 5.69 Inverted 2.40 þ0.06 -0.12 0.45 þ0.08 -0.03 Á Á Á -1.75 þ0.89 -1.22
Recent results from T2K favor maximal CP violation, δ CP ≈ -π=2[3], while recent measurements from NOvA disfavor CPviolating scenarios[4].
A single PMT imploded during refilling of the detector, creating a shockwave and chain reaction, destroying additional PMTs. Acrylic covers were mounted to all the PMTs to prevent similar accidents in the future.
Because the two mass-ordering scenarios are not nested hypotheses, taking the square root of Δχ 2 I:O:-N:O: to estimate the significance by invoking Wilks' theorem is not recommended[55].