Atmospheric neutrino oscillation analysis with improved event reconstruction in Super-Kamiokande IV

for the normal hierarchy, disfavoring the inverted hierarchy at 74% assuming oscillations at the best ﬁt of the analysis.


Introduction
Neutrino oscillations have been confirmed by a variety of experiments using both natural and artificial sources. At present the data are well described assuming mixing among all three active neutrinos with the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) formalism [1,2]. Though most of its parameters have been experimentally measured [3], the ordering of the mass states with the largest splitting (known as the mass hierarchy), the octant of the atmospheric mixing angle θ 23 , and the value of its CP-violating phase are unknown. These unresolved issues have been at the forefront of results from the T2K [4], NOvA [5], and Super-Kamiokande [6] (Super-K, SK) experiments and are the focus of next-generation experiments planned in the USA [7], China [8], and Japan [9].
The atmospheric neutrino flux provides a source of both neutrinos and antineutrinos with a wide variety of energies and path lengths suitable to probe each of these open questions. Matter effects influence the oscillations of neutrinos passing through the Earth and produce a resonance enhancement of the oscillation probability that depends on the values of θ 13 , θ 23 , and δ CP . Importantly, this resonance occurs only for neutrinos if the hierarchy is normal (largest splitting is between the two heaviest states) and only for antineutrinos if it is inverted. The hierarchy signal manifests most strongly as an upward-going excess of electron neutrino or antineutrino events with multi-GeV energies, a region where interactions frequently have multi-particle final states that complicate particle identification and where the flux is orders of magnitude lower than its peak. Accordingly, recent Super-Kamiokande results [6] have been limited by both a lack of statistics and event miscategorization in the signal region.
A new event reconstruction algorithm has been developed based on a maximum likelihood method designed to extract specific event topologies and determine the best set of kinematic parameters. The new algorithm demonstrates improved reconstruction performance across a variety of metrics including vertex resolution, particle identification probability, and momentum resolution and does so across a larger volume of the Super-K detector than its predecessor. As a result, it is expected to improve Super-K's sensitivity to the remaining parameters of the PMNS mixing paradigm.
This paper presents an analysis of 3118.5 days of SK-IV atmospheric data for a 253.9 kton·year exposure with a 32% larger fiducial volume (FV) than previous Super-K analyses. Section 2 briefly reviews atmospheric neutrino oscillations before the detector is introduced in Sect. 3. In Sect. 4 the atmospheric neutrino data reduction and event categorization are discussed. The new reconstruction algorithm and its performance in comparison to the conventional reconstruction are detailed in Sect. 5. Based on this performance improvement the expansion of the FV and corresponding systematic errors are presented in Sects. 6 and 7, respectively. An analysis of the atmospheric neutrino data by themselves is presented in Sect. 8.1 and is followed by an analysis in which θ 13 is constrained by the world data in Sect. 8.2, before concluding in Sect. 9.

Neutrino oscillations
Neutrino oscillations are a result of the neutrino mass eigenstates differing from their weakinteraction (flavor) eigenstates. For neutrinos in vacuum, their oscillation probability in the standard three-flavor paradigm can be described by six parameters, nearly all of which have been confirmed to PTEP 2019, 053F01 M. Jiang et al. have non-zero values based on the results of reactor, atmospheric, solar, and long-baseline neutrino experiments [3]. For brevity the details of the PMNS oscillation formalism will be omitted here; more details about the treatment of neutrino oscillation in Super-K can be found in Ref. [6].
In vacuum the leading terms in the oscillations of atmospheric ν e and ν μ can be written as: P(ν e → ν e ) ∼ = 1 − sin 2 2θ 13 sin 2 1.27 m 2 32 L E P(ν μ → ν μ ) ∼ = 1 − 4 cos 2 θ 13 sin 2 θ 23 (1 − cos 2 θ 13 sin 2 θ 23 ) × sin 2 1.27 m 2 32 L E P(ν μ ↔ ν e ) ∼ = sin 2 θ 23 sin 2 2θ 13 sin 2 1.27 m 2 32 L E (1) using the assumption that m 2 31 ≈ m 2 32 since m 2 21 is known to be small in comparison. In these equations the neutrino travel length is represented by L (km) and its energy by E (GeV). The unit of m 2 is eV 2 . When neutrinos propagate in matter, the oscillation probabilities are modified due to interactions with the electrons in the Earth. Indeed the forward scattering amplitude of ν e differs from those of the other flavors, producing an effective potential felt only by this species. This phenomenon is known as the matter effect or the MSW effect of neutrino oscillations [10,11]. When neutrinos travel through homogeneous matter, the oscillation parameters sin 2 θ 13 and m 2 32 in Eq. 1 are replaced by their matter-equivalents [12]: where A CC = ±2 √ 2G F N e E ν and the electron density N e is assumed to be constant. Here G F represents the Fermi constant, and the sign of A CC is positive for neutrinos and negative for antineutrinos.
Resonant enhancement of these matter variables occurs when A CC / m 2 32 = cos 2θ 13 and is dependent upon the sign of m 2 32 and whether the neutrino is a particle (A CC > 0) or an antiparticle (A CC < 0). Note that the enhancement only occurs for neutrinos if the hierarchy is normal ( m 2 32 > 0) whereas it only happens for antineutrinos if it is inverted ( m 2 32 < 0). Since atmospheric neutrinos generally experience a varying matter profile, and hence electron density N e changes as they travel through the Earth, they experience a variety of matter effects and similar enhancements to the oscillation probability. The calculation of oscillation probability in this analysis takes such variation on matter density into consideration, with a simplified version of the preliminary reference Earth model (PREM) (cf. Ref. [6]). Typically neutrinos with a few to ten GeV of energy that travel through the Earth's core experience matter effects the most strongly. From the equations above it can be seen that atmospheric neutrinos are sensitive to sin 2 2θ 13

The Super-Kamiokande detector
Super-Kamiokande is a cylindrical water Cherenkov detector with a total volume of 50 kilotons and is located inside the Kamioka mine in Gifu Prefecture, Japan. It is optically separated into two regions, an inner detector (ID), which forms the primary neutrino target and has a total volume of 32 kilotons, and an outer detector (OD), which is a two-meter cylindrical shell surrounding the ID and is used primarily as a veto. The ID is instrumented with more than 11 000 inward-facing 20-inch photomultiplier tubes (PMTs), representing 40% photocathode coverage of the target volume. A total of 1885 8-inch PMTs coupled to wavelength-shifting plates are mounted on the inner surface of the OD, while its outer surface is covered with reflective sheets to increase light collection.
Since the start of operations in 1996, Super-Kamiokande has gone through four data-taking periods, SK-I, -II, -III, and -IV. The present work focuses on the latter, which started in September 2008 and ended on 31 May 31 2018, when the detector began upgrade work for its next phase. Though the topologies of the detector in the SK-III and -IV periods are the same, at the start of SK-IV the frontend electronics were upgraded to a system based on ASIC that uses a high-speed charge-to-time converter [13]. After this upgrade Super-Kamiokande is able to collect all PMT hits above threshold without incurring any dead time. More detailed descriptions of the detector and its electronics are presented in Refs. [13][14][15].
Atmospheric neutrino interactions in the detector are simulated using the Honda et al. flux calculation [16] and the NEUT [17] neutrino interaction software (version 5.4.0). Particles emerging from the interactions are tracked through a simulation of the detector based on GEANT3 [18]. The present work uses an updated version of NEUT relative to the previous atmospheric neutrino analysis (cf. Ref. [6]). In particular, charged-current quasi-elastic (CCQE) interactions in the new simulation are on the local Fermi-gas model of Nieves [19,20], assuming an axial mass M A = 1.05 GeV/c 2 and using the random phase approximation correction. Further, deep inelastic scattering (DIS) is modeled using the GRV98 parton distribution function [21] and utilizing CKM matrix elements for the calculation of structure functions. Corrections for low q 2 scattering have been updated to those of Bodek and Yang [22], where q 2 denotes the square of the transferred four-momentum.

Event sample and reduction
The current analysis utilizes atmospheric neutrino data collected during the SK-IV period with a total livetime of 3118.45 days. As in previous Super-K analyses, the atmospheric neutrino data are separated into three general categories, fully contained (FC), partially contained (PC), and upwardgoing muons (Up-μ).
For events classified as FC and PC, the neutrino interacts within the fiducial volume, defined as the region located more than 2 m from the ID wall in previous analyses. Events with no activity in the outer detector are classified as FC. If energy deposition in the OD is observed, generally from a high-energy muon exiting the ID, the event is classified as PC. Muons created by neutrino interactions in the rock around SK or in the OD water and traveling upward through the detector form the Up-μ sample.
At trigger level the Super-Kamiokande event sample consists mainly of downward-going cosmic ray muons and low-energy radioactivity from contaminants in the water such as radon. To remove these backgrounds, reduction processes specific to the three main sample categories are applied to the data, details of which can be found in Ref. [23].
As an example, five steps of data reduction criteria are used for an FC sample. In the first and second reduction steps, criteria on the total charge collected by ID and the number of hits on OD 5/39 PTEP 2019, 053F01 M. Jiang et al. are applied to remove most cosmic ray muons. In the third reduction step, a dedicated cosmic ray muon reconstruction algorithm is used to apply a cut on OD hits near the entrance or exit point of a reconstructed track. The fourth step is designed to remove so-called "flasher" events caused by internal corona discharge of a PMT. Failing tubes often produce multiple events with similar topologies as these discharges repeat. Typical flasher events have a broad distribution in time among hit PMTs and spatially similar hit distributions. Cuts based on the hit timing distribution and the charge pattern correlation with other events that passed the third step reduction are applied to remove flasher events. The fifth reduction rejects the remaining cosmic ray muons, flasher events, and electronic noise further using a more precise fitter. Finally, the reconstructed vertex is required to be within a specified fiducial volume to further reduce non-neutrino background events and to prevent the reconstruction performance from deteriorating near the detector wall. These steps are particularly relevant to the new reconstruction algorithm considered here and will be discussed in Sect. 6 in detail. An event's visible energy (E vis ), which is defined as the energy of an electromagnetic shower producing the same amount of Cherenkov light as observed in the event, is also required to be larger than 30 MeV to remove low-energy backgrounds. All events passing the reduction are scanned by eye to determine the level of background contamination in the final analysis sample and to estimate the uncertainty inherent in the reduction process.
After reduction, the FC data are subdivided based upon the number of observed Cherenkov rings, the particle ID (PID) of the most energetic ring, and the number of observed electrons from muon decays into combinations of single-or multi-ring, electron-like (e-like) or muon-like (μ-like), and sub-GeV (E vis < 1330.0 MeV) or multi-GeV (E vis > 1330.0 MeV) categories, as well as the number of decay electrons (0, 1, 2 or more). In the present study, all FC events have been reconstructed using the new algorithm, which will be introduced in the next section.
The PC and Up-μ samples also have their own reduction processes, which are optimized for the topology of OD activity of neutrino events. These two samples are reconstructed using the preexisting algorithm and are divided into "stopping" and "through-going" subsamples based on the estimated muon stopping point for PC and Up-μ events. The "through-going" events in the Upμ sample are further divided into "showering" and "non-showering" based on whether the event induces an electromagnetic shower while traversing the ID. Although Super-K also works as the far detector for the T2K experiment [24] and detects accelerator-generated neutrinos, they are excluded based on event time information in this analysis. After all the selections there is a total of 13 FC analysis samples, two PC samples and three Up-μ samples.

FC event reconstruction
The Super-K event reconstruction algorithms determine an event's physical properties such as the interaction vertex, number of particles, the particle types and momenta based on PMT hit information.
The conventional event reconstruction algorithm, APFit, was introduced at the beginning of Super-K, contributing to the discovery of atmospheric neutrino oscillation and has been used in both the K2K and T2K experiments. APFit is a single-iteration fitter based on the time and charge information of hit PMTs [25]. The interaction vertex is reconstructed based on hit timing information after accounting for the photon time of flight. Then the direction of the first ring found, usually also the most energetic ring, is determined based on the observed charge distribution with respect to the interaction vertex. Additional ring candidates are found using a Hough-transform-based method and selected by a likelihood function optimized to reject spurious ring candidates. After determining the number of Cherenkov rings in the event, the particle type of each ring is determined based upon the Cherenkov ring pattern and opening angle. Rings from electrons tend to have rough edges produced by the light from their electromagnetic showers, while rings from muons or charged pions predominantly produce crisp edges. In the last step of the fit, the momentum of each ring is evaluated based on the observed charge inside a cone with a half angle of 70 • drawn along the line connecting the interaction vertex to each ring's center. Corrections and adjustments are made to account for charge sharing between overlapping rings. It is worth noting that the hit time information is only used during the first step to find the initial vertex candidate in APFit.
The new algorithm, named fiTQun, employs a maximum likelihood method to reconstruct particle types and determine kinematics in the detector simultaneously. The algorithm is based on methods developed for the MiniBooNE experiment [26], but has been developed from scratch for Super-K with additional features such as multi-ring reconstruction for events with multiple final-state particles. Compared to APFit, fiTQun uses more information, including information from PMT hits outside of the expected Cherenkov cone and hit timing information, during the fitting procedure. For a given event fiTQun's fit procedure will run multiple times to determine the best kinematic parameters for each possible particle configuration hypothesis, while APFit only fits those parameters once. The remarkable evolution of computing power since the start of Super-Kamiokande has enabled fiTQun to achieve higher reconstruction precision on similar timescales as APFit used to be. FiTQun has already been used in the T2K analyses using an expanded fiducial volume due to its improved resolution of reconstructed quantities and particle identification [27].

Likelihood function
An event topology hypothesis (e.g., single-ring e-like) together with its associated kinematic parameters θ , which include the vertex position, particle creation times, the azimuthal and zenith angles of the particle directions, as well as their momenta, are considered in the likelihood function during a fit. Based on the observed charge and hit time of each PMT, fiTQun constructs the following likelihood function for a given hypothesis to estimate the kinematic variables: In this equation, the index j runs over all PMTs that did not register a hit ("unhit" PMTs) and for each of these the probability that it does not register a hit given the fitting hypothesis ( , θ) is calculated as P j (unhit| , θ). PMTs that did register a hit are indexed with i. For such PMTs the likelihood density for observing a charge q i under the fitting hypothesis is represented by the charge likelihood f q (q i | , θ). The likelihood density of producing a hit at the observed time t i is defined similarly as f t (t i | , θ). Since the processes of particle and optical photon propagation are decoupled from the response of the PMT and the electronics, the likelihood can be rewritten in terms of the expected number of photoelectrons produced at the ith PMT given the hypothesis (the predicted charge), μ i ( , θ), as P j (unhit|μ j ) and f q (q i |μ i ) are properties of the PMT and the electronics, and therefore do not explicitly depend on the process of Cherenkov photon emission and propagation in water.
In the calculation of the predicted charge, the contributions from direct light and light that has scattered or been reflected (indirect light) are considered separately and summed to form the final μ i . The predicted charge from direct light reaching a PMT is calculated by integrating the Cherenkov emission profile along the track while correcting for the distance from the track to the PMT, the light transmission in water, and the PMT angular acceptance. Charge produced by indirect light reaching the PMT is predicted by integrating the product of the direct light emission profile and a scattering function that has been generated in advance based on simulation and incorporates effects arising from the relative position of the PMT and light source. For events with multiple Cherenkov rings the predicted charge of each ring is first calculated separately and then the sum from all the rings is used to calculate the total expected μ i . The final charge likelihood f q (q i |μ i ) is obtained by comparing the observed charge in a PMT against the prediction assuming photoelectrons generated according to Poisson statistics.
The time likelihood term can be expressed as f t (t i |t exp i , , p, μ i ), where p is the momentum for the topology and t exp i is the expected hit time. The latter is defined as the arrival time of unscattered photons emitted at the track midpoint and traveling directly to the PMT as where x and t are the vertex and creation time of the particle, respectively, d is the particle direction, R PMT i is the position of the ith PMT, and s mid represents half of the track length. Here c and c n represent the light velocity in vacuum and in water, respectively. The time likelihood depends on the predicted charge since a hit is recorded by the first photon arriving at a PMT, which leads to a narrower distribution of hit time for higher numbers of incident photons and hence for more predicted charge. The track length of a particle s, which is determined by the topology and momentum p, also affects the shape of the time likelihood since not all photons are generated at the track midpoint. Contributions to the likelihood from direct and indirect photon hits are calculated separately and then merged according to their relative intensities to obtain the final time likelihood. The time likelihoods are determined using particle gun simulations. For multi-particle hypotheses, the time likelihood is calculated ring-by-ring and then merged to a final likelihood function assuming that the photons from a particle with earlier t exp i always arrive earlier than the photons from any other particles with later t exp i values. Once the likelihood function is defined, the best set of kinematic parametersθ for a given event topology hypothesis is defined as that which maximizes L( , θ). The best estimate for the particle content of a given event is determined by comparison of L( ,θ) among all hypotheses, .

Fitting procedure
The fiTQun reconstruction process can be divided into four steps. The first step, vertex pre-fitting, roughly estimates the interaction vertex based on the PMT timing information. During the second step, clusters of PMT hits in time are identified as candidates for particle activities. Thereafter, the single-ring reconstruction, which performs fits assuming event topologies with only a single light-producting particle, and multi-ring reconstruction, which fits using hypotheses with multiple particles, are performed in sequence. During these fits the negative log likelihood − ln L( , θ) is minimized with respect to θ, using the MINUIT [28] package.

Vertex pre-fitting
The vertex pre-fitter is a fast algorithm that uses only the hit time information from PMTs around the primary event trigger to estimate an initial vertex position assuming all observed light was emitted from a common point. This is done by searching for the vertex position x and time t that maximize the goodness function where The pre-fitter is executed several times while gradually shrinking the grid size and σ to achieve high precision efficiently. The fitted vertex from this step, called the pre-fit vertex, is just a rough estimation and will be fitted again with higher precision during minimization of L( , θ).

Hit clustering
Events in Super-K are defined by detector activity in an O(10 μs) time window around an event trigger but may contain multiple subevents representing clusters of PMT hits separated in time from the primary trigger. A common example is muon decay in which the muon produces the primary event trigger and its decay produces additional delayed detector hits that form a subevent. A hit clustering algorithm is used to search for activity around the primary trigger and locates any additional subevents for further fitting with the more precise reconstruction methods discussed below.
The hit clustering algorithm starts by searching in time for subevent activity around the event trigger using a peak-finding algorithm. The vertex goodness G(x, t) from Eq. 7 is scanned in t with the vertex position fixed to the pre-fit vertex to search for additional peaks from particle activity in the detector. An example of the goodness distribution from a muon decay event is shown in Fig. 1. The two dominant peaks in the distribution correspond to the parent muon and its Michel electron, respectively. To avoid counting peaks created by scattering or reflection processes within a cluster of hit activity, peaks are required to be above a minimum threshold F(t), which is defined as: where M represents all local maxima of goodness function G(x, t). The value of the time constant γ is 25 ns when t < t i and 70 ns otherwise. An offset η = 9 is added to the threshold function to suppress the effect of dark hits. The function F(t) is the blue curve in Fig. 1. The minimum goodness between any two peaks must be lower than a second threshold, 0.6 × F(t), which is shown as a green dashed curve. Under these criteria only the peaks labeled with a triangle are identified as candidates for further fitting. During this procedure all vertex positions are assumed to lie close to the pre-fit vertex when the peak-finding algorithm runs. This assumption is broken when the primary particle travels a significant 9  distance from the interaction vertex, as for high-momenta muons. Therefore, the vertex pre-fitting and peak-finding algorithm are rerun after masking the hits caused by the primary particle to improve decay electron reconstruction efficiency. The vertex position used in the goodness function x is then close to the vertex of the secondary particle.
A time window defined as −180 ns < T res i < 800 ns around each peak found is defined to contain its associated hits. Afterwards the vertex pre-fitter and peak finder will be run once again in each time window using only the hits within each window. Peaks remaining after this step become the final candidates for full event reconstruction.

Single-ring reconstruction
The most basic reconstruction that is applied to the time windows defined in the previous step is the single-ring fitter. This fitter poses single-particle hypotheses for the likelihood function in Eq. 4. Three types of single-ring hypothesis are considered in fiTQun, that of an electron, a muon, and a charged pion. For each hypothesis the kinematic parameters of the event, including the interaction vertex, the particle momentum, and its direction are varied to maximize the likelihood function against the observation. Particle identification (PID) is based on the best-fit likelihood values for each of these hypotheses. Electrons and muons, for instance, are separated by cutting on ln(L e /L μ ), the logarithm of the likelihood ratio between the best-fit electron and muon hypotheses. The distribution of this variable is shown for the FC atmospheric neutrino data and Monte Carlo (MC) simulation result with sub-GeV and multi-GeV energies in Fig. 2. In both cases there is a clear separation of the likelihood variable between electron-like (e-like) and muon-like (μ-like) events. Figure 3 shows the vertex resolution as a function of visible energy for the FC single-ring chargedcurrent quasi-elastic (CCQE) event sample in the atmospheric neutrino MC whose true interaction vertex lies within the conventional fiducial volume, the region located more than 2 m from the ID wall. As the visible energy increases from 100 MeV to 1330 MeV the vertex resolution of fiTQun for CCQE ν e events is stable at 20.6 cm. For the CCQE ν μ events, the vertex resolution of fiTQun,   which changes from 29.2 cm to 15.9 cm, is better than that of APFit. In the same energy range, the momentum resolution of fiTQun for CCQE ν e events improves from 5.39% to 2.58% as the visible energy increases, as shown in the top-left plot of Fig. 4. For CCQE ν μ events the momentum resolution is constant across the energy range, being lower than 2.5% for fiTQun and slightly worse for APFit. Similarly, fiTQun shows an improved ability to discriminate between electrons and muons with less than a 1% misidentification rate for visible energies less than 1330 MeV (Fig. 5). Other metrics used to measure fit quality are generally stable and their typical values can be found in Table 1. In general fiTQun performs as well as or better than APFit.

Multi-ring reconstruction
In atmospheric neutrino analyses it is essential to reconstruct events with multiple light-producing particles since a large fraction of events with multi-GeV energies have multi-particle final states. To save computing time, the fiTQun multi-ring fitter is applied only to the time window around the primary event trigger and not to any additional subevents. The process of reconstructing a multi-ring event starts by performing an iterative search for an additional ring on top of any existing rings from a previous fit result. The hypothesis in the likelihood from Eq. 5 is updated to include a new ring and the likelihood is minimized again allowing the kinematic parameters for the rings to vary. Three hypotheses, including e-like, μ-like, and π + -like rings, are tested. The result of the updated fit for each of the particle hypotheses is compared to the original result to determine the validity of the added ring. If the likelihood improves with the new ring it is accepted and the processes to search for further rings is repeated. The cycle of adding, fitting, and examining new rings iterates until either a newly added ring fails the likelihood criterion or six rings are found in the event. Figure 6 shows the difference in the best-fit likelihoods from the two-ring and one-ring hypotheses. Based on an optimization performed in MC, a cut at 9.35 (11.83) on this likelihood ratio is used to separate single-  and multi-ring events when the first ring hypothesis is e-like (μ-like). However, ring candidates that have an angular separation of less than 20 • from the most energetic ring are discarded as spurious since MC studies indicate that such candidates are typically due to particle scattering rather than a new particle. When this occurs the two rings are merged and refitted as one for all particle hypotheses while keeping all other rings fixed. This procedure is repeated for all rings in the event in descending order of their energy. Table 2 shows the performance of the ring counting in both APFit and fiTQun using atmospheric neutrino MC events. In this comparison only charged particles with energies more than 30 MeV above the Cherenkov threshold in the final states without any requirement on their angular separation with other particles are considered as true ring candidates. Compared to APFit, fiTQun shows a greater ability to reconstruct multi-ring events, while the ability to correctly identify single-ring topologies is the same in both algorithms. Although fiTQun tends to have more fake rings than APFit, the improvements largely outweigh this slight downside.      Figure 7 shows the PID likelihood variable distribution of the most energetic ring in fully contained multi-ring events. Due to the overlap of Cherenkov photons from multiple particles, the separation between e-like and μ-like is not as good as that for single-ring events.
The performance of the fitter's kinematic reconstruction on multi-ring events is checked by studying the invariant mass of the π 0 in charged-current single-pion (CC1π 0 ) events. Figure 8 shows the reconstructed invariant mass calculated using the second and third most energetic rings for events with three e-like rings (eee events) for both reconstruction algorithms. Shaded histograms show true CCν e 1π 0 events. Both fitters have a peak in the invariant mass distribution near the π 0 mass, but the fiTQun peak is larger, sharper, and has fewer backgrounds, indicating improved reconstruction efficiency over APFit. The FWHM of the π 0 mass peak selected by fiTQun is 38 MeV, while this value is 42 MeV for APFit. The event rate and the purity of the CCν e 1π 0 events with 85 < m < 185 MeV (indicated by the arrows in the figure) are 278.5 and 52.1% for fiTQun, which are much higher than the values for APFit: 175.9 and 42.9%. Similarly, CCν μ 1π 0 events from the three-ring sample whose leading ring is μ-like and passing this selection are also studied and summarized in Table 3. The fiTQun algorithm demonstrates a higher efficiency and purity for multi-ring event reconstruction. 15   Another important performance indicator for the multi-ring event reconstruction algorithm is the separation of the neutrino and antineutrino components of the atmospheric neutrino sample, since at multi-GeV energies they have the most sensitivity to the mass hierarchy. A two-stage likelihood method has been developed to purify the neutrino and antineutrino components for multi-ring events. The first stage of the separation is designed to extract CC ν e +ν e interactions based on a likelihood selection as in Ref. [29] using the APFit algorithm. Four variables, including the visible energy fraction of the most energetic ring, the visible energy fraction of the most energetic charged pionlike ring, the number of decay electrons, and the distance to the farthest decay electron normalized by the total visible energy, are used as the inputs to the likelihood function. However, in the present study the inputs have been replaced by the equivalent variables from the fiTQun reconstruction. Events that pass this selection are classified as "multi-ring e-like" while those that fail are termed "multi-ring other" as shown in Fig. 9. Both are used in the oscillation analysis discussed below.  Neutrino and antineutrino interactions are separated from the multi-ring e-like sample during the second stage of the selection. The method is similar to the one discussed in Ref. [6] and uses the number of decay electrons, the number of reconstructed rings, and the event's transverse momentum as reconstructed by fiTQun. Figure 10 shows the final likelihood distribution used in SK-IV. The efficiency and purity for selecting e-like events and identifying true CCν e (ν e ) events asν e -like (ν e -like) are summarized in Table 4. FiTQun shows improved performance in each of these metrics that will translate to better sensitivity to sin 2 θ 13 and the mass hierarchy.

Fiducial volume expansion
In previous analyses using the APFit algorithm the fiducial volume (FV) was defined as the region located 200 cm from the wall of the inner detector. With this definition approximately 30% of the Super-K inner detector mass is not used for atmospheric neutrino analysis. On the other hand, the sensitivity to several open questions is currently limited by a lack of statistics. The analysis sensitivity can be increased by taking advantage of fiTQun's improved reconstruction performance to expand the fiducial volume while keeping signal purities high and backgrounds low.
The event selection and categorization in this study are similar to those one used in Ref. [6], but all cut variables used to define the FC samples are reconstructed by fiTQun. Partially contained and Up-μ events, on the other hand, are defined using APFit in the same manner as in the previous publications.

MC study for fiducial volume expansion
Though a larger fiducial volume provides higher statistics, the purity of signal interactions may decrease if the reconstruction performance deteriorates for events close to the ID wall. Events that are too close to the wall and producing particles that travel toward it will only illuminate a few PMTs and as a result will have poorly imaged Cherenkov rings. Further, as the FV cut is moved closer to the wall, the possibility of events with interaction vertices outside of the ID wall being reconstructed within the FV increases. Although the OD is designed to veto such events, as introduced in Sect. 4, events that originate in the 55 cm dead region between the ID and OD optical barriers or interact very near the OD wall may not produce enough light to trigger the OD cut and may remain in the 17  analysis sample. Such "entering events" are a type of misreconstructed background in the oscillation analysis with potentially large systematic errors. Figure 11 shows the purity of several components of the FC sub-GeV μ-like 0 decay-e sample as a function of the cut on the reconstructed distance from the event vertex to the nearest ID wall (dwall) as an example. The fraction of signal events, ν μ CC andν μ CC interactions (red) for this sample, remains stable as the cut position is decreased to smaller values until it reaches 50 cm, where the background from entering ν events increases dramatically. All 13 FC event categories exhibit similar behavior, indicating that a FV cut at 50 cm is an acceptable limit due to entering ν events. The purity of the 13 FC samples defined using this cut and the conventional FV definition are summarized in Table 5 (cf. Table II of Ref. [6]).
An MC study was performed to investigate the effect of decreases in purity against improved statistical power of an expanded FV. Figure 12 shows the sensitivity to the normal mass hierarchy assuming different FV cut values. Here systematic errors are assumed to be the same as those from the standard SK analysis, which uses a cut at 200 cm. The rise in sensitivity as the FV is expanded indicates that an FV cut at 50 cm can be expected to provide improved hierarchy sensitivity as long as systematic errors in the detector region outside of the conventional cut boundary are stable. A more precise sensitivity study with updated systematic errors for this region is presented below.

Background estimation
In addition to entering backgrounds and events with incorrectly assigned PID, non-neutrino (non-ν) backgrounds from cosmic ray muons and flashing PMTs are potential background sources when the FV is expanded.
Cosmic ray muons must traverse the OD before reaching the ID and are therefore mostly rejected by cuts on activity in the OD. Any remaining events are rejected by the FV cut, since their true vertex position, i.e., the entrance point of the muon, is expected to be on the wall of the tank. However, when the FV cut is near the wall cosmic ray background events might be introduced into the FC sample due to the vertex resolution of the reconstruction algorithm.  Flasher events are caused by the internal electrical discharge of a PMT and represent another kind of non-neutrino background. The reconstructed vertex for such events is also expected to be on the wall of the detector, so, as was the case for cosmic ray muons, more flasher backgrounds might exist in the FC sample when the FV is enlarged.
In order to observe possible excesses in the data due to such backgrounds, which are not modeled in the atmospheric neutrino MC, the distribution of the distance from the reconstructed vertex to the nearest wall of FC events is compared to MC in Fig. 13. Though not shown here, the distributions of the events' momentum, direction, and particle ID have been checked for both the conventional FV region (dwall > 200 cm) and new region (50 cm < dwall < 200 cm), respectively. No evidence for the presence of a large unmodeled component to the data is observed. 19   To evaluate the background rate more precisely all FC events within the fiducial volume were scanned by eye using a graphical event display program. In total, one flasher event and 24 cosmic ray muon events were found in the FC sample with dwall > 50 cm. Most of the cosmic ray muon backgrounds were reconstructed as multi-ring downward-going μ-like events with more than 5 GeV/c of momentum. Eye scans of these events indicated that their vertices were incorrectly reconstructed within the new fiducial volume due to the presence of an erroneous second ring found by the fitter. Both the rate and type of such misreconstructions are consistent with those found in independent 20   cosmic ray muon samples. Using the dwall distribution of the cosmic ray muon events, the expected cosmic ray background in the final analysis sample is statistically removed. This is done under the assumption that the FC event reduction process does not depend strongly on the muon kinematics. By comparison of the event classification of this control sample with that observed in the eyescanned background, this assumption has been validated for all particle directions and momenta. As an example, Fig. 14 shows the dwall distribution for multi-ring μ-like events. More than 20 000 cosmic ray muon events (solid line) are shown here and their dwall distribution has been normalized to the eye-scanned background events (points) in the FC sample. The shapes of these two samples are in good agreement. Due to their limited contamination in the final FC sample, the effect of flasher PMTs is treated with a systematic error on their rate in the oscillation analysis presented below. 21    Expanding the FV opens the possibility that events interacting within the ID but close to its wall may produce particles that escape into the OD. Poor modeling of the response of the OD can thereby potentially introduce biases and relative inefficiencies in the FC sample. Figure 15 shows the distribution of OD hits used to define the FC sample after the reduction processes for data with dwall > 50 cm. The detection inefficiency due to the cut on OD hits for both data and MC are confirmed to be consistent, around 0.2%. Based on this result in conjunction with the stability of the reconstruction algorithm, its systematic errors, and sample purities across the detector, the fiducial volume definition is expanded from its conventional value to dwall > 50 cm (expanded FV) in the analysis presented below.
Zenith angle distributions for each analysis sample using the expanded FV are shown in Fig. 16. Their event rates since the start of SK-IV have been stable at 8.3 (2.2) FC events per day in the conventional FV (new region), 0.73 PC events per day, and 1.49 Up-μ events per day, as shown in Fig. 17.

Systematic error
While systematic errors related to the atmospheric neutrino flux and cross section model are the same as those used in Ref. [6], systematic errors on the event selection have been updated for FC events using the expanded fiducial volume outlined above. Since fiTQun is not used to reconstruct PC or Up-μ events in the present work, the errors for those samples are the same as in the previous publications.
The uncertainty on the absolute energy scale is estimated by comparison of data and MC across three calibration samples spanning momenta up to 10 GeV/c: electrons from cosmic ray muon decay, atmospheric neutrino events producing single π 0 s from neutral current interactions, and stopping cosmic ray muons. The difference between data and MC for the calibration samples reconstructed within the expanded FV is shown in Fig. 18. The calculation of the total systematic error follows that of Ref. [6] and is the sum in quadrature of the absolute energy scale error, which is defined as the largest data-MC difference across all samples, and momenta with the average of the time variation of these samples throughout SK-IV. The uncertainty from the up/down asymmetry of the detector is 23  measured with Michel electrons from cosmic ray muons by comparison of their momenta for data and MC as a function of zenith angle. The difference between data and MC at the most deviated direction is taken as the uncertainty from the detector asymmetry. As summarized in Table 6, the uncertainty of the energy scale on the conventional (expanded) FV is 2.17% (2.02%) using fiTQun reconstruction, which is similar to the value with APFit: 2.1% [6].
Lacking other control samples that span the same energies and event topologies as the atmospheric neutrino data, systematic errors on the event selection including the particle identification (PID), estimated number of rings (ring counting), and the two-stage separation of the multi-GeV multi-ring e-like event samples described above are evaluated using the atmospheric neutrino data by comparing the likelihood distributions of the atmospheric neutrino data and MC. The MC distributions are shifted and smeared during a fit to data distribution. The fractional change in the MC after the cut used in the likelihood selection before and after fitting is defined as the systematic error.
Systematic errors are estimated for each of the FC analysis samples for both the conventional fiducial volume and the additional region formed between the conventional and expanded FV boundaries. The results are summarized in Table 7. Errors for the individual samples subject to a particular uncertainty, such as the particle ID for single-ring events, are assumed to be fully correlated (or fully anti-correlated). Anti-correlations are indicated by negative numbers in the table. No correlation is assumed between the systematic error categories. The sizes of the systematic errors evaluated using fiTQun and APFit are consistent within the conventional fiducial volume, with improvements seen in fiTQun for many samples. For many error categories the fiTQun reconstruction shows consistent errors between the conventional fiducial volume and the new region (50 cm < dwall < 200 cm). In some cases however, such as the ring counting uncertainty for multi-GeV single-ring events, larger systematic errors are observed. At present these are not large enough to offset the benefit to the oscillation sensitivity from including events in the new region in the analyses. A quantitative estimation of the sensitivity is presented in the next section. The numbers for the different FV regions are merged to form the final response function used in the analysis, which will be described in the next section. Systematic errors and their sizes at the best-fit point of the analysis are presented in Table A1.

Atmospheric neutrino oscillation analysis
The Super-K atmospheric neutrino data are separately fitted against both the normal and inverted hierarchy hypotheses. Since θ 13 has been measured and is well constrained by reactor experiments, two types of analyses are done: one with θ 13 as a free fitting parameter and the other with it constrained 24 Table 7. Summary of systematic errors related to the ring counting, particle identification, and multi-GeV multi-ring e-like separation. The sign of the number represents the correlation among samples, with negative numbers denoting fully anti-correlated samples.  to 0.0210±0.0011 [3]. Though the constraint on this parameter is much stronger than can be expected from Super-K alone, the appearance of upward-going electron events characteristic of the mass hierarchy signal is driven by θ 13 and therefore measurements with atmospheric neutrinos represent an important test of the analysis' hierarchy preference. The solar mixing parameters, m 2 21 and sin 2 θ 12 , on the other hand, have little impact on the oscillations of neutrinos in the energy range considered here and will be fixed in the analysis as described below.   Table 9. Values of oscillation parameters fixed in the analysis and their systematic errors. Note that sin 2 θ 13 is only fixed in the "θ 13 constrained"analyses described in Sect. 8.2.

Particle ID of the brightest ring in multi-ring events
Parameter Value m 2 21 (7.53 ± 0.18) × 10 −5 eV 2 sin 2 θ 12 0.307 ± 0.013 sin 2 θ 13 0.0210 ± 0.0011 Data are fitted to the MC expectation with a binned χ 2 method assuming Poisson statistics. The effects of systematic errors are regarded as scaling factors on the MC in each bin [30], with different error sources assumed to be independent. The definition of χ 2 used in the fits is where In this equation O n is the observed number of events in the nth analysis bin. Similarly, E n,0 represents the nominal MC expectation in that bin and E n is the expectation after incorporating the effect of the systematic errors. The coefficient f i n describes the fractional change in the bin's MC under a 1σ i variation of the ith systematic error source. The sum over ( i /σ i ) 2 penalizes the χ 2 for adjusting the systematic errors when bringing the MC into agreement with the data. For each set of oscillation parameters tested, the MC are fitted to the data by minimizing χ 2 over the error parameters i . All fits are performed over 515 analysis bins in SK-IV.

Analysis with unconstrained θ 13
The fit with θ 13 unconstrained by reactor measurements is performed over four parameters, | m 2 32,31 |, sin 2 θ 13 , sin 2 θ 23 , and δ CP with a total of 71 systematic errors. Fits to the normal hierarchy use m 2 32 while those for the inverted hierarchy use m 2 31 . The agreement between the data and MC is evaluated using Eq. 10 at each point in the grid of points shown in Table 8, while keeping the solar mixing parameters, m 2 12 and sin 2 θ 12 , fixed to the values in Table 9. Systematic errors representing the uncertainty in these parameters are included as Eq. 10 in the analysis. All systematic errors have been profiled over in the figures presented below. The best-fit parameters for each hierarchy are defined as the parameter set returning the smallest χ 2 value. Between the two hierarchy fits the one with smallest minimum χ 2 value is taken as the global best fit and is taken to indicate the mass hierarchy preference. 26 . 20. Constraints on neutrino oscillation parameters from the Super-K atmospheric neutrino data with no assumed constraint on θ 13 . The solid blue and dashed orange lines denote the normal and inverted hierarchy fit results, respectively. The latter has been offset from the former by the difference in their minimum χ 2 values.
Results and discussion Figure 19 shows the expected lower limit of sin 2 θ 13 at 1σ as a function of the true value of sin 2 θ 23 for the fiTQun-based analysis samples with both the conventional and expanded FV as well as an analysis with APFit-reconstructed samples with the conventional FV. The sensitivities for the two algorithms using the same conventional FV are similar since the fraction of events common between them is larger than 97%. The expected numbers of atmospheric neutrino events selected by APFit PTEP 2019, 053F01 M. Jiang et al. Table 11. Parameter estimates for each analysis and hierarchy hypothesis considered. Here NH (IH) refers to the normal (inverted) hierarchy fit. The terms "free" and "constrained" refer to fits without and with a constraint on sin 2 θ 13 , respectively. For sin 2 θ 23 parameter ranges are shown for both octants, with the best-fitting octant enclosed in a box. The expected absolute χ 2 value for the constrained fit is 523.6 and the probability for obtaining a larger value is 0.244 for NH. However, since the contamination of ν μ events in the e-like samples most sensitive to θ 13 is lower for the fiTQun-selected sample as shown in Table 10, there is a slight boost in its sensitivity. Expanding the FV naturally incorporates more signal events in the analysis and the sensitivity improvement is clear.
One dimensional allowed regions for θ 13 , | m 2 32,31 |, sin 2 θ 23 , and δ CP from the fit to the data are shown in Fig. 20. The normal hierarchy hypothesis yields a better data-MC agreement than the inverted hierarchy hypothesis with χ 2 NH,min − χ 2 IH,min = −1.81. The 1σ allowed region for sin 2 θ 13 is from 0.003 to 0.033 (from 0.001 to 0.023) for the normal (inverted) hierarchy fit, which is consistent with the globally preferred value. The point at sin 2 θ 13 = 0.0 is disfavored at approximately 1.8σ (1.2 σ ) for the normal (inverted) fit. A summary of the best-fit information and parameter constraints is presented in Table 11.
As discussed in Sect. 2, the determination of sin 2 θ 13 and the mass hierarchy using atmospheric neutrinos is achieved using the upward-going electron neutrino appearance at multi-GeV energies. Figure 21 shows the up-down asymmetry of the multi-GeV single-and multi-ring electron-like samples, where the asymmetry is defined as the ratio of the difference of the number of upward-going and downward-going events relative to their sum. Here upward-going (downward-going) events are events whose zenith angle satisfy cos θ z < −0.4 (cos θ z > 0.4). Excesses between a few and ten GeV in the multi-GeV e-like ν e and the multi-ring e-like ν e samples drive the normal hierarchy preference.
The best-fit atmospheric mixing parameters from the normal hierarchy fit are m 2 32 = 2.63 +0.10 −0.21 × 10 −3 eV 2 and sin 2 θ 23 = 0.588 +0.030 −0.062 for the second octant (best-fit) and 0.425 +0.051 −0.034 for the first octant. Maximal mixing (sin 2 θ 23 = 0.5) is weakly disfavored at around 1σ significance. Besides the data excesses in the upward-going regions of the single-ring e-like ν e sample and multi-ring other sample, the data deficits in the multi-GeV μ-like sample also contribute to this preference. It should be noted that the preference for sin 2 θ 23 is coupled to the θ 13 measurement, since both parameters feature in the ν μ → ν e oscillation probability. The jagged nature of the | m 2 32 | curve is due to the limited statistics at low energies in both data and MC. In this region, the oscillation probability is  roughly proportional to sin 2 (1.27 m 2 32 L/E) and changes wildly for small changes in the neutrino energy, E, or the neutrino travel length, L. Since there is no sufficient statistics to fully sample and average over these oscillations, a rapid change in the χ 2 is observed in the figure for small changes in m 2 32 . The best-fit value of δ CP is 3.84 radians for both the normal and inverted hierarchy fits. Comparing with the least preferred parameter value of 0.8 radians, more electron neutrino appearances in the sub-GeV e-like samples are expected to be observed due to ν μ → ν e oscillations. In the multi-GeV region, δ CP similarly modulates the amount of electron neutrino appearance but it is subdominant to the effects of θ 13 . In the inverted hierarchy fit, the data are less sensitive to θ 13 and δ CP since only antineutrino events experience matter effects in the Earth while there are more neutrino events than antineutrino events.

Analysis with constrained θ 13
Since atmospheric neutrinos come from all directions, carry a wide range of energies, and often produce particles that are invisible to Super-K it is not possible in general to fully reconstruct the neutrino kinematics. As a result atmospheric neutrinos themselves do not typically have oscillation parameter sensitivity at the same level as that in long-baseline or reactor experiments, such that the introduction of constraints from those measurements can improve sensitivity to the mass hierarchy and δ CP . Indeed, Eq. 2 shows that the size of the MSW resonance relies on the value of sin 2 θ 13 directly. Since reactor neutrino experiments constrain this parameter more precisely than the analysis in the 29   previous section, the analysis presented here is performed with the value of sin 2 θ 13 constrained to 0.0210 ± 0.0011 [3]. During the fit the parameter is held at its central value and a systematic error that imparts the impact of its uncertainty is added to the analysis. The fit uses the same data samples, binning, and parameter grid for all other oscillation parameters as the analysis presented above.

Results and discussion
The expected sensitivities to the mass hierarchy and θ 23 octant for the true normal hierarchy are illustrated in Figs. 22 and 23 for the fiTQun-reconstructed analysis samples using both the conventional and expanded FV and for the APFit-based analysis with the conventional FV. Improved sensitivity is found with the fiTQun-based analyses regardless of the fiducial volume used for the same reasons as produced the improved sensitivity to θ 13 in the previous section.    25. Constraints on the atmospheric mixing parameters using the SK-IV atmospheric neutrino and the expanded FV. The solid blue (dashed orange) line shows the 90% CL for the normal (inverted) hierarchy. The star denotes the best-fit value, which is at the same point for the normal and inverted hierarchies, as shown in Table 11. In each contour sin 2 θ 13 is constrained to be 0.0210 ± 0.0011. The contours have both been drawn relative to the global best fit. Figure 24 shows the χ 2 value as a function of the atmospheric neutrino mixing parameters and δ CP in the θ 13 -constrained fit. As in the unconstrained fit the data prefer the normal mass hierarchy with χ 2 = −2.45. The 1σ allowed region for | m 2 32 | is from 2.41-2.75 ×10 −3 eV 2 (2.36-2.67 ×10 −3 eV 2 ) for the normal (inverted) hierarchy fit, which is consistent with the result of the unconstrained fit. The contour for the allowed region with 90% CL is shown in Fig. 25. The result of the θ 13 -constrained fit shows a weak preference for the first octant of θ 23 , while it is in good agreement with other measurements, as shown in Fig. 26. Adding the constraint on sin 2 θ 13 , the preferred value of θ 23 shifts from the second to the first octant, though both are allowed at 1σ . Indeed, the difference in χ 2 between the minimum value of each octant is χ 2 1st−2nd = −0.73 (0.13) for the constrained (unconstrained) fit. This result can be explained by degeneracies in the appearance probability that arise for certain combinations of sin 2 2θ 13 and sin 2 θ 23 as shown in Eq. 1; an increase in one can be 31 26. Constraints on the atmospheric mixing parameters at the 90% CL. The fiTQun result (solid red) is taken from the analysis with the normal mass hierarchy hypothesis and sin 2 θ 13 constrained to be 0.0210 ± 0.0011. Constraints from the Super-K I-IV with APFit (dashed cyan) [6], T2K (dotted yellow) [4], NOvA (dashed green) [31], IceCube (dashed black) [32], and MINOS+ (dashed blue) [33] experiments are also shown. Fig. 27. Distributions of the difference in best-fit χ 2 values between first octant and second octant fits to pseudo data sets used in the generation of the CL O s value for the SK θ 13 constrained analysis. In the cyan (orange) histogram the pseudo data have been generated assuming sin 2 θ 23 = 0.4 (sin 2 θ 23 = 0.6). Shaded portions of the histograms denote the fraction of pseudo data sets with more extreme values than those observed in the data, χ 2 data = −0.73.
compensated by a decrease in the other to describe the data. The same upward-going event excesses (deficits) in the e-like (μ-like) samples that drive the hierarchy, θ 13 , and octant preferences in the unconstrained analysis are not strong enough to support the increased excesses (deficits) expected for higher values of θ 13 and the second octant. Since the θ 13 is now fixed, the data can be accommodated by reducing the expected number of e-like events and increasing the number of μ-like events by moving θ 23 to the first octant. The best-fit value of δ CP is 3.14 (4.89) for the normal (inverted) hierarchy hypothesis, with a tighter constraint compared with the unconstrained fit. Parameter values and their 1σ errors are summarized in Table 11.
The CL s method [34] is used to address the significance of the octant and mass hierarchy preferences observed in the data. As in the previous Super-K study [6] MC ensembles assuming different parameter combinations of the octant and hierarchy were generated with statistical and systematic error variations applied to each pseudo data set. Statistical variations have been applied assuming 32 28. Distributions of the difference in best-fit χ 2 values between normal and inverted hierarchy fits to pseudo data sets used in the generation of the CL H s value for the SK θ 13 constrained analysis. In the cyan (orange) histogram the pseudo data have been generated assuming the normal (inverted) hierarchy at the analysis best fit shown in Table 11. Shaded portions of the histograms denote the fraction of pseudo data sets with more extreme values than those observed in the data, χ 2 data = −2.45. Table 12. Summary of the probability of observing a χ 2 preference for the NH more extreme than that observed in the data assuming an IH, p 0 (IH), and CL H s values for a range of assumed parameters. The first row shows the true θ 23 used to generate MC ensembles used in the calculations. Other oscillation parameters are taken from the analysis best fit. The value of θ 13 is constrained to 0.0210 ± 0.0011.
where p 0 (O i ) represents the p-value with respect to obtaining a difference in the best-fit χ 2 values between both octant hypotheses smaller (larger) than the value of data assuming that the true value of sin 2 θ 23 is O 2 = 0.6 (O 1 = 0.4). These pseudo MC distributions of these hypotheses are shown as the orange and cyan shaded histograms in Fig. 27, respectively. The CL O s value is found to be just 0.098 and, accordingly, the data show no strong preference for the octant.
Similarly, the parameter for the mass hierarchy is defined as where p 0 (IH) and p 0 (NH) are p-values for obtaining a difference in the χ 2 values of the best-fit mass hierarchies more extreme than those of the data assuming a true IH and NH, respectively. Figure 28 shows the distribution for the mass hierarchy determination. Due to the large uncertainty on θ 23   the preference is 26.0%. A smaller value of sin 2 θ 23 predicts fewer electron neutrino appearances, which results in smaller p 0 (IH) and larger CL H s values due to the multi-GeV e-like sample event excesses seen in the data.

Conclusion
A new reconstruction algorithm based on a maximum likelihood method has been developed for Super-K. Compared to the conventional reconstruction algorithm, the new algorithm shows improved performance in a variety of metrics including event vertex resolution, particle momentum resolution, and particle identification. The new algorithm has further demonstrated reliable performance over a larger volume of the detector than the previous volume. Accordingly, in the present analysis the fiducial volume definition has been expanded to include all events reconstructed more than 50 cm from any ID wall. This represents a 32% increase in the number of events available for analysis relative to the 200 cm threshold used in previous SK analyses.
Using the new algorithm with its expanded fiducial volume definition an analysis of a 253.9 kton·year exposure of the SK-IV atmospheric data has yielded oscillation parameter estimates consistent with both previous Super-K measurements and results from other experiments. Assuming a normal mass hierarchy, the constraints on the atmospheric mixing parameters are sin 2 θ 23 = 0.425 +0.046 −0.037 (0.600 +0.013 −0.030 ) for the first (second) octant and m 2 32 = 2.53 +0.22 −0.12 , with δ CP = 3.14 +2.67 −1. 35 . The data show a weak preference for the normal mass hierarchy, disfavoring the inverted mass hierarchy at 74.0% assuming oscillation parameters at the analysis best-fit point. No strong preference for the θ 23 octant is observed. Figure 29 shows the expected sensitivity to mass hierarchy as a function of livetime for both reconstruction algorithms and their fiducial volumes. The expected improvement in sensitivity with the new algorithm becomes more apparent as data are accumulated, even assuming the conventional FV. While the new reconstruction has only been applied to the 3118.5 day SK-IV fully contained data set, future efforts will expand this study to include other Super-K samples and run periods, which constitute an additional 2800 days of livetime.

Appendix A. Systematic uncertainties
The systematic errors uses in the atmospheric neutrino analysis within the expanded FV and with a constrained θ 13 are summarized in Table A1. Table A1. Systematic error used in this analysis. The second column shows the best-fit value of the systematic error parameter in % and the third column shows the estimated 1σ error size in %. The fit result is within the expanded FV and with θ 13 constrained to 0.0210 ± 0.0011. The systematic errors in ring counting, particle identification, and multi-ring likelihood selection within the conventional fiducial volume (dwall > 200 cm) and within the new region (50 cm < dwall < 200 cm) have been merged.