Atmospheric Neutrino Oscillation Analysis With Improved Event Reconstruction in Super-Kamiokande IV

TL;DR

This paper presents an improved event reconstruction algorithm for Super-Kamiokande IV, enhancing neutrino detection and analysis, leading to a weak preference for the normal neutrino mass hierarchy based on 253.9 kton-year data.

Contribution

The paper introduces a new maximum likelihood-based reconstruction algorithm that improves event analysis and sensitivity to neutrino mass hierarchy in Super-Kamiokande IV.

Findings

Enhanced event reconstruction capabilities.

Increased sensitivity to neutrino mass hierarchy.

Weak preference for the normal hierarchy at 74% confidence.

Abstract

A new event reconstruction algorithm based on a maximum likelihood method has been developed for Super-Kamiokande. Its improved kinematic and particle identification capabilities enable the analysis of atmospheric neutrino data in a detector volume 32% larger than previous analyses and increases sensitivity to the neutrino mass hierarchy. Analysis of a 253.9 kton-year exposure of the Super-Kamiokande IV atmospheric neutrino data has yielded a weak preference for the normal hierarchy, disfavoring the inverted hierarchy at 74% assuming oscillations at the best fit of the analysis.

Tables

Table 1. Table 1: Summary of the basic performance of the APFit and fiTQun reconstruction algorithms on the fully contained CCQE single ring event sample with visible energy of 1 GeV.

Reconstruction	fiTQun	APFit
True CCQE ${\nu }_e$ sample
Vertex Resolution	20.64 cm	24.86 cm
Direction Resolution	1.48∘	1.68∘
Momentum Bias	0.43%	0.63%
Momentum Resolution	2.90%	3.56%
Mis-PID rate	0.02%	0.50%
True CCQE ${\nu }_{\mu }$ sample
Vertex Resolution	15.83 cm	17.32 cm
Direction Resolution	1.00∘	1.28∘
Momentum Bias	-0.18%	0.54%
Momentum Resolution	2.26%	2.60%
Mis-PID rate	0.05%	0.91%

Table 2. Table 2: Summary table of ring counting performance of both APFit and fiTQun on fully contained atmospheric neutrino events. Columns denote the number of reconstructed rings and rows the number of true rings. The true number of rings is counted using only particles with energy more than 30 MeV above Cherenkov threshold in the final states. The terms 1R, 2R, and 3R correspond to one ring, two rings, and three or more rings.

True Number	fiTQun Reconstruction			APFit Reconstruction
of Rings	1R	2R	$\ge$ 3R	1R	2R	$\ge$ 3R
True 1R	95.0%	4.64%	0.41%	95.9%	3.85%	0.29%
True 2R	27.8%	66.7%	5.56%	42.5%	52.8%	4.63%
True $\ge$3R	7.04%	25.5%	67.5%	20.2%	33.0%	46.8%

Table 3. Table 3: Summary table of performance of CC1 π 0 superscript 𝜋 0 \pi^{0} events selection in three ring events. The terms e ​ e ​ e 𝑒 𝑒 𝑒 eee and μ ​ e ​ e 𝜇 𝑒 𝑒 \mu ee represent the reconstructed PID of the three rings in order of momentum. The event rate column shows the number of events passing the invariant mass cut to select π 0 superscript 𝜋 0 \pi^{0} . The purity column shows the purity of target events in the passed invariant mass region.

Topology	fiTQun selection		APFit selection
	Event rate	Purity	Event rate	Purity
$𝒆𝒆𝒆$ events (target: ${\nu }_e$CC1${\pi }^0$)
Sub-GeV	278.5	52.1%	175.9	42.9%
Multi-GeV	112.1	47.0%	32.8	36.2%
$𝝁𝒆𝒆$ events (target: ${\nu }_{\mu }$CC1${\pi }^0$)
Sub-GeV	384.5	54.6%	201.2	38.0%
Multi-GeV	143.5	64.6%	51.6	32.0%

Table 4. Table 4: Summary of the Multi-ring e-like event selection. The efficiency represents the fraction of events passing the selection among all true charged-current ν e subscript 𝜈 𝑒 \nu_{e} events. Similarly, the purity is the fraction true CC ν e subscript 𝜈 𝑒 \nu_{e} out of all events passing the selection.

	fiTQun	APFit
First stage
Multi-ring $e$-like events
Efficiency	75.7%	69.7%
Purity	77.8%	69.5%
Second stage
Multi-ring True CC${\nu }_e$ events
Efficiency	56.8%	53.6%
Purity	58.8%	52.6%
Multi-ring True CC${\overline{\nu }}_e$ events
Efficiency	68.4%	70.9%
Purity	30.0%	25.9%

Table 5. Table 5: Sample purity broken down by neutrino flavor assuming neutrino oscillations with Δ ​ m 32 2 = 2.52 × 10 − 3 ​ eV 2 Δ subscript superscript 𝑚 2 32 2.52 superscript 10 3 superscript eV 2 \Delta m^{2}_{32}=2.52\times 10^{-3}\mbox{eV}^{2} and sin 2 ⁡ θ 23 = 0.51 superscript 2 subscript 𝜃 23 0.51 \sin^{2}\theta_{23}=0.51 with a FV cut of 200 cm. The number in the parentheses is the result within the new region, where the distance to nearest wall is between 50 cm and 200 cm.

Sample	CC ${\nu }_e$	CC ${\overline{\nu }}_e$	CC ${\nu }_{\mu }+{\overline{\nu }}_{\mu }$	CC ${\nu }_{\tau }$	NC
Fully Contained (FC) Sub-GeV
e-like, Single-ring
0 decay-e	0.728	0.242	0.001	0.000	0.028
	(0.702)	(0.227)	(0.025)	(0.001)	(0.045)
1 decay-e	0.907	0.020	0.033	0.001	0.040
	(0.712)	(0.015)	(0.208)	(0.001)	(0.063)
$\mu$-like, Single-ring
0 decay-e	0.010	0.004	0.795	0.001	0.191
	(0.034)	(0.011)	(0.805)	(0.001)	(0.150)
1 decay-e	0.000	0.000	0.974	0.000	0.026
	(0.001)	(0.000)	(0.968)	(0.000)	(0.031)
2 decay-e	0.000	0.000	0.984	0.000	0.016
	(0.000)	(0.000)	(0.980)	(0.000)	(0.019)
${\pi }^0$-like
Two-ring	0.051	0.016	0.011	0.000	0.922
	(0.109)	(0.036)	(0.018)	(0.000)	(0.837)
Fully Contained (FC) Multi-GeV
Single-ring
${\nu }_e$-like	0.726	0.077	0.058	0.027	0.113
	(0.748)	(0.066)	(0.064)	(0.016)	(0.105)
${\overline{\nu }}_e$-like	0.553	0.379	0.003	0.008	0.056
	(0.566)	(0.371)	(0.003)	(0.007)	(0.053)
$\mu$-like	0.000	0.000	0.996	0.003	0.001
	(0.000)	(0.000)	(0.995)	(0.004)	(0.001)
Multi-ring
${\nu }_e$-like	0.588	0.117	0.054	0.036	0.204
	(0.609)	(0.112)	(0.059)	(0.032)	(0.188)
${\overline{\nu }}_e$-like	0.526	0.300	0.021	0.020	0.134
	(0.541)	(0.301)	(0.023)	(0.016)	(0.118)
$\mu$-like	0.010	0.001	0.959	0.004	0.026
	(0.016)	(0.002)	(0.946)	(0.005)	(0.031)
Other	0.283	0.026	0.342	0.053	0.295
	(0.302)	(0.032)	(0.342)	(0.051)	(0.274)

Table 6. Table 6: Summary of uncertainty from the energy scale. The energy scale error is the quadratic-sum of absolute energy scale and the time variation. The absolute energy scale value is obtained from the most discrepant sample, i.e. Michel electrons from cosmic ray muons in this study.

Reconstruction	fiTQun	fiTQun	APFitAbe:2017osc
Dwall range	$>$50 cm	$>$ 200 cm	$>$ 200 cm
Energy calibration error	2.02%	2.17%	2.1%
Absolute energy scale	1.92%	2.09%
Time variation	0.62%	0.59%
Up/down asymmetry energy calibration error	0.67%	0.58%	0.4%

Table 7. Table 7: Summary of systematic errors related to the ring counting, particle identification and multi-GeV multi-ring e 𝑒 e -like separation. The sign of the number represents the correlation between these systematic errors, with negative numbers denoting fully anti-correlated samples.

Reconstruction	fiTQun	fiTQun	APFit
Dwall range	50 cm$\sim$200 cm	$>$ 200 cm	$>$ 200 cm
Ring counting
Sub-GeV , Single-ring
$e$-like, p 400 MeV	1.94%	1.20%	1.6%
$e$-like, p $>$ 400 MeV	0.59%	0.48%	1.0%
$\mu$-like, p 400 MeV	1.08%	0.42%	-3.0%
$\mu$-like, p $>$ 400 MeV	1.25%	1.21%	0.6%
Sub-GeV, Multi-ring
$e$-like	-3.58%	-2.39%	-1.9%
$\mu$-like, $E_{\text{vis}}>$ 600MeV	-2.32%	-1.66%	2.3%
Multi-GeV Single-ring
$e$-like	8.61%	1.21%	1.0%
$\mu$-like	2.65%	-2.33%	-1.2%
Multi-GeV Multi-ring
$e$-like	-4.21%	-0.62%	-0.9%
$\mu$-like	-3.07%	0.72%	2.4%
Particle ID of single-ring events
Sub-GeV
$e$-like	0.99%	0.36%	-0.28%
$\mu$-like	-0.89%	-0.37%	0.22%
Multi-GeV
$e$-like	0.23%	0.06%	-0.35%
$\mu$-like	-0.21%	-0.06%	0.35%
Particle ID of the brightest ring in multi-ring events
Sub-GeV
$e$-like	-3.19%	-0.72%	4.19%
$\mu$-like	1.31%	0.31%	-1.56%
Multi-GeV
$e$-like	1.94%	1.10%	3.33%
$\mu$-like	-1.06%	-0.66%	-1.56%
Multi-GeV Multi-ring $e$-like - other separation
$e$-like	-0.88%	-0.67%	3.0%
other	0.50%	0.53%	-3.4%
Multi-GeV Multi-ring ${\nu }_e$-like - ${\overline{\nu }}_e$-like separation
${\nu }_e$-like	-3.64%	-2.33%	6.82%
$\overline{{\nu }_e}$-like	4.51%	2.10%	-6.04%

Table 8. Table 8: Parameter grid used in the fit.

Parameter	Range	Scan points
$\left|\Delta m_{32,31}^2\right|$	$1.5\times {10}^{-3}\sim 3.5\times {10}^{-3}$	40 points
${\sin }^2{\theta }_{13}$	$0.0\sim 0.06$	25 points
${\sin }^2{\theta }_{23}$	$0.3\sim 0.725$	35 points
${\delta }_{CP}$	$0.0\sim 2\pi$	19 points

Table 9. Table 9: Values of oscillation parameters fixed in the analysis and their systematic errors. Note that sin 2 ⁡ θ 13 superscript 2 subscript 𝜃 13 \sin^{2}\theta_{13} is only fixed in the “ θ 13 subscript 𝜃 13 \theta_{13} constrained”analyses described in Section  8.2 .

Parameter	Value
$\Delta m_{21}^2$	$\left(7.53\pm 0.18\right)\times {10}^{-5}{\text{eV}}^2$
${\sin }^2{\theta }_{12}$	$0.307\pm 0.013$
${\sin }^2{\theta }_{13}$	$0.0210\pm 0.0011$

Table 10. Table 10: ν μ + ν ¯ μ subscript 𝜈 𝜇 subscript ¯ 𝜈 𝜇 \nu_{\mu}+\bar{\nu}_{\mu} CC backgrounds in multi-GeV e-like sample.

Reconstruction	fiTQun	fiTQun	APFit
Dwall range	50 $\sim$ 200 cm	$>200$ cm	$>$ 200 cm
Single-ring
${\nu }_e$-like	0.064	0.058	0.109
${\overline{\nu }}_e$-like	0.003	0.003	0.009
Multi-ring
${\nu }_e$-like	0.059	0.054	0.121
${\overline{\nu }}_e$-like	0.023	0.021	0.040

Table 11. Table 11: Parameter estimates for each analysis and hierarchy hypothesis considered. Here NH (IH) refers to the normal (inverted) hierarchy fit. The terms “Unconstrained” and “Constrained” refer to fits without and with a constraint on sin 2 ​ θ 13 superscript sin 2 subscript 𝜃 13 \mbox{sin}^{2}\theta_{13} , respectively. For sin 2 ⁡ θ 23 superscript 2 subscript 𝜃 23 \sin^{2}\theta_{23} parameter ranges are shown for both octants, with the best-fit octant enclosed in a box. The expected absolute χ 2 superscript 𝜒 2 \chi^{2} value for the constrained fit is 523.6 and the probability for obtaining larger value is 0.244 for the NH.

	${\theta }_{13}$ Free		${\theta }_{13}$ Constrained
Hierarchy	NH	IH	NH	IH
${\chi }^2$	576.3	578.1	576.5	579.0
${\text{sin}}^2{\theta }_{13}$	${0.008}_{-0.005}^{+0.025}$	${0.008}_{-0.007}^{+0.015}$	–	–
${\text{sin}}^2{\theta }_{23}$ ($1^{st}$ oct.)	${0.425}_{-0.034}^{+0.051}$	${0.425}_{-0.027}^{+0.075}$	\collectbox \BOXCONTENT ${0.425}_{-0.037}^{+0.046}$	\collectbox \BOXCONTENT ${0.425}_{-0.036}^{+0.055}$
${\text{sin}}^2{\theta }_{23}$ ($2^{nd}$ oct.)	\collectbox \BOXCONTENT ${0.588}_{-0.062}^{+0.030}$	\collectbox \BOXCONTENT ${0.575}_{-0.075}^{+0.034}$	${0.600}_{-0.030}^{+0.013}$	${0.588}_{-0.037}^{+0.022}$
$\left|\Delta m_{32,31}^2\right|$ [$\times {10}^{-3}$ eV2]	${2.63}_{-0.21}^{+0.10}$	${2.53}_{-0.08}^{+0.14}$	${2.53}_{-0.12}^{+0.22}$	${2.53}_{-0.31}^{+0.14}$
${\delta }_{CP}$	${3.84}_{-2.14}^{+2.00}$	${4.19}_{-4.19}^{+2.09}$	${3.14}_{-1.35}^{+2.67}$	${4.89}_{-3.46}^{+1.51}$

Table 12. Table 12: Normal hierarchy significance summarized in terms of the probability of observing a χ 2 superscript 𝜒 2 \chi^{2} preference for the NH more extreme than that observed in the data assuming an IH, p 0 ​ ( IH ) subscript 𝑝 0 IH p_{0}(\mbox{IH}) , and CL s H subscript superscript CL H 𝑠 \mbox{CL}^{\mbox{H}}_{s} values for a range of assumed parameters. The first row shows the true θ 23 subscript 𝜃 23 \theta_{23} used to generate MC ensembles used in the calculations. Other oscillation parameters are taken from the analysis’ best-fit. The value of θ 13 subscript 𝜃 13 \theta_{13} is constrained to 0.0210 ± plus-or-minus \pm 0.0011.

True sin${}_{}{}^2\theta _{23}^{}$	0.4	0.425	0.5	0.6
$p_0\left(\text{IH}\right)$	0.025	0.033	0.065	0.072
${\text{CL}}_s^{\text{H}}$	0.308	0.260	0.229	0.143

Table 13. Table 13: Systematic error used in this analysis. The second column shows the best-fit value of the systematic error parameter in percent and the third column shows the estimated 1 ​ σ 1 𝜎 1{\sigma} error size in percent. The fit result within expanded FV and θ 13 subscript 𝜃 13 \theta_{13} constrained to 0.0210 ± 0.0011 plus-or-minus 0.0210 0.0011 0.0210\pm 0.0011 . The same applies below. The systematic errors about ring counting, particle identification and Multi-ring likelihood selection within conventional fiducial volume (dwall > > 200 cm) and within new region (200 cm > > dwall > > 50 cm) have been merged.

Systematic Error			Fit Value (%)	$\sigma$ (%)
Flux normalization	111Uncertainty decreases linearly with $\log E_{\nu }$ from 25 %(0.1 GeV) to 7 %(1 GeV).		$7.3$	$25$
	$E_{\nu }>1\mathrm{GeV}$222Uncertainty is 7 % up to 10 GeV, linearly increases with $\log E_{\nu }$ from 7 %(10 GeV) to 12 %(100 GeV) and then to 20 %(1 TeV)		$9.6$	$15$
$\left({\nu }_{\mu }+{\overline{\nu }}_{\mu }\right)/\left({\nu }_e+{\overline{\nu }}_e\right)$			$1.3$	$2$
			$2.7$	$3$
	$E_{\nu }>10\mathrm{GeV}$333Uncertainty linearly increases with $\log E_{\nu }$ from 5 %(30 GeV) to 30 %(1 TeV).		$6.6$	$5$
${\overline{\nu }}_e/{\nu }_e$			$-1.1$	$5$
			$-2.7$	$5$
	$E_{\nu }>10\mathrm{GeV}$444Uncertainty linearly increases with $\log E_{\nu }$ from 8 %(100 GeV) to 20 %(1 TeV).		$-0.34$	$8$
${\overline{\nu }}_{\mu }/{\nu }_{\mu }$			$0.36$	$2$
			$3.24$	$6$
	$E_{\nu }>10\mathrm{GeV}$555Uncertainty linearly increases with $\log E_{\nu }$ from 6 %(50 GeV) to 40 %(1 TeV).		$9.2$	$15$
Up/down ratio		$e$-like	$0.079$	$0.1$
		$\mu$-like	$0.24$	$0.3$
		0-decay $\mu$-like	$0.87$	$1.1$
	$>400\mathrm{MeV}$	$e$-like	$0.63$	$0.8$
		$\mu$-like	$0.40$	$0.5$
		0-decay $\mu$-like	$1.34$	$1.7$
	Multi-GeV	$e$-like	$0.55$	$0.7$
		$\mu$-like	$0.16$	$0.2$
	Multi-ring Sub-GeV	$e$-like	$0.32$	$0.4$
		$\mu$-like	$0.16$	$0.2$
	Multi-ring Multi-GeV	$e$-like	$0.24$	$0.3$
		$\mu$-like	$0.16$	$0.2$
	PC		$0.16$	$0.2$
Horizontal/vertical ratio		$e$-like	$-0.023$	$0.1$
		$\mu$-like	$-0.023$	$0.1$
		0-decay $\mu$-like	$-0.069$	$0.3$
	$>400\mathrm{MeV}$	$e$-like	$-0.32$	$1.4$
		$\mu$-like	$-0.44$	$1.9$
		0-decay $\mu$-like	$-0.32$	$1.4$
	Multi-GeV	$e$-like	$-0.74$	$3.2$
		$\mu$-like	$-0.53$	$2.3$
	Multi-ring Sub-GeV	$e$-like	$-0.32$	$1.4$
		$\mu$-like	$-0.30$	$1.3$
	Multi-ring Multi-GeV	$e$-like	$-0.64$	$2.8$
		$\mu$-like	$-0.35$	$1.5$
	PC		$-0.39$	$1.7$
K/$\pi$ ratio in flux calculation666Uncertainty increases linearly from 5$\%$ to 20$\%$ between 100GeV and 1TeV.			$-5.3$	$10$
Neutrino path length			$-2.3$	$10$
Sample-by-sample	FC Multi-GeV		$-0.22$	$5$
	PC + Stopping UP-$\mu$		$2.46$	$5$
Matter effects			$0.06$	$6.8$
Solar Activity			$0.31$	$7$
$M_A$ in QE			$-4.79$	$10$
Single $\pi$ Production, Axial Coupling			$-7.03$	$10$
Single $\pi$ Production, $C_{A5}$			$9.33$	$10$
Single $\pi$ Production, BKG			$-3.03$	$10$
CCQE cross section777Difference from the Nieves Nieves:2004wx model is set to 1.0			$7.37$	$10$
CCQE $\overline{\nu }/\nu$ ratio777Difference from the Nieves Nieves:2004wx model is set to 1.0			$4.81$	$10$
CCQE $\mu /e$ ratio777Difference from the Nieves Nieves:2004wx model is set to 1.0			$-2.74$	$10$
DIS cross section			$5.72$	$10$
DIS model comparisons888Difference from CKMT CKMT94 parametrization is set to 1.0			$4.84$	$10$
DIS $Q^2$ distribution (high W)999Difference from GRV98 Gluck:1998xa is set to 1.0			$-5.70$	$10$
DIS $Q^2$ distribution (low W)999Difference from GRV98 Gluck:1998xa is set to 1.0			$-9.21$	$10$
Coherent $\pi$ production			$-9.41$	$100$
NC/CC			$-0.10$	$20$
${\nu }_{\tau }$ cross section			$9.08$	$25$
Single $\pi$ production, ${\pi }^0/{\pi }^{\pm }$			$1.75$	$40$
Single $\pi$ production, ${\overline{\nu }}_i/{\nu }_i$ (i=$e,\mu$)101010Difference from the HernandezHernandez07 model is set to 1.0			$2.30$	$10$
NC fraction from hadron simulation			$-5.07$	$10$
${\pi }^{+}$ decay uncertainty Sub-GeV 1-ring
	$e$-like 0-decay		$-0.050$	$0.6$
	$\mu$-like 0-decay		$-0.066$	$0.8$
	$e$-like 1-decay		$-0.34$	$4.1$
	$\mu$-like 1-decay		$-0.025$	$0.3$
	$\mu$-like 2-decay		$-0.47$	$5.7$
Final state and seconday interactions111111Ref. Pinzon:2017fsi
	1-ring		$-6.83$	$10$
	Multi-ring		$6.93$	$10$
Meson exchange current121212Difference from NEUT without model from Nieves:2004wx is set to 1.0			$1.08$	$10$
$\Delta m_{21}^2$ Olive:2016xmw			$0.0$	$0.0002$
${\sin }^2\left({\theta }_{12}\right)$ Olive:2016xmw			$0.05$	$1.3$
${\sin }^2\left({\theta }_{13}\right)$ Olive:2016xmw			$0.00$	$0.11$
FC reduction			$0.27$	$1.3$
PC reduction			$-0.26$	$1$
FC/PC separation			$0.000$	$0.02$
PC stopping/through-going separation (bottom)			$-1.75$	$6.8$
PC stopping/through-going separation (barrel)			$1.13$	$8.5$
PC stopping/through-going separation (top)			$-29.0$	$40$
Non-$\nu$ background	Sub-GeV $\mu$-like		$-0.001$	$0.02$
	Multi-GeV $\mu$-like		$-0.002$	$0.07$
	Sub-GeV 1-ring 0-decay $\mu$-like		$-0.001$	$0.02$
	PC		$-0.015$	$0.49$
	Sub-GeV $e$-like (flasher event)		$-0.000$	$0.03$
	Multi-GeV $e$-like (flasher event)		$-0.001$	$0.07$
Fiducial Volume			$-0.40$	$2.00$
Ring separation		$e$-like	$-0.44$	$-0.86$
		$\mu$-like	$0.43$	$0.84$
	$>400\mathrm{MeV}$	$e$-like	$0.19$	$0.37$
		$\mu$-like	$1.03$	$2$
	Multi-GeV	$e$-like	$1.32$	$2.57$
		$\mu$-like	$-0.74$	$-1.44$
	Multi-ring Sub-GeV	$e$-like	$-0.42$	$-0.82$
		$\mu$-like	$0.42$	$0.82$
	Multi-ring Multi-GeV	$e$-like	$-0.66$	$-1.28$
		$\mu$-like	$0.59$	$1.14$
Particle identification (1 ring)
	Sub-GeV	$e$-like	$-0.11$	$0.52$
		$\mu$-like	$0.10$	$-0.49$
	Multi-GeV	$e$-like	$-0.02$	$0.09$
		$\mu$-like	$0.02$	$-0.09$
Particle identification (multi-ring)
	Sub-GeV	$e$-like	$0.29$	$-1.33$
		$\mu$-like	$-0.12$	$0.53$
	Multi-GeV	$e$-like	$-0.27$	$1.25$
		$\mu$-like	$0.16$	$-0.73$
Multi-ring likelihood selection	Multi-ring $e$-like	${\nu }_e$,${\overline{\nu }}_e$	$-0.42$	$-0.71$
	Multi-ring Other		$0.32$	$0.52$
Energy calibration			$-0.24$	$2.02$
Up/down asymmetry energy calibration			$-0.53$	$0.67$
UP-$\mu$ reduction	Stopping		$0.048$	$0.5$
	Through-going		$0.029$	$0.3$
UP-$\mu$ stopping/through-going separation			$0.018$	$0.6$
Energy cut for stopping UP-$\mu$			$0.30$	$1.7$
Path length cut for through-going UP-$\mu$			$-0.33$	$1.5$
Through-going UP-$\mu$ showering separation			$3.26$	$3$
Background subtraction for UP-$\mu$
	Stopping131313The uncertainties in BG subtraction for upward-going muons are only for the most horizontal bin, .		$-4.9$	$11$
	Non-showering131313The uncertainties in BG subtraction for upward-going muons are only for the most horizontal bin, .		$2.3$	$17$
	Showering131313The uncertainties in BG subtraction for upward-going muons are only for the most horizontal bin, .		$-6.55$	$24$
${\nu }_e/{\overline{\nu }}_e$ Separation			$1.89$	$2.57$
Sub-GeV 2-ring ${\pi }^0$			$-0.039$	$1.03$
Decay-e tagging			$0.29$	$0.7$

Equations

$$P({\nu }_e\to {\nu }_e)$$

$$P({\nu }_{\mu }\to {\nu }_{\mu })$$

$$P({\nu }_{\mu }\leftrightarrow {\nu }_e)$$

$${\sin }^{2}2{\theta }_{13,M}=\frac{{\sin }^{2}2{\theta }_{13}}{(\cos 2{\theta }_{13}-A_{CC}/\Delta m_{32}^2{)}^2+{\sin }^{2}2{\theta }_{13}}$$

$$\Delta m_{32,M}^2=\Delta m_{32}^2\sqrt{(\cos 2{\theta }_{13}-A_{CC}/\Delta m_{32}^2{)}^2+{\sin }^{2}2{\theta }_{13}}$$

$$L(\Gamma ,\theta )$$

$$t_i^{\text{exp}}=t+s_{\text{mid}}/c+|{\mathbf{R}}_i^{\text{PMT}}-\mathbf{x}-s_{\text{mid}}\mathbf{d}|/c_n,$$

$$G(\mathbf{x},t)$$

$$T_i^{\text{res}}$$

$$F(t):=0.25\underset{i\in M}{\mathrm{argmax}}\{\frac{G(\mathbf{x},t_i)}{(1+((t-t_i)/\gamma {)}^2}\}+\eta$$

$${\chi }^2=2\sum _n\left(E_n-{\mathcal{O}}_n+{\mathcal{O}}_n\ln \frac{{\mathcal{O}}_n}{E_n}\right)+\sum _i{\left(\frac{{ϵ}_i}{{\sigma }_i}\right)}^2,$$

$$E_n$$

$$\text{mbox}{CL}_s^{\text{mbox}O}=\frac{p_0(O_2)}{1-p_0(O_1)},$$

$$\text{mbox}{CL}_s^{\text{mbox}H}=\frac{p_0(\text{mbox}IH)}{1-p_0(\text{mbox}NH)},$$

Full text

1 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 Pontecorvo:1967fh ; Maki:1962mu . Though most of its parameters have been experimentally measuredOlive:2016xmw , the ordering of the mass states with the largest splitting (known as the mass hierarchy), the octant of the atmospheric mixing angle $\theta_{23}$, and the value of its CP-violating phase are unknown. These unresolved issues have been at the forefront of results from the T2K Abe:2017uxa , NOvA Adamson:2017gxd and Super-Kamiokande Abe:2017osc (Super-K, SK) experiments and are the focus of next-generation experiments planned in the U.S. Acciarri:2016crz , China An:2015jdp , and Japan Abe:2018uyc .

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 value of $\theta_{13}$, $\theta_{23}$, and $\delta_{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 Abe:2017osc 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 kiloton$\cdot$year exposure with a 32% larger fiducial volume (FV) that previous Super-K analyses. Section 2 briefly reviews atmospheric neutrino oscillations before the detector is introduced in Section 3. In Section 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 Section 5. Based on this performance improvement the expansion of the FV and corresponding systematic errors are presented in Sections 6 and  7, respectively. An analysis of the atmospheric neutrino data by themselves is presented in Section 8.1 and is followed by an analysis in which $\theta_{13}$ is constrained by the world data in Section 8.2, before concluding in Section 9.

2 Neutrino Oscillations

Neutrino oscillations are a result of the neutrino mass eigenstates differing from their weak-interaction (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 have non-zero values based on the results of reactor, atmospheric, solar, and long-baseline neutrino experiments Olive:2016xmw . For brevity the details of the PMNS oscillation formalism will be omitted here, more details about the treatment of neutrino oscillation in Super-K are at Abe:2017osc .

In vacuum the leading terms in the oscillations of atmospheric $\nu_{e}$ and $\nu_{\mu}$ in can be written as:

[TABLE]

using the assumption that $\Delta m_{31}^{2}\approx\Delta m_{32}^{2}$ since $\Delta m_{21}^{2}$ 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 $\Delta m^{2}$ is eV2. Atmospheric neutrinos are sensitive to $\sin^{2}2\theta_{13}$ and the mass splitting $\Delta m_{32}^{2}$ via $\nu_{\mu}\leftrightarrow\nu_{e}$ oscillations, which depend on $\sin^{2}\theta_{23}$.

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 $\nu_{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 Wolfenstein:1978msw ; Mikheev:1985msw . When neutrinos travel through homogeneous matter, the oscillation parameters $\sin^{2}\theta_{13}$ and $\Delta m_{32}^{2}$ in Equation 1 are replaced by their matter-equivalents Giunti:1997MSW ,

[TABLE]

[TABLE]

where $A_{CC}=\pm 2\sqrt{2}G_{F}N_{e}E_{\nu}$ and the electron density $N_{e}$ is assumed to be constant. Here $G_{F}$ represents the Fermi constant.

Resonant enhancement of these matter variables occurs when $A_{CC}/\Delta m^{2}_{32}=\cos 2\theta_{13}$ and is dependent upon the sign of $\Delta m_{32}^{2}$ 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 ($\Delta m_{32}^{2}>0$) whereas it only happens for antineutrinos if it is inverted ($\Delta m_{32}^{2}<0$). Though atmospheric neutrinos generally experience a varying matter profile, and hence electron density $N_{e}$ changes as they travel through the earth, this phenomena is nonetheless present. 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) (c.f. Abe:2017osc ). Typically neutrinos with a few to ten GeV of energy that travel through the Earth’s core experience matter effects the most strongly.

3 The Super-Kamiokande Detector

Super-Kamiokande is a cylindrical water Cherenkov detector with a total volume of 50 kilotons and 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 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 1,885 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 May 31${}^{\text{st}}$ 2018 when the detector began upgrade work for its next phase. Though the topology of the detector in the SK-III and -IV periods are the same, at the start of SK-IV the front end electronics were upgraded to a system based on ASIC that uses a high-speed charge-to-time converter Nishino:2009zu . 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, its electronics are presented in Fukuda:2002uc ; Nishino:2009zu ; Abe:2013gga .

Atmospheric neutrino interactions in the detector are simulated using the Honda et. al flux calculation Honda:2011nf and the NEUT Hayato:2002sd neutrino interaction software (version 5.4.0). Particles emerging from the interactions are tracked through a simulation of the detector based on GEANT3 Brun:1994aa . The present work uses an updated version of NEUT relative to the previous atmospheric neutrino analysis (c.f. Abe:2017osc ). In particular, charged-current quasi-elastic (CCQE) interactions in the new simulation are on the local Fermi-Gas model of Nieves Nieves:2013wx ; Nieves:2011wx , assuming an axial mass $M_{A}=1.05~{}\mbox{GeV/c}^{2}$ and using the random phase approximation correction. Further, deep inelastic scattering (DIS) is modeled using the GRV98 parton distribution function Gluck:1998xa 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 Bodek:2005de , where $q^{2}$ denotes the square of the transferred four-momentum.

4 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 upward-going muons (Up-$\mu$).

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-$\mu$ 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, whose detail can be found in Ashie:2005osc .

As an example, five steps of data reduction criteria are used for FC sample. In the first and second reduction steps, criteria on the total charge collected by ID and the number of hits on OD 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 Section 6 in detail. An event’s visible energy ($E_{\mbox{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 to the reduction process.

After reduction, the FC data are sub-divided 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 ($\mu$-like), and sub-GeV ($E_{\mbox{vis}}<1330.0$ MeV) or multi-GeV ($E_{\mbox{vis}}>1330.0$ MeV) categories, as well as the number of decay electrons (0, 1, or 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-$\mu$ 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 pre-existing algorithm and are divided into “stopping” and “through-going” subsamples based on the estimated muon stopping point for PC and Up-$\mu$ events. The “through-going” events in the Up-$\mu$ sample are further divided in to “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 T2K experiment Abe:2013nka and detects accelerator-generated neutrinos, they are excluded based on event time information in this analysis. After all selections there are a total of 13 FC analysis samples, two PC samples and three Up-$\mu$ samples.

5 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 Shiozawa:1999ap . The interaction vertex is reconstructed based on hit timing information after accounting for the photon time of flight. Then the direction of the first found ring, 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 predominately produce crisp edges. In the last step of the fit, the momentum of each ring is evaluated based on the observed charge inside cone with a half angle of $70^{\circ}$ 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 experimentPat:2009mini , 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 time scales as APFit used to be. FiTQun has already been used in the T2K analyses using an expanded fiducial volume due its improved resolution of reconstructed quantities and particle identification Abe:2018wpn .

5.1 Likelihood Function

An event topology hypothesis $\Gamma$ (e.g. single-ring $e$-like) together with its associated kinematic parameters $\theta$, 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 the hypothesis to estimate the kinematic variables:

[TABLE]

In this equation, the index $j$ runs over all PMTs which 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 ($\Gamma,\theta$) is calculated as $P_{j}(\text{unhit}|\Gamma,\theta)$. PMTs which 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 charge likelihood $f_{q}(q_{i}|\Gamma,\theta)$. The likelihood density of producing a hit at the observed time $t_{i}$ is defined similarly as $f_{t}(t_{i}|\Gamma,\theta)$.

Since the processes of particle and optical photon propagation are decoupled from the response of the PMT and the electronics, the charge likelihood can be rewritten in terms of the expected number of photoelectrons produced at the $i$-th PMT given the hypothesis (the predicted charge), $\mu_{i}(\Gamma,\theta)$, as

[TABLE]

$P_{j}(\text{unhit}|\mu_{j})$ and $f_{q}(q_{i}|\mu_{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 $\mu_{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 $\mu_{i}$. The final charge likelihood $f_{q}(q_{i}|\mu_{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^{\text{exp}}_{i},\Gamma,p,\mu_{i})$, where $p$ is the momentum for the topology $\Gamma$ and $t^{\text{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

[TABLE]

where $\bm{x}$ and $t$ are the vertex and creation time of the particle, respectively, $\bm{d}$ is the particle direction, $\bm{R}_{i}^{\text{PMT}}$ is the position of the $i$-th PMT and $s_{\text{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 a narrower distribution of hit time for higher numbers of incident photons and hence, more predicted charge. The track length of a particle $s$, which is determined by the topology $\Gamma$ and momentum $p$, also affects the shape of 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 the photons from a particle with earlier $t^{\text{exp}}_{i}$ always arrive earlier than the photons from any other particles with later $t^{\text{exp}}_{i}$ values.

Once the likelihood function is defined, the best set of kinematic parameters $\hat{\theta}$ for a given event topology hypothesis is defined as that which maximizes $L(\Gamma,\theta)$. The best estimate for the particle content of a given event is determined by comparison of $L(\Gamma,\hat{\theta})$ among all hypotheses, $\Gamma$.

5.2 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 candidate for particle activity. 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(\Gamma,\theta)$ is minimized with respect to $\theta$, using the MINUITJames:1994vla package.

5.2.1 Vertex pre-fitting

The vertex pre-fitter is a fast algorithm which 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 $\bm{x}$ and time $t$ which maximize the goodness function

[TABLE]

where

[TABLE]

is the residual hit time calculated on the assumption of a point-like light source and subtracting the photon time of flight. The value of parameter $\sigma$ determines the precision of pre-fitting. Here $\bm{R}_{i}^{\text{PMT}}$ is the position of the $i^{th}$ PMT. When the vertex and time get close to their true values by doing a grid search, the $T_{i}^{\text{res}}$ distribute near zero, which results in a large value of the goodness. The pre-fitter is executed several times with gradually shrinking the grid size and $\sigma$ to achieve high precision efficiently. The fitted vertex from this step, called pre-fit vertex, is just a rough estimation and will be fitted again with higher precision during minimization of $L(\Gamma,\theta)$.

5.2.2 Hit clustering

Events in Super-K are defined by detector activity in an $O(10~{}\mu$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(\bm{x},t)$ from Equation 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 Figure 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:

[TABLE]

where $M$ represents all local maxima of goodness function $G(\bm{x},t)$. The value of time constant $\gamma$ is 25 ns when $t<t_{i}$ and 70 ns otherwise. An offset $\eta$ = 9 is added to the threshold function to suppress the effect of dark hits. The function $F(t)$ is the blue curve in Figure 1. The minimum goodness between any two peaks must be lower than a second threshold, $0.6\times F(t)$, which is shown as green dashed curve. Under these criteria only the peaks labelled 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 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 $\bm{x}$ is then close to the vertex of the secondary particle.

A time window defined as -180 ns $<T_{i}^{\text{res}}<$ 800 ns around each found peak 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.

5.2.3 Single-ring reconstruction

The most basic reconstruction which 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 Equation 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 example, are separated by cutting on ln($L_{e}/L_{\mu}$), the logarithm of the likelihood ratio between the best-fit electron and muon hypotheses. The distribution of this variable is shown for FC atmospheric neutrino data and Monte Carlo (MC) simulation result with sub-GeV and multi-GeV energies in Figure 2. In both cases there is a clear separation of the likelihood variable between electron-like ($e$-like) and muon-like ($\mu$-like) events.

Figure 3 shows the vertex resolution as the function of visible energy for the FC single-ring charged-current 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 $\nu_{e}$ events is stable at 20.64 cm , while the resolution of APFit varies from 34.63 cm to 25.26 cm. For the CCQE $\nu_{\mu}$ events, the vertex resolution of fiTQun, which changes from 29.24 cm to 15.89 cm, is better than that of APFit, which varies from 34.21 cm to 18.16 cm. In the same energy range, the momentum resolution of fiTQun for CCQE $\nu_{e}$ events improves from 5.39% to 2.58% as the visible energy increases, while the resolution of APFit changes from 7.04% to 3.32%, as shown in the top left plot of Figure 4. For CCQE numu 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% mis-identification rate for visible energies less than 1330 MeV (Figure 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 or better than APFit.

5.2.4 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 Equation 5 is updated to include a new ring and minimized again allowing the kinematic parameters for the rings to vary. Three hypotheses, including $e$-like, $\mu$-like and $\pi^{+}$-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 2-ring and a 1-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 event when the first ring hypothesis is $e$-like ($\mu$-like). However, ring candidates that have an angular separation of less than $20^{\circ}$ from the most energetic ring are discarded as spurious since MC studies indicate such candidates are typically due to the particle scattering rather than a new particle. When this occurs the two rings are merged and refit 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 energy more than 30 MeV above 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 that of 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 $\mu$-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 $\pi^{0}$s in charged-current single pion (CC1$\pi^{0}$) events. Figure 8 shows the reconstructed invariant mass calculated using the second and the third most energetic rings for events with three $e$-like rings ($eee$ events) for both reconstruction algorithms. Shaded histograms show true CC$\nu_{e}1\pi^{0}$ events. Both fitters have a peak in the invariant mass distribution near the $\pi^{0}$ mass, but the fiTQun peak is larger, sharper, and has fewer backgrounds indicating improved reconstruction efficiency over APFit. The FWHM of the $\pi^{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$\nu_{e}1\pi^{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 those values of APFit: 175.9 and $42.9\%$. Similarly, CC$\nu_{\mu}1\pi^{0}$ events from the three ring sample whose leading ring is $\mu$-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.

Another important performance indicator for multi-ring event reconstruction algorithm is the separation of 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 $\nu_{e}+\bar{\nu_{e}}$ interactions based on a likelihood selection as in Wendell:2010md using the APFit algorithm. 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 Figure 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 Abe:2017osc and uses the number of decay electrons, the number of reconstructed rings, electrons, 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 $\bar{\nu}_{e}$ ($\nu_{e}$) events as $\bar{\nu}_{e}$-like ($\nu_{e}$-like) is summarized in Table 4. FiTQun shows improved performance in each of these metrics which will translate to better sensitivity to sin${}^{2}\theta_{13}$ and the mass hierarchy.

6 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 is similar to the one used in Abe:2017osc , but all cut variables used to define the FC samples are reconstructed by fiTQun. Partially contained and Up-$\mu$ events on the other hand, are defined using APFit in the same manner as in the previous publications.

6.1 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 Section 4, events which 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 analysis sample. Such “entering events” are a type of mis-reconstructed background in the oscillation analysis with potentially large systematic errors.

Figure 11 shows the purity of several components of the FC sub-GeV $\mu$-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, $\nu_{\mu}$CC and $\bar{\nu}_{\mu}$ 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 $\nu$ events increases dramatically. All of 13 FC event categories exhibit similar behavior, indicating that a FV cut at 50 cm is an acceptable limit due to entering $\nu$ events. The purity of the 13 FC samples defined using this cut and the conventional FV definition are summarized in Table 5 (c.f. Table II of Ref Abe:2017osc ).

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 a 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.

6.2 Background estimation

In addition to entering backgrounds and events with incorrectly assigned PID, non-neutrino (non-$\nu$) 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, that is the entrance point of the muon, is expected to be on the wall of 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 Figure 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.

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 observed in the FC sample with dwall $>$ 50 cm. Most of the cosmic ray muon backgrounds were reconstructed as multi-ring downward-going $\mu$-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 mis-reconstructions is consistent with that found in independent 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 eye-scanned background, this assumption has been validated for all particle directions and momenta. As an example, Figure 14 shows the dwall distribution for multi-ring $\mu$-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 shape 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.

Expanding the FV opens the possibility that events interacting within the ID but close to its wall may produce particles which 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 Figure 16. Their event rates since the start of SK-IV of the experiment 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-$\mu$ events per day, as shown in Figure 17.

7 Systematic error

In order to use the fiTQun reconstruction algorithm for physics analysis the corresponding systematic errors must be evaluated. While systematic errors related to the atmospheric neutrino flux and cross section model are the same as those used in Abe:2017osc , 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-$\mu$ 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 control samples spanning momenta up to 10 GeV/c: electrons from cosmic ray muon decay, atmospheric neutrino events producing single $\pi^{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 Figure 18. The total systematic error follows that of Abe:2017osc and is the sum in quadrature of the absolute energy scale, which is defined as the largest data-MC difference across all samples, and momenta with the average time variation of these samples throughout SK-IV. The uncertainty from the up/down asymmetry of the detector is 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 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% Abe:2017osc .

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.

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 next the 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 next section. Systematic errors and their sizes at the best fit point of the analysis are presented in Tables LABEL:tab:sysa.

8 Atmospheric Neutrino Oscillation Analysis

The Super-K atmospheric neutrino data are separately fit against both the normal and inverted hierarchy hypotheses. Since $\theta_{13}$ has been measured and is well constrained by reactor experiments, two types of analyses are done: one with $\theta_{13}$ as a free fitting parameter and the other with it constrained to $0.0210\pm 0.0011$ Olive:2016xmw . 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 $\theta_{13}$ and therefore measurements with atmospheric neutrinos represent an important test of the analysis’ hierarchy preference. The solar mixing parameters, $\Delta m^{2}_{12}$ and $\sin^{2}\theta_{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.

Data are fit to the MC expectation with a binned $\chi^{2}$ method assuming Poisson statistics. The effect of systematic errors are regarded as scaling factors on the MC in each bin Fogli:2002pt , with different error sources assumed to be independent. The definition of $\chi^{2}$ used in the fits is

[TABLE]

where,

[TABLE]

In this equation $\mathcal{O}_{n}$ is the observed number of events in the $n^{th}$ 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{\sigma_{i}}$ variation of the $i^{th}$ systematic error source. The sum over $(\epsilon_{i}/\sigma_{i})^{2}$ penalizes the $\chi^{2}$ for adjusting the systematic errors when bringing the MC into agreement with the data. For each set of oscillation parameters tested, the MC fit to the data by minimizing $\chi^{2}$ over the error parameters $\epsilon_{i}$. All fits are performed over 515 analysis bins in SK-IV.

8.1 Analysis with unconstrained $\theta_{13}$

The fit with $\theta_{13}$ unconstrained by reactor measurements is performed over four parameters, $|\Delta m^{2}_{32,31}|$, $\sin^{2}\theta_{13}$,$\sin^{2}\theta_{23}$ and $\delta_{CP}$ with a total of 71 systematic errors. Fits to the normal hierarchy use $\Delta m^{2}_{32}$ while those for the inverted hierarchy use $\Delta m^{2}_{31}$. The agreement between the data and MC is evaluated using Equation 10 at each point in the grid of points shown in Table 8, while keeping the solar mixing parameters, $\Delta m^{2}_{12}$ and $\sin^{2}\theta_{12}$, fixed to the values in Table 9. Systematic errors representing the uncertainty in these parameters are included in the analysis. The best-fit parameters for each hierarchy are defined as the parameter set returning the smallest $\chi^{2}$ value. Between the two hierarchy fits the one with smallest minimum $\chi^{2}$ value is taken as the global best-fit and is taken to indicate the mass hierarchy preference.

Results and Discussion

Figure 19 shows the expected sensitivity to reject a zero value of $\sin^{2}\theta_{13}$ as a function of the true value of sin${}^{2}\theta_{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 number of atmospheric neutrino events selected by APFit and fiTQun are 24188.7 and 24184.3, respectively, while the number of the events passed both the cut of the two algorithms is 23520.2. The 3% difference is mainly due to the vertex resolution near FV boundary.

However, since the contamination of $\nu_{\mu}$ events in the e-like samples most sensitive to $\theta_{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 $\theta_{13}$, $|\Delta m^{2}_{32,31}|$, $\sin^{2}\theta_{23}$ and $\delta_{CP}$ from the fit to the data are shown in Figure 20. The normal hierarchy hypothesis yields better data-MC agreement than the inverted hierarchy hypothesis with $\chi^{2}_{\text{NH,min}}-\chi^{2}_{\text{IH,min}}=-1.81$. The 1$\sigma$ allowed region for $\sin^{2}\theta_{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}\theta_{13}=0.0$ is disfavored at approximately $1.8{\sigma}$ (1.2 $\sigma$) for normal (inverted) fit. A summary of the best-fit information and parameter constraints is presented in Table 11.

As discussed in Section 2, the determination of $\sin^{2}\theta_{13}$ and the mass hierarchy using atmospheric neutrinos is achieved using 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 downward-going and downward-going events relative to their sum. Here downward-going (upward-going) events are events whose zenith angle satisfy $\mbox{cos}\theta_{z}<-0.4$ ($\mbox{cos}\theta_{z}>0.4$). Excesses between a few and ten GeV in the multi-GeV e-like $\nu_{e}$ and the Multi-Ring e-like $\nu_{e}$ samples, drive the normal hierarchy preference.

The best-fit atmospheric mixing parameters from the normal hierarchy fit are $\Delta m^{2}_{32}=2.63^{+0.10}_{-0.21}\times 10^{-3}\mbox{eV}^{2}$ and $\sin^{2}\theta_{23}=0.588^{+0.030}_{-0.062}$ for second octant (best-fit) and $0.425^{+0.051}_{-0.034}$ for first octant. Maximal mixing ($\sin^{2}\theta_{23}=0.5$) is weakly disfavored at around $1{\sigma}$ significance. Besides the data excesses in the upward-going regions of the single-ring e-like $\nu_{e}$ sample and multi-ring other sample, the data deficits in the multi-GeV $\mu$-like sample also contribute to this preference. It should be noted that the preference for $\sin^{2}\theta_{23}$ is coupled to the $\theta_{13}$ measurement, since both parameters feature in the $\nu_{\mu}\rightarrow\nu_{e}$ oscillation probability.

The best-fit value of $\delta_{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 appearance in the sub-GeV e-like samples are expected to be observed due to $\nu_{\mu}\rightarrow\nu_{e}$ oscillations. In the multi-GeV region, $\delta_{CP}$ similarly modulates the amount of electron neutrino appearance but it is subdominant to the effects of $\theta_{13}$. In the inverted hierarchy fit, the value of $\theta_{13}$ can be more easily adjusted to bring the MC prediction into better agreement with data since only antineutrino events experience matter effects in the earth. Since there are more neutrino events than antineutrino events in the atmospheric sample, the freedom to move $\theta_{13}$ ultimately leads to a weaker constraint on $\delta_{CP}$.

8.2 Analysis with constrained $\theta_{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 and $\delta_{CP}$. Indeed, Equation 2 shows that the size of MSW resonance relies on the value of $\mbox{sin}^{2}\theta_{13}$ directly. Since reactor neutrino experiments constrain this parameter more precisely than the analysis in the previous section, the analysis presented here constrains the value of $\sin^{2}\theta_{13}$ to $0.0210\pm 0.0011$ Olive:2016xmw . 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 sensitivity to the mass hierarchy and $\theta_{23}$ octant for true normal hierarchy are illustrated in Figure 22 and Figure 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 $\theta_{13}$ in the previous section.

Figure 24 shows the $\chi^{2}$ value as a function of the atmospheric neutrino mixing parameters and $\delta_{CP}$ in the $\theta_{13}$-constrained fit. As in the unconstrained fit the data prefer the normal mass hierarchy with $\Delta\chi^{2}=-2.45$. The 1$\sigma$ allowed region for $|\Delta m^{2}_{32}|$ is from 2.41 to 2.75 $\times 10^{-3}\mbox{eV}^{2}$ (2.36 to 2.67 $\times 10^{-3}\mbox{eV}^{2}$) for the normal (inverted) hierarchy fit, which is consistent with the result of unconstrained fit. The contour for the allowed region with 90% CL is shown in Figure 25. The preference for the second octant of $\theta_{23}$ in the unconstrained fit is changed. Adding the constraint on $\sin^{2}\theta_{13}$ the preferred value of $\theta_{23}$ shifts from the second to first octant, though both are allowed at 1$\sigma$. Indeed the difference in $\chi^{2}$ between the minimum value of each octant is $\Delta\chi^{2}_{1^{st}-2^{nd}}$ = -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\theta_{13}$ and $\sin^{2}\theta_{23}$ as shown in Equation 1; an increase in one can be compensated by a decrease in the other to describe the data. The same upward-going event excesses (deficits) in the $e$-like ($\mu$-like) samples that drive the hierarchy, $\theta_{13}$, and octant preferences in the unconstrained analysis are not strong enough to support the increased excess (deficits) expected for higher values of $\theta_{13}$ and the second octant. Since the $\theta_{13}$ is now fixed, the data can be accommodated by reducing the expected number of $e$-like events and increasing the number of $\mu$-like events by moving $\theta_{23}$ to the first octant.

The best-fit value of $\delta_{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{\sigma}$ errors are summarized in Table 11.

The $\mbox{CL}_{s}$ method Read:2002hq is used to address the significance of the octant and mass hierarchy preferences observed in the data. As in the previous Super-K study Abe:2017osc 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 the current detector exposure and systematic errors have been varied as Gaussian parameters with widths specified by their uncertainties. The $\mbox{CL}_{s}$ parameter for the octant study, $\mbox{CL}^{\mbox{O}}_{s}$, is defined as:

[TABLE]

where $p_{0}(O_{i})$ represents the $p$-value with respect to obtaining a difference in the best-fit $\chi^{2}$ values between both octant hypotheses smaller (larger) than the value of data assuming the true value of $\mbox{sin}^{2}\theta_{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 Figure 26, respectively. The $\mbox{CL}^{\mbox{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

[TABLE]

where $p_{0}(\mbox{IH})$ and $p_{0}(\mbox{NH})$ are p-values for obtaining a difference in the $\chi^{2}$ of the best-fit mass hierarchies more extreme than that of the data assuming a true IH and NH, respectively.

Figure 27 shows the distribution for the mass hierarchy determination. Due to the large uncertainty on $\theta_{23}$, MC ensembles have been generated with different assumed values $\theta_{23}$ and with the other oscillation parameters fixed to their best-fit values. The $p$-values and $\mbox{CL}^{\mbox{H}}_{s}$ values for hierarchy test are summarized in Table 12. As sin${}^{2}\theta_{23}$ ranges from 0.4 to 0.6 the observed $\mbox{CL}^{\mbox{H}}_{s}$ values decrease from a 30.8% C.L. preference for the inverted mass hierarchy to 14.3%. At the analysis best-fit value the preference is 26.0%. A smaller value of $\mbox{sin}^{2}\theta_{23}$ predicts less electron neutrino appearance, which results in a smaller $p_{0}(\mbox{IH})$ values and larger $\mbox{CL}^{\mbox{H}}_{s}$ due to the multi-GeV $e$-like samples’ event excesses seen in data.

9 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 the normal mass hierarchy constraints on the atmospheric mixing parameters are $\mbox{sin}^{2}\theta_{23}=0.425^{+0.046}_{-0.037}$ ($0.600^{+0.013}_{-0.030}$) for first (second) octant and $\Delta m^{2}_{32}=2.53^{+0.22}_{-0.12}$, with $\delta_{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 $\theta_{23}$ octant is observed.

Figure 28 shows the expected sensitivity to mass hierarchy as the 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 is 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 2,800 days of livetime.

10 Acknowledgments

We gratefully acknowledge the cooperation of the Kamioka Mining and Smelting Company. The Super-Kamiokande experiment has been built and operated from funding by the Japanese Ministry of Education, Culture, Sports, Science and Technology, the U.S. Department of Energy, and the U.S. National Science Foundation. Some of us have been supported by funds from the National Research Foundation of Korea NRF-2009-0083526 (KNRC) funded by the Ministry of Science, ICT, and Future Planning, the European Union H2020 RISE-GA641540-SKPLUS, the Japan Society for the Promotion of Science, the National Natural Science Foundation of China under Grants No. 11235006, the National Science and Engineering Research Council (NSERC) of Canada, the Scinet and Westgrid consortia of Compute Canada, and the National Science Centre of Poland (2015/17/N/ST2/04064, 2015/18/E/ST2/00758).

11 Systematic Uncertainties