PHYSICAL REVIEW D 97, 072001 (2018)

Atmospheric neutrino oscillation analysis with external constraints
in Super-Kamiokande I-IV
K. Abe,1,38 C. Bronner,1 Y. Haga,1 Y. Hayato,1,38 M. Ikeda,1 K. Iyogi,1 J. Kameda,1,38 Y. Kato,1 Y. Kishimoto,1,38 Ll. Marti,1
M. Miura,1,38 S. Moriyama,1,38 M. Nakahata,1,38 T. Nakajima,1 Y. Nakano,1 S. Nakayama,1,38 Y. Okajima,1 A. Orii,1
G. Pronost,1 H. Sekiya,1,38 M. Shiozawa,1,38 Y. Sonoda,1 A. Takeda,1,38 A. Takenaka,1 H. Tanaka,1 S. Tasaka,1
T. Tomura,1,38 R. Akutsu,2 T. Irvine,2 T. Kajita,2,38 I. Kametani,2 K. Kaneyuki,2,38,* Y. Nishimura,2 K. Okumura,2,38
E. Richard,2 K. M. Tsui,2 L. Labarga,3 P. Fernandez,3 F. d. M. Blaszczyk,4 J. Gustafson,4 C. Kachulis,4 E. Kearns,4,38
J. L. Raaf,4 J. L. Stone,4,38 L. R. Sulak,4 S. Berkman,5 S. Tobayama,5 M. Goldhaber,6,* G. Carminati,7 M. Elnimr,7W. R. Kropp,7 S. Mine,7 S. Locke,7 A. Renshaw,7 M. B. Smy,7,38 H. W. Sobel,7,38 V. Takhistov,7,† P. Weatherly,7
K. S. Ganezer,8 B. L. Hartfiel,8 J. Hill,8 N. Hong,9 J. Y. Kim,9 I. T. Lim,9 R. G. Park,9 T. Akiri,10 A. Himmel,10 Z. Li,10
E. O’Sullivan,10 K. Scholberg,10,38 C. W. Walter,10,38 T. Wongjirad,10 T. Ishizuka,11 T. Nakamura,12 J. S. Jang,13 K. Choi,14J. G. Learned,14 S. Matsuno,14 S. N. Smith,14 J. Amey,15 R. P. Litchfield,15 W. Y. Ma,15 Y. Uchida,15 M. O. Wascko,15
S. Cao,16 M. Friend,16 T. Hasegawa,16 T. Ishida,16 T. Ishii,16 T. Kobayashi,16 T. Nakadaira,16 K. Nakamura,16,38 Y. Oyama,16
K. Sakashita,16 T. Sekiguchi,16 T. Tsukamoto,16 KE. Abe,17 M. Hasegawa,17 A. T. Suzuki,17 Y. Takeuchi,17,38 T. Yano,17
T. Hayashino,18 S. Hirota,18 K. Huang,18 K. Ieki,18 M. Jiang,18 T. Kikawa,18 KE. Nakamura,18 T. Nakaya,18,38 N. D. Patel,18
K. Suzuki,18 S. Takahashi,18,38 R. A. Wendell,18,38 L. H. V. Anthony,19 N. McCauley,19 A. Pritchard,19 Y. Fukuda,20
Y. Itow,21,22 G. Mitsuka,21 M. Murase,21 F. Muto,21 T. Suzuki,21 P. Mijakowski,23 K. Frankiewicz,23 J. Hignight,24
J. Imber,24 C. K. Jung,24 X. Li,24 J. L. Palomino,24 G. Santucci,24 C. Vilela,24 M. J. Wilking,24 C. Yanagisawa,24,‡ S. Ito,25
D. Fukuda,25 H. Ishino,25 T. Kayano,25 A. Kibayashi,25 Y. Koshio,25,38 T. Mori,25 H. Nagata,25 M. Sakuda,25 C. Xu,25
Y. Kuno,26 D. Wark,27,33 F. Di Lodovico,28 B. Richards,28 R. Tacik,29,42 S. B. Kim,30 A. Cole,31 L. Thompson,31
H. Okazawa,32 Y. Choi,34 K. Ito,35 K. Nishijima,35 M. Koshiba,36 Y. Totsuka,36,* Y. Suda,37 M. Yokoyama,37,38
R. G. Calland,38 M. Hartz,38 K. Martens,38 B. Quilain,38 C. Simpson,38,27 Y. Suzuki,38 M. R. Vagins,38,7 D. Hamabe,39
M. Kuze,39 T. Yoshida,39 M. Ishitsuka,40 J. F. Martin,41 C. M. Nantais,41 P. de Perio,41 H. A. Tanaka,41 A. Konaka,42
S. Chen,43 L. Wan,43 Y. Zhang,43 R. J. Wilkes,44 and A. Minamino45(Super-Kamiokande Collaboration)
1

Kamioka Observatory, Institute for Cosmic Ray Research, University of Tokyo, Kamioka,
Gifu 506-1205, Japan
2
Research Center for Cosmic Neutrinos, Institute for Cosmic Ray Research, University of Tokyo, Kashiwa,
Chiba 277-8582, Japan3
Department of Theoretical Physics, University Autonoma Madrid, 28049 Madrid, Spain
4
Department of Physics, Boston University, Boston, Massachusetts 02215, USA
5
Department of Physics and Astronomy, University of British Columbia, Vancouver, British Columbia
V6T1Z4, Canada
6
Physics Department, Brookhaven National Laboratory, Upton, New York 11973, USA
7
Department of Physics and Astronomy, University of California, Irvine, Irvine,
California 92697-4575, USA
8
Department of Physics, California State University, Dominguez Hills, Carson, California 90747, USA
9
Department of Physics, Chonnam National University, Kwangju 500-757, Korea
10
Department of Physics, Duke University, Durham North Carolina 27708, USA
11
Junior College, Fukuoka Institute of Technology, Fukuoka, Fukuoka 811-0295, Japan
12
Department of Physics, Gifu University, Gifu, Gifu 501-1193, Japan
13
GIST College, Gwangju Institute of Science and Technology, Gwangju 500-712, Korea
14
Department of Physics and Astronomy, University of Hawaii, Honolulu, Hawaii 96822, USA
15
Department of Physics, Imperial College London, London SW7 2AZ, United Kingdom
16High Energy Accelerator Research Organization (KEK), Tsukuba, Ibaraki 305-0801, Japan17
Department of Physics, Kobe University, Kobe, Hyogo 657-8501, Japan18
Department of Physics, Kyoto University, Kyoto, Kyoto 606-8502, Japan
19
Department of Physics, University of Liverpool, Liverpool L69 7ZE, United Kingdom
20
Department of Physics, Miyagi University of Education, Sendai, Miyagi 980-0845, Japan
21
Institute for Space-Earth Enviromental Research, Nagoya University, Nagoya, Aichi 464-8602, Japan
22
Kobayashi-Maskawa Institute for the Origin of Particles and the Universe, Nagoya University,
Nagoya, Aichi 464-8602, Japan

2470-0010=2018=97(7)=072001(23)072001-1

Published by the American Physical Society

K. ABE et al.

PHYS. REV. D 97, 072001 (2018)23

National Centre For Nuclear Research, 00-681 Warsaw, Poland
Department of Physics and Astronomy, State University of New York at Stony Brook,
New York 11794-3800, USA
25
Department of Physics, Okayama University, Okayama, Okayama 700-8530, Japan
26
Department of Physics, Osaka University, Toyonaka, Osaka 560-0043, Japan
27
Department of Physics, Oxford University, Oxford OX1 3PU, United Kingdom
28
School of Physics and Astronomy, Queen Mary University of London, London E1 4NS, United Kingdom
29
Department of Physics, University of Regina, 3737 Wascana Parkway, Regina,
Saskatchewan S4SOA2, Canada
30
Department of Physics, Seoul National University, Seoul 151-742, Korea
31
Department of Physics and Astronomy, University of Sheffield, S10 2TN Sheffield, United Kingdom
32
Department of Informatics in Social Welfare, Shizuoka University of Welfare, Yaizu,
Shizuoka 425-8611, Japan
33
STFC, Rutherford Appleton Laboratory, Harwell Oxford, and Daresbury Laboratory,
Warrington OX11 0QX, United Kingdom
34
Department of Physics, Sungkyunkwan University, Suwon 440-746, Korea
35
Department of Physics, Tokai University, Hiratsuka, Kanagawa 259-1292, Japan
36
The University of Tokyo, Bunkyo, Tokyo 113-0033, Japan
37
Department of Physics, University of Tokyo, Bunkyo, Tokyo 113-0033, Japan
38
Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes
for Advanced Study, University of Tokyo, Kashiwa, Chiba 277-8583, Japan
39
Department of Physics, Tokyo Institute of Technology, Meguro, Tokyo 152-8551, Japan40
Department of Physics, Faculty of Science and Technology, Tokyo University of Science,
Noda, Chiba 278-8510, Japan
41
Department of Physics, University of Toronto, Ontario M5S 1A7, Canada
42
TRIUMF, 4004 Wesbrook Mall, Vancouver, British Columbia V6T2A3, Canada
43
Department of Engineering Physics, Tsinghua University, Beijing 100084, China
44
Department of Physics, University of Washington, Seattle, Washington 98195-1560, USA
45
Faculty of Engineering, Yokohama National University, Yokohama 240-8501, Japan
24(Received 25 October 2017; published 3 April 2018)
An analysis of atmospheric neutrino data from all four run periods of Super-Kamiokande optimized forsensitivity to the neutrino mass hierarchy is presented. Confidence intervals for Δm232, sin2 θ23 , sin2 θ13 andδCP are presented for normal neutrino mass hierarchy and inverted neutrino mass hierarchy hypotheses,
based on atmospheric neutrino data alone. Additional constraints from reactor data on θ13 and frompublished binned T2K data on muon neutrino disappearance and electron neutrino appearance are added tothe atmospheric neutrino fit to give enhanced constraints on the above parameters. Over the range ofparameters allowed at 90% confidence level, the normal mass hierarchy is favored by between 91.9% and94.5% based on the combined Super-Kamiokande plus T2K result.
DOI: 10.1103/PhysRevD.97.072001

I. INTRODUCTION
The principal goal of contemporary neutrino oscillation
experiments is to fully test the three-neutrino mixing*

Deceased.
Also at Department of Physics and Astronomy, UCLA, CA
90095-1547, USA.
‡
Also at BMCC/CUNY, Science Department, New York,USA.
†

Published by the American Physical Society under the terms of
the Creative Commons Attribution 4.0 International license.
Further distribution of this work must maintain attribution to
the author(s) and the published article’s title, journal citation,
and DOI. Funded by SCOAP3.

paradigm based on the Pontecorvo-Maki-NakagawaSakata (PMNS) matrix [1,2]. This paradigm is characterized by three mixing angles, two mass splittings, and one
CP-violating phase. Some neutrino mixing parameters
have been experimentally determined, such as the magnitude of the two mass splittings, the ordering of the massstates with the smallest splitting, and the values of themixing angles. In particular, measurements by reactor
antineutrino [3–5] experiments and T2K [6] have established that the mixing angle θ13 is small but non-zero andthey have precisely measured its value. There remain
unknown parameters in the PMNS formalism, most notably
the ordering of the mass states with the largest splitting,which is mathematically expressed as the sign of Δm231 , and072001-2

ATMOSPHERIC NEUTRINO OSCILLATION ANALYSIS …
is commonly referred to as the neutrino mass hierarchy.
Although it is known that muon and tau neutrino mixing
is nearly maximal, i.e. θ23 is near π=4, it is not known if
θ23 takes exactly that value, or is slightly larger orslightly smaller [7,8]. With all three neutrino flavors andmass states mixing, it is possible to measure the unknown
CP-violating phase δCP and perhaps conclude thatneutrinos and antineutrinos have different oscillation
probabilities, if it is found that δCP is neither 0 nor π.
The value of δCP is considered to be unknown, although
the T2K and NOvA long-baseline experiments, and theresults published in this paper, are beginning to constrain
it [8,9].
Due to the presence of neutrinos and antineutrinos, theeffects of matter on neutrino oscillations, and the wide
variety of energies and pathlengths spanned, atmospheric
neutrinos are sensitive to the unknown parameters of the
PMNS formalism. The measurement of the mass hierarchy is driven by an expected hierarchy-dependent,
upward-going excess of either electron neutrino or antineutrino interactions driven by θ13 -induced matter effects
between two and ten GeV. In order to take advantage ofthis phenomenon, sign selection of neutrino interactions
and sufficient statistics are necessary. It should be notedthat the size of this event excess is a function of θ23 , andas will be discussed below, constraints on this parameter
improve sensitivity to the hierarchy. Determining themass hierarchy and measuring θ23 play an importantrole in interpreting any neutrino versus antineutrino
oscillation difference and thereby establishing CP
violation.
In this paper we analyze 328 kiloton · years of SuperKamiokande (Super-K) atmospheric data. The sensitivityof our experiment is not sufficient to definitively resolve theunknown parameters. In particular we are limited by low
statistics and difficult event classification in the highenergy hierarchy-sensitive sample. Nevertheless, we analyze the atmospheric neutrino data in a manner optimized
for sensitivity to the mass hierarchy and report our best
estimates and confidence intervals. We present results withand without constraints from external experiments. In
Sec. II atmospheric neutrino oscillations are reviewed
before discussing the Super-K detector and data set inSec. III. An analysis of the atmospheric neutrino data bythemselves is then presented in Sec. IV and followed by ananalysis incorporating constraints from external measurements in Sec. V. These results are interpreted in Sec. VI
before concluding in Sec. VII.
II. OSCILLATIONS
A. In vacuum
Neutrinos oscillate because the neutrino eigenstates of
the weak interaction are different from the neutrino mass
eigenstates. The flavor eigenstates να are related to the mass
eigenstates νi by

PHYS. REV. D 97, 072001 (2018)jνα i ¼

3
Xi

U α;i jνi i;

ð1Þwhere U is the 3 × 3 Pontecorvo-Maki-Nakagawa-Sakata(PMNS) matrix [1,2]
0
10
1
1
0
0
0 s13 e−iδCP
c13
B
CB
C
U ¼ @ 0 c23 s23 A@
0
1
0
A
iδ
CP
0 −s23 c23
0
c13−s13 e10
c12 s12 0CBð2Þ
× @ −s12 c12 0 A:
001

Here cij ¼ cos θij ; sij ¼ sin θij . Propagation of these statesaccording to their vacuum Hamiltonians leads to thestandard oscillation formula for relativistic neutrinos in
vacuumX
Pðνα → νβ Þ ¼ δαβ − 4 ℜðU αi Uβi Uαj U βj Þsin2 Δij
i>jX 2 ℑðUαi Uβi Uαj U βj Þ sin 2Δij ;

ð3Þ

i>j

where
Δij ¼

1.27Δm2ij ðeV2 ÞLðkmÞ
EðGeVÞ

and the sign before the second summation is positive forneutrinos and negative for antineutrinos. Neutrino oscillations in vacuum are thus fully described by 6 parameters:
the 3 mixing angles θ13 , θ12 , θ23 , the two mass splittings
Δm221 ; Δm231 , and the CP-violating phase δCP . Data fromreactor, atmospheric, solar, and long-baseline neutrino
experiments indicate that nearly all of these parameters
have nonzero values [10]. Currently the sign of Δm231 andthe value of δCP are unknown. Note that throughout thispaper the indices of the mass splittings present the neutrinomass states in descending order from left to right regardlessof the hierarchy assumption.
The unoscillated atmospheric neutrino flux consists ofelectron- and muon-flavored neutrinos and antineutrinos.
Since ντ charged current interactions are either kinematically disallowed or suppressed compared to νμ and νecharged current (CC) interactions over the energy range
considered in the analysis below, the atmospheric data are
predominantly described by the νμ and νe survival probabilities and the νμ ↔ νe oscillation probability. For suffi1.27Δm2 L

ciently small L=E, sin2 ð E 12 Þ ≪ 1 and so the Δm212
terms in Eq. (3) can be ignored and the approximation
Δm231 ≈ Δm232 applied. Under these assumptions, the dominant νe and νμ oscillation probabilities become:

072001-3

K. ABE et al.

PHYS. REV. D 97, 072001 (2018)
1.27Δm231 LPðνe → νe Þ ≅ 1 − sin 2θ13 sin
E
2

TABLE I. Model of the Earth used in the analysis, a simplified
version of the PREM.

2Pðνμ → νμ Þ ≅ 1 − 4cos2 θ13 sin2 θ23 ð1 − cos2 θ13 sin2 θ23 Þ2
2 1.27Δm31 L
× sin
E


2
2
2
2 1.27Δm31 L
Pðνμ ↔ νe Þ ≅ sin θ23 sin 2θ13 sin:
ð4ÞERegionRmin (km)Rmax (km)density (g=cm3 )inner coreouter coremantle
crust

0
1220
3480
5701

1220
3480
5701
6371

13.0
11.3
5.0
3.3

X¼

k

B. In matter
When neutrinos travel through matter, the effective
Hamiltonian is modified from its vacuum form due tothe difference in the forward scattering amplitudes of νe andνμ;τ (presented here in the mass eigenstate basis):
0 m2

0

1

B 2E
B
Hmatter ¼ B 0
@
0

m22
2E

0

1

0
a
C
C
B
†
0 CþU @0
A
0
m23
0

2E

0 0

1

C
0 0 AU;

ð5Þ

0 0

pﬃﬃﬃ
where a ¼  2GF N e , Gf is the Fermi constant, N e is theelectron density, U is the PMNS matrix, and the plus
(minus) sign is for neutrinos (antineutrinos). For constantdensity matter the resulting oscillation probabilities can be
written using effective mixing parameters. In particular, thePðνμ ↔ νe Þ probability from Eq. (4) can be rewritten byreplacing Δm231 and θ13 by their matter-effective parametersΔm231;M ¼ Δm231
sin2 2θM
13 ¼

XY 2EHmatter − M2j I

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
sin2 2θ13 þ ðΓ − cos 2θ13 Þ2

sin2 2θ13;
sin2 2θ13 þ ðΓ − cos 2θ13 Þ2

ð6Þ

ð7Þ

where Γ ¼ aE=Δm231 . In this form it can be seen that forneutrino energies, matter densities, and Δm231 such thatΓ ∼ cos 2θ13 the effective mixing angle becomes maximal.
This resonant enhancement of the oscillation probability
depends on the sign of the mass hierarchy and occurs for
either neutrinos or antineutrinos through the sign of thematter potential a.
In general atmospheric neutrinos do not traverse constantdensity matter as they travel through the earth, but suchresonant oscillations are nonetheless present. The analyses
presented below use exact three-flavor oscillation probabilities computed including matter effects for varying
matter profiles. Following Ref. [11], the matrix X, whoserow vectors are the propagated mass eigenvectors, can be
written as:j≠k

M2k − M 2j



M2k L
; ð8Þ
exp −i2E

where the M2i =2E are the eigenvalues of the constantdensity matter Hamiltonian Hmatter , and I is the identity
matrix. The oscillation probability can then be written as:Pðνα → νβ Þ ¼ jðUXU† Þαβ j2 :

ð9ÞThe eigenvalues M 2i =2E have been found as Eqs. (21) and(22) of Ref. [11].
An atmospheric neutrino can pass through various
densities of matter on its way to the detector. The
Earth’s atmosphere is modeled as vacuum, and the Earth
as a sphere of radius 6371 km, with a spherical density
profile which is a simplified version of the preliminaryreference Earth model (PREM) [12], as shown in Table I.
The use of the full PREM model with 82 layers provides
no perceptible change in the sensitivity of the SuperKamiokande analysis, so the simplified matter profile isadopted to reduce computation times. To calculate theoscillation probability of a neutrino with energy E produced at a height h above the surface of the Earth, the pathfrom the detector to the neutrino production location istraced through N steps across the atmosphere and different
regions of the Earth’s interior (Fig. 1). Note that because theEarth is modeled as spherically symmetric, this path is afunction of only the production height and zenith angle; it isindependent of azimuthal angle. The oscillation probability
for a given neutrino is calculated by stepping along its path:
 N

Y

Pνα →νβ ðE; h; cos θzenith Þ ¼ UXðLi ; ρi ; EÞU†
i

αβ

2

;

ð10Þ

where Li and ρi are the length and density of the ith step.
Figure 2 shows the νμ survival and νμ → νe transition
probabilities for neutrinos and antineutrinos assuming thenormal mass hierarchy. Resonant oscillation effects areclear in both channels for upward-going neutrinos withenergies between two and ten GeV. In this region matter
effects suppress the disappearance of νμ while enhancing
the appearance of νe . The discontinuity in the oscillation
probabilities for cosine zenith angles steeper than −0.9072001-4

ATMOSPHERIC NEUTRINO OSCILLATION ANALYSIS …

FIG. 1. The propagation of two neutrinos through the simplified model of the Earth used in the analysis below. Both νA and νBare produced in the atmosphere. νA then experiences 6 oscillation
steps (air → crust → mantle → outer core → mantle → crust),
while νB experiences 4 oscillation steps (air → crust → mantle →crust).

corresponds to neutrinos crossing both the outer core andmantle regions of the Earth. For shallower zenith angles the
distortion in the νμ survival probability and the resonant
feature in the νe appearance probability are caused bymatter effects in the mantle region. Note that none of thesefeatures appear in the antineutrino plots. If the inverted
hierarchy were assumed instead, the roles of neutrinos andantineutrinos switch completely and the discontinuities andresonance effects appear with nearly the same magnitude
but in the antinuetrino plots.
III. THE SUPER-KAMIOKANDE DETECTOR
Super-Kamiokande is a cylindrical 50-kiloton water
Cherenkov detector, located inside the Kamioka mine in
Gifu, Japan. An inner detector (ID) volume is viewed bymore than 11,000 inward-facing 20-inch photomultipliertubes (PMTs) and contains a 32-kiloton target volume. The
outer detector, which is defined by the two meter-thick
cylindrical shell surrounding the ID, is lined with reflective
Tyvek to increase light collection to 1,885 outward-facing
eight-inch PMTs mounted on the shell’s inner surface.
Since the start of operations in 1996, Super-Kamiokande
has gone through four data taking periods, SK-I, -II, -III,and -IV.
Though the basic configuration the detector is similar
across the phases there are a few important differences. At
the start of the SK-IV period in 2008 the front-end
electronics were upgraded to a system with an ASIC based

PHYS. REV. D 97, 072001 (2018)

on a high-speed charge-to-time converter [13]. The new
system allows for the loss-less data acquisition of all PMT
hits above threshold and has improved the tagging efficiency of delayed Michel electrons from muon decay from73% in SK-III to 88%.
Further, following a period of detector maintenance andupgrades at the end of SK-I (1996-2001), the implosion of asingle PMT at the bottom of the detector on November 12,
2001, created a shock wave and chain reaction that went on
to destroy 6,665 ID and 1,027 OD PMTs. The detector was
rebuilt the following year with nearly half of the photocathode coverage (19%) in the ID (5,137 PMTs) and thefull complement of OD PMTs for the SK-II period (20022005). Since that time all ID PMTs have been encased infiber-reinforced plastic shells with 1.0 cm thick acrylic
covers to prevent further chain reactions. This resulted in an
increased threshold of 7.0 MeV in SK-II compared to
5.0 MeV in SK-I. In 2006 the detector underwent a second
upgrade in which the remaining ID PMTs were replaced
and additional optical barriers were added to the top andbottom portions of the OD to improve separation with its
barrel region. Both SK-III (2006-2008) and SK-IV (2008present) were operated with the full 40% photocathodecoverage in the ID.
Neutrino interactions which produce charged particles
above the Cherenkov threshold in water are reconstructedbased on the observed ring patterns projected on thedetector walls. Photomultiplier timing information is used
to reconstruct the initial interaction vertex after correcting
for the photon time of flight. Particles are divided into two
broad categories based upon their Cherenkov ring patternand opening angle. Rings from particles which produce
electromagnetic showers, such as electrons and photons,
tend to have rough edges due to the many overlapping rings
from particles in the shower and are labeled e-like or
showering. Muons and charged pions on the other hand,which do not form showers, produce Cherenkov rings withcrisp edges. Such rings are labeled μ-like or non-showering. The event reconstruction assigns momenta to eachreconstructed ring in an event based on the observed
number of photons in the ring. Particles with higher
momenta produce brighter Cherenkov rings. Similarly,
particle directions are inferred based on the shape of theirring pattern. Since the neutrino itself is unobserved, energy
and direction variables for use in the oscillation analysis
described below are based on the properties of theirdaughter particles.
More detailed descriptions of the detector and its
electronics can be found in [13–15].
A. Detector calibration
Over the 20 year history of the experiment changes in the
run conditions have been unavoidable. Seasonal changes inprecipitation and the expansion of underground activities at
the Kamioka site have variable impact on the quality and072001-5

K. ABE et al.

PHYS. REV. D 97, 072001 (2018)

FIG. 2. Oscillation probabilities for neutrinos (upper panels) and antineutrinos (lower panels) as a function of energy and zenith angle
assuming a normal mass hierarchy. Matter effects in the Earth produce the distortions in the neutrino figures between two and ten GeV,
which are not present in the antineutrino figures. Distortions in the νμ survival probability and enhancements in the νe appearance
probability occur primarily in angular regions corresponding to neutrino propagation across both the outer core and mantle regions(cosine zenith < −0.9) and propagation through the mantle and crust (−0.9 < cosine zenith < −0.45). For an inverted hierarchy thematter effects appear in the antineutrino figures instead. Here the oscillation parameters are taken to be Δm232 ¼ 2.5 × 10−3 eV2 ,sin2 θ23 ¼ 0.5, sin2 θ13 ¼ 0.0219, and δCP ¼ 0.

quantity of underground water available to fill the detector
and maintain its temperature. These changes impact thewater transparency and subsequent performance of thedetector and therefore must be corrected through calibrations. Since neutrino oscillations are a function of theneutrino energy, a thorough understanding of the detector
energy scale is important for precision measurements.
At the same time the range of energies of interest to
atmospheric neutrino analysis spans from tens of MeV totens of TeV, eliminating the possibility of calibration
through radioactive isotopes. Accordingly, the energy scaleis calibrated using natural sidebands covering a variety ofenergies. Neutral pions reconstructed from atmospheric

neutrino interactions provide a calibration point via the π 0momentum and stopping cosmic ray muons of various
momenta are used to measure photoelectron production as
a function of muon track length (Cherenkov angle) for
multi-GeV (sub-GeV) energies. Here the muon track lengthis estimated using the distance between the entering vertex
and the position of the electron produced in its subsequent
decay. The energy spectrum of these Michel electrons
additionally serves as a low energy calibration point.
Figure 3 shows the absolute energy scale measurement
using each of these samples.
In the oscillation analysis the absolute energy scale
uncertainty is conservatively taken to be the value of the072001-6

ATMOSPHERIC NEUTRINO OSCILLATION ANALYSIS …

PHYS. REV. D 97, 072001 (2018)( MC - Data ) / Data [%]

4
SK-I

SK-II

SK-III

SK-IV

2

0

-2

-4

Decay electron
Neutral current π0
Sub-GeV stopping μ
Multi-GeV stopping μ

103

102

104

102

103

104

102

103

104

102

103

104

Momentum [MeV]

FIG. 3. Absolute energy scale measurements for each SK period. Vertical error bars denote the statistical uncertainty and horizontal
error bars the momentum range spanned by each analysis.

most discrepant sample from this study in each run period.
The total systematic error is assigned taking this value
summed in quadrature with the time variation of the energy
scale, which is measured using the variation in the averagereconstructed momentum of Michel electrons and thevariation in the stopping muon momentum divided byrange. An example of the latter showing the energy scalestability since SK-I appears in Fig. 4. Note the SK-III
period was subject to poor and volatile water transparency
conditions, resulting in a comparatively turbulent energy
scale. The stability seen in the SK-IV period is a result of
improvements in the water purification system and incorrections for the time variation of the PMT response.
The total energy scale uncertainty in each period has been
estimated as 3.3% in SK-I, 2.8% in SK-II, 2.4% in SK-III,and 2.1% in SK-IV.
B. Sample selection
The current analysis utilizes atmospheric neutrino datacollected during each of the SK run periods andDeviation from mean

0.04

SK-I

SK-II

SK-III

SK-IV

0.02

0

-0.02

-0.04
2000.12.31

2005.12.31

2010.12.31

2015.12.31

FIG. 4. Energy scale stability measured as a function of datesince the start of SK operations. The energy scale is taken as theaverage of the reconstructed momentum divided by range of
stopping cosmic ray muon data in each bin. The vertical axis
shows the deviation of this parameter from the mean value foreach SK period separately. Error bars are statistical.

corresponds to a total livetime of 5,326 days, 2,519 ofwhich are from SK-IV. Super-Kamiokande’s atmospheric
neutrino data are separated into three broad categories, fully
contained (FC), partially contained (PC) and upward-going
muons (Up-μ) that are further sub-divided into the final
analysis samples. Fully contained events have a reconstructed vertex within the 22.5 kton fiducial volume,
defined as the region located more than 2 m from the
ID wall, and with no activity in the OD. The FC data are
sub-divided based upon the number of observed Cherenkov
rings, the particle ID (PID) of the most energetic ring,and visible energy or momentum into combinations of
single- or multi-ring, electron-like (e-like) or muon-like(μ-like), and sub-GeV (E < 1330.0 MeV) or multi-GeV(E > 1330.0 MeV). Additional selections are made basedon the number of observed electrons from muon decays andthe likelihood of containing a π 0 . For the SK-I, -II, and -III
data periods the latter selection is based on [16] and for SKIV it is performed using the improved algorithm presented
in [6]. After all selections there are a total of 14 FC analysis
samples. Events with a fiducial vertex but with energy
deposition in the OD are classified as PC. Based on theenergy deposition within the OD, PC events are further
classified into “stopping” and “through-going” subsamples.
The Up-μ sample is composed of upward-going muon
events produced by neutrino interactions in the rocksurrounding SK or in the OD water. Accordingly, light
deposition in both the OD and ID is expected and thesample is divided into “through-going” and “stopping”
subsamples for events that cross or stop within the ID,
respectively. Through-going events with energy deposition
consistent with radiative losses are separated into a “showering” subsample. The 19 analysis samples defined for eachof the SK run periods are summarized in Table II. Zenith
angle distributions of each sample are shown in Fig. 5.
Distributions of the true neutrino energy for the FC, PC,
and Up-μ event categories appear in Fig. 6. Their event
rates over the lifetime of the experiment have been stable at
8.3 FC events per day, 0.73 PC events per day, and 1.49 Upμ events per day, as shown in Fig. 7. In total 27, 505 μ-like

072001-7

K. ABE et al.

PHYS. REV. D 97, 072001 (2018)

TABLE II. Sample purity broken down by neutrino flavor assuming neutrino oscillations with Δm232 ¼ 2.4 × 10−3 eV2 andsin2 θ23 ¼ 0.5. The data and MC columns refer to the total number of observed and expected events, respectively, including oscillations
but before fitting, for the full 328 kiloton-year exposure. Sub-GeV multi-ring interactions are not used in the present analysis. The
numbers of observed and expected events in this sample are enclosed in parenthesis.
SampleEnergy bins

Fully Contained (FC) Sub-GeV
e-like, Single-ring
0 decay-e
5 e momentum
1 decay-e
5 e momentum
μ-like, Single-ring
0 decay-e
5 μ momentum
1 decay-e
5 μ momentum
2 decay-e
5 μ momentum0
π -likeSingle-ring
5 e momentum
Two-ring
5 π 0 momentum
Multi-ringcos θz bins

CC νe

CC ν¯eCC νμ þ ν¯μ

CC ντ

NC

Data

MC

10 in ½−1; 1
single bin

0.717
0.805

0.248
0.019

0.002
0.108

0.000
0.001

0.033
0.067

10294
1174

10266.1
1150.710 in ½−1; 1
10 in ½−1; 1
single bin

0.041
0.001
0.000

0.013
0.000
0.000

0.759
0.972
0.979

0.001
0.000
0.001

0.186
0.027
0.020

2843
8011
687

2824.38008.7687.0

single binsingle bin

0.096
0.067
0.294

0.033
0.025
0.047

0.015
0.011
0.342

0.000
0.000
0.000

0.856
0.897
0.318

578
1720
(1682)

571.8
1728.4(1624.2)

Fully Contained (FC) Multi-GeV
Single-ringνe -like
4 e momentumν̄e -like
4 e momentum
μ-like2 μ momentum
Multi-ringνe -like
3 visible energy
3 visible energyν̄e -like
μ-like4 visible energy
Other
4 visible energy

10 in ½−1; 1
10 in ½−1; 1
10 in ½−1; 1

0.621
0.546
0.003

0.090
0.372
0.001

0.100
0.009
0.992

0.033
0.010
0.002

0.156
0.063
0.002

705
2142
2565

671.3
2193.72573.810 in ½−1; 1
10 in ½−1; 1
10 in ½−1; 1
10 in ½−1; 1

0.557
0.531
0.027
0.275

0.102
0.270
0.004
0.029

0.117
0.041
0.913
0.348

0.040
0.022
0.005
0.049

0.184
0.136
0.051
0.299

907
745
2310
1808

915.5
773.8
2294.01772.6

Partially Contained (PC)
Stopping
2 visible energy
Through-going
4 visible energy

10 in ½−1; 1
10 in ½−1; 1

0.084
0.006

0.032
0.003

0.829
0.978

0.010
0.007

0.045
0.006

566
2801

570.0
2889.910 in ½−1; 0

0.008

0.003

0.986

0.000

0.003

1456.41448.910 in ½−1; 0
10 in ½−1; 0

0.002
0.001

0.001
0.000

0.996
0.998

0.000
0.000

0.001
0.001

5035.31231.04900.41305.0

Upward-going Muons (Up-μ)Stopping
3 visible energy
Through-goingNon-showering
single bin
Showering
single bin

and 20,946 e-like data events are used in the analysis.
Though events classified as sub-GeV multi-ring interactions are present in the data, they are a small fractionof the available events and provide little additional oscillation sensitivity. As a result they are excluded from thepresent analysis.
As outlined in Sec. II the primary handle for distinguishing the normal from the inverted mass hierarchy is whether
neutrinos or antineutrinos undergo resonant oscillations as
they traverse the earth. The effect of resonant oscillations
would manifest most prominently as an excess of upwardgoing e-like events at OðGeVÞ energies driven by νμ → νeoscillations, so extracting the mass hierarchy requires
separation of νe from ν̄e interactions. As the SK detector
is insensitive to the charge sign of particles traversing thedetector, charged-current (CC) neutrino interactions andantineutrino interactions cannot be differentiated on an
event-by-event basis. Instead this separation is done statistically. It should be noted that due to their larger crosssection and higher flux, more than twice as many neutrino
interactions are expected in the data. Further, while
hierarchy-sensitive matter effects are also present in theνμ → νμ channel, attempts to similarly separate the μ-likedata yielded no significant change in sensitivity and are notconsidered here.
Between two and ten GeV, in addition to charged-current
quasielastic interactions, single-pion (1π) production via Δ
resonance excitation and deep inelastic scattering (DIS)
processes are significant. In the case of the former, theoutgoing π − in antineutrino reactions, such asν̄e þ n → eþ nπ − , will often capture on a 16O nucleus
leaving the positron as the only Cherenkov light-emitting
particle. Neutrino interactions, on the other hand, areaccompanied by a π þ , such as in νe þ n → e− nπ þ , wherethe π þ does not capture in this manner and can therefore
survive long enough to produce a delayed electron throughits decay chain. For CC ν̄e interactions in which the π −has captured there will be no such decay electrons.

072001-8

ATMOSPHERIC NEUTRINO OSCILLATION ANALYSIS …
Sub-GeV e-like 1-dcy e

400

PHYS. REV. D 97, 072001 (2018)

Sub-GeV μ-like 0-dcy e

Sub-GeV e-like 0-dcy e

Sub-GeV μ-like 1-dcy e

400

Up-μ Stopping

400

1000

1000
200

200

200

500

0
1.3

0
1.3

0
1.3

0
1.3

0
1.5

0.7

0.7

0.7

0.7

0.5

Multi-GeV e-like ν e

Sub-GeV π0-like Single-Ring

Multi-GeV e-like ν e

400Up-μ Non-showering

PC Stop
100

1000

50

500

100

200

200

0
2

0
1.5

0
1.3

0
1.5

0
1.3

0

0.5

0.7

0.5

0.7

Sub-GeV μ-like 2-dcy e

200

Multi-Ring e-like ν eMulti-Ring e-like ν e

Up-μ Showering

PC Thru

200

500
200

100

100
100

0
1.3

0
1.3

0
1.5

0
1.3

0
1.5

0.7

0.7

0.5

0.7

0.5

Sub-GeV π0-like Two-Ring

Multi-GeV μ-like

Multi-Ring μ-like

400

-1

-0.5

0

Multi-Ring Other

400
200

500
200

200

0
1.3

0
1.3

0.7

500

0
1.3

0.7

1000

-1

0

momentum (MeV)

1

0
1.3

0.7

-1

cos zenith

0

1

0.7

-1

0

cos zenith

1

cos zenith

cos zenith

FIG. 5. Data and MC comparisons for the entire Super-K data divided into 19 analysis samples. Samples with more than one zenith
angle bin (cf. Table II) are shown as zenith angle distributions (second through fifth column) and other samples are shown asreconstructed momentum distributions (first column). Lines denote the best fit MC assuming the normal hierarchy. Narrow panels below
each distribution show the ratio of the data to this MC. In all panels the error bars represent the statistical uncertainty. In this projection
each bin contains events of all energies, which obscures the difference between the hierarchies. If the inverted hierarchy MC were also
drawn it would lie on top of the normal hierarchy line and for this reason it is not shown here. Figure 10 provides a better projection forcomparing the hierarchies.

Accordingly, an antineutrino enriched subsample isextracted from the single-ring multi-GeV e-like sample
by additionally requiring there are no decay electronspresent. This cut defines the single-ring multi-GeV ν̄e-likesample and its rejected events form the single-ring
multi-GeV νe-like sample. After this selection the fractionsof charged-current electron neutrino and antineutrino
events in the νe -like sample are 62.1% and 9.0%,
250

2500

Events/5000 Days

1400
1200

FC Sub-GeV ( e+ e)

1000

FC Multi-GeV ( e+ e)

Up-μ Stop

FC Sub-GeV ( μ+ μ)

Up-μ Through200

FC Multi-GeV ( μ+ μ)

2000

Up-μ Showering

PC Stop

150

800

1500

600

1000

100

500

50

PC Through

400
200
0
10 1

1

10

102

103

104

Neutrino Energy (GeV)

FIG. 6.

105

0
10 1

1

10

102

103

104

Neutrino Energy (GeV)

105

0
10 1

1

10

102

104

Neutrino Energy (GeV)

True Super-K atmospheric neutrino energy spectra from simulation without oscillations.

072001-9

103

105

K. ABE et al.

PHYS. REV. D 97, 072001 (2018)

SK-I

SK-II

SK-III

180

SK-IV

SK-IV Data
Total MC
Charged Current νeCharged Current νe160
140
Number of Events

Rate [ events / day ]

10

1

FC
PC
Up-μ
10-1
0

1000

120
100
80
60
40
20

2000

3000
4000
5000
Days since April 1 st , 1996

6000

0
-0.6

7000

FIG. 7. Final event rates as a function of time since the start of
SK operations. The error bars are statistical. Circles denote thefully contained event rate and upward-facing (downward-facing)triangles show the partially contained (upward-going muon)
event rates.

respectively. For the ν̄e -like sample the fractions are 54.6%
and 37.2%.
At these energies, events with more than one reconstructed ring are often DIS interactions, which produce
both multiple charged pions and nuclear fragments. In
order to purify the neutrino and antineutrino components ofthe multi-ring samples a two-stage likelihood method has
been developed. Due to the presence of several lightproducing particles the Cherenkov ring produced by theleading lepton is often obscured, resulting in degraded PID
performance and accordingly, significant NC and νμ induced backgrounds in multi-ring events whose mostenergetic ring is e-like. The first stage of the separation
is designed to extract and purify CC νe þ ν̄e interactions
from this base sample. To perform this selection a likelihood function, detailed in a previous publication [16], isbuilt from the PID variable of the event’s most energetic
ring, the fraction of the event’s total momentum it carries,
the number of decay electrons, and the largest distance to a
decay electron vertex from the primary event vertex. The
efficiency of this method for selecting true CC νe þ ν̄eevents is 72.7% and results in a sample that is 73.0% purein these interactions. Separate likelihoods are prepared foreach of the run periods and yield similar efficiencies andpurities. Events that pass this selection are classified as“multi-ring e-like” while those that fail are termed “multiring other.” Though the multi-ring other sample has notbeen used in previous Super-K oscillation analyses it isintroduced here since its νe component offers somehierarchy sensitivity and its oscillation-induced ντ andNC components provide additional constraints on related
systematic uncertainties.
The second stage of the separation process focuses on
separating samples enriched in neutrino and antineutrino
interactions from the multi-ring e-like data. A second

-0.4

-0.2
0
0.2
0.4
0.6
0.8
Multi-Ring νe - νe Separation Parameter1

1.2

FIG. 8. Likelihood distribution used to separate SK-IV Multiring e-like events into the neutrinolike and antineutrinolike
samples. Error bars represent the statistical uncertainty of thedata. Events with negative (positive) likelihood values aredesignated ν̄e -like (νe -like).

likelihood method is introduced based on three variables,
the number of reconstructed rings, the number of decay
electrons, and the event’s transverse momentum. For
charged-current interactions the conservation of charge
implies the total charge of the recoiling hadronic system
must be positive to balance the negative charge of the outgoing lepton. The total charge carried by hadrons emerging
from antineutrino interactions, on the other hand, will be
zero or negative. As a result, the charged pion multiplicity,and hence number of visible Cherenkov rings, in neutrinoinduced events is expected to exceed that from antineutrino
events. This difference is enhanced by the propensity for π −to capture in water. In combination these two effects
suggest that more electrons from the π decay chain are
expected in ν interactions. Due to the V-A structure of the
weak interaction, the angular distribution of the leading
lepton from ν̄ interactions is more forward than those fromν processes. As a result, the transverse momentum of thesystem is expected to be smaller for the former. Since there
is no direct knowledge of an incoming atmospheric
neutrino’s direction the transverse momentum of each
event is defined relative to the direction of the mostenergetic ring. The final likelihood is defined over fivevisible energy bins, 1.33–2.5 GeV, 2.5–5.0 GeV, 5.0–
10.0 GeV, 10.0–20 GeV and > 20 GeV for each SK run
period. Figure 8 shows the combined likelihood distribution used in SK-IV. The efficiency for identifying true CC
ν̄e (νe ) events as ν̄e -like is 71.5% (47.1%).
C. Simulation
The simulation of atmospheric neutrinos is performed
following the flux calculation of Honda et al. [17] andusing the NEUT [18] simulation software (version 5.3.6) togenerate neutrino interactions for tracking in a GEANT3

072001-10

ATMOSPHERIC NEUTRINO OSCILLATION ANALYSIS …

PHYS. REV. D 97, 072001 (2018)[19]-based simulation of the Super-K detector [15]. Several
improvements to NEUT have been made since the previous
version used for atmospheric neutrino analysis (cf. [20]).
Charged-current quasi-elastic interactions are simulatedusing the Llewellyn-Smith formalism [21] with nucleons
distributed according to the Smith-Moniz relativistic Fermi
gas [22] assuming an axial mass MA ¼ 1.21 GeV=c2 andform factors from [23]. Interactions on correlated pairs of
nucleons, so-called meson exchange currents (MEC), havebeen included following the model of Nieves [24]. Pionproduction processes are simulated using the Rein-Sehgal
model [25] with Graczyk form factors [26]. Since the MEC
simulation includes delta absorption processes, the pionless
Δ decay process, Δ þ N 0 → N 00 þ N 0 , in NEUT’s previous
pion production model has been removed in the present
version.
NEUT’s cascade model is used in the detector simulation
to treat the hadronic interactions of pions with nuclei in thedetector. The cross sections underlying the model, including charge exchange, absorption, inelastic scattering, andhadron production processes, have been tuned using a fit toexternal pion scattering data as described in Ref. [27] (c.f.
Table IV). Uncertainties from that fit have been propagatedas systematic uncertainties in the present analysis.
Differences in the expected number of pions in the final
state between the NEUT prediction and measurements fromthe CHORUS experiment [28] are considered as an additional source of systematic uncertainty affecting the event
selection presented above.

On ¼

IV. ATMOSPHERIC NEUTRINO ANALYSIS
Three fits, each incorporating a different degree ofexternal information, are performed to estimate oscillation
parameters. In the first and least-constrained fit, the SuperK atmospheric neutrino data are fit allowing θ13 to vary as afree parameter. The second fit similarly uses only atmospheric neutrino data, but assumes θ13 to be the average ofseveral reactor neutrino disappearance experiment measurements, sin2 θ13 ¼ 0.0219  0.0012 [29]. Finally, theT2K samples discussed in Sec. V are fit alongside theatmospheric neutrino data under the same assumption. In
each of these fits the data are fit against both the normal andinverted hierarchy hypotheses.
Data are fit to the MC using a binned χ 2 method built
assuming Poisson statistics and incorporating systematic
errors as scaling factors on the MC in each bin [30]: X 2X
Oϵi
χ2 ¼ 2En − On þ On ln n þ; ð11Þ
E
σ
n
i
n
i


X
X
i
En ¼
En;j 1 þ
f n;j ϵi
j

i

ð12Þ

ð13Þ

j

In this equation En;j represents the MC expectation in the
nth analysis bin for the jth SK period. Similarly, On;j is thecorresponding data in that bin and f in;j is a coefficientdescribing the fractional change in the bin’s MC under a 1σ ivariation of the ith systematic error source. Systematic
errors penalize the χ 2 based on their corresponding fitting
parameters, ϵi . Solving the system of equations defined by
the requirement ∂χ 2 =∂ϵi ¼ 0 for each systematic error
brings the data and MC into the best agreement allowed by
the systematic errors. This minimization in the systematic
error parameters is repeated over a grid of oscillation
parameters and the parameter set returning the smallest
value of χ 2 is taken as the best fit.
The fit is performed over 520 analysis bins for each of
the SK periods and a total of 155 systematic error sources.
In addition, a systematic error on the presence of mesonexchange currents has been added to the analysis where thedifference between the NEUT model with and without
MEC is taken as the 1σ uncertainty. Further, the single pion
production error of previous analyses has been broken
down into three parts following the parametrization of
Ref. [26]. Systematic errors and their sizes at the best fit
point of the analysis are presented in Tables VII, VIII,
and IX.
When the atmospheric data are studied without external
constraints the fit is performed over four parameters. The
agreement between the data and MC is evaluated using
Eq. (11) at each point in the grid spanned by 0.0 ≤sin2 θ13 ≤ 0.10 (15 points), 0.3 ≤ sin2 θ23 ≤ 0.7 (25 points),
1.0 × 10−3 ≤ jΔm232;31 j ≤ 5.0 × 10−3 eV2 (51 points), and
0.0 ≤ δCP ≤ 2π (19 points). The solar mixing parameters
are set to the values in Table III but their uncertainties aretreated as a source of systematic error in the analysis. For
the normal (inverted) hierarchy fit the fitting parameter is
Δm232 (Δm231 ). Independent fits are performed for thenormal and inverted hierarchies and the grid point returning
the smallest value of χ 2 is termed the best fit for each. The
smallest of these is taken as the global best fit.
Further, the compatibility of the atmospheric neutrino
data with oscillations subject to matter effects in the Earth isevaluated by performing the same fits with sin2 θ13
TABLE III. Values of oscillation parameters fixed in theanalysis and their systematic errors. Note that sin2 θ13 is only
fixed in the “θ13 constrained”analyses described in Sec. V.
Parameterwhere,

X
On;j :

Δm221
2

sin θ12sin2 θ13

072001-11

Valueð7.53  0.18Þ × 10−5 eV20.304  0.014
0.0219  0.0012

K. ABE et al.

PHYS. REV. D 97, 072001 (2018)

constrained to 0.0219  0.0012 (discussed below) andintroducing an additional scaling parameter on the electron
density in Eq. (5), α. This parameter is allowed to range in
20 steps from 0.0 to 1.9, with α ¼ 1.0 corresponding to thestandard electron density in the Earth.
A. Results and discussion

20

20

15

15

Δ χ2

Δ χ2

Figure 9 shows one-dimensional allowed regions forjΔm232;31 j, sin2 θ23 , θ13 and δCP . In each plot the curve isdrawn such that the χ 2 for each point on the horizontal axis
is the smallest value among all parameter sets including thatpoint. When the atmospheric neutrino data are fit bythemselves with no constraint on θ13 , the normal hierarchy
hypothesis yields better data-MC agreement than the
inverted hierarchy hypothesis with χ 2NH;min − χ 2IH;min ¼−3.48. The preferred value of sin2 θ13 is 0.018(0.008)
assuming the former (latter). Though both differ from theglobally preferred value of 0.0219 the constraints are weak
and include this value at the 1σ level. In the normal
hierarchy fit the point at sin2 θ13 ¼ 0.0 is disfavored at
approximately 2σ indicating the data have a weak preference for nonzero values. A summary of the best fit information and parameter constraints is presented in Table V.
The data’s preference for both nonzero sin2 θ13 and thenormal mass hierarchy suggest the presence of upwardgoing electron neutrino appearance at multi-GeV energies
driven by matter effects in the Earth (cf. Fig. 2). Figure 10
shows the up-down asymmetry of the multi-GeV singleand multiring electronlike analysis samples. Here theasymmetry is defined as N U − N D =N U þ N D , whereN U ðN D Þ are the number of events whose zenith angle
satisfy cos θz < −0.4 (cos θz > 0.4). Small excesses seen
between a few and ten GeV in the multi-GeV e-like νe andthe multiring e-like νe and ν̄e samples drive these
preferences.

10

10

99%

5

99%

5

95%
90%

95%
90%

68%

0
0.001

68%

0.002

0.003

0.004

0
0.2

0.005

0.4

0.6

0.8

4

6

sin2 θ23| Δ m232 | , | Δ m231 | eV2
20

20

Inverted Hierarchy

10

10

99%

99%

5

5
95%
90%

95%
90%

68%

0
0

Normal Hierarchy

15

Δ χ2

Δ χ2

15

0.02

0.04

68%

0
0

0.06

2

δCP

sin2 θ13

FIG. 9. Constraints on neutrino oscillation parameters from the Super-K atmospheric neutrino data fit with no external constraints.
Orange lines denote the inverted hierarchy result, which has been offset from the normal hierarchy result, shown in blue, by thedifference in their minimum χ 2 values.

072001-12

ATMOSPHERIC NEUTRINO OSCILLATION ANALYSIS …
The normal hierarchy fits to the atmospheric mixing−3eV2 andparameters yield Δm232 ¼ 2.50þ0.13−0.31 × 10
þ0.036sin2 θ23 ¼ 0.587−0.069 . However, the Super-K data show
a weak preference for the second octant of θ23 , disfavoring
maximal mixing (sin2 θ23 ¼ 0.5) at around 1σ significance.
This preference is driven by data excesses (deficits) at
multi-GeV energies in the upward-going regions of thesingle-ring e-like νe (μ-like) and multi-ring other samples.
These features are consistent with expectations fromνμ → νe oscillations driven by non-zero θ13 .
The best fit value of δCP is found to be 4.18 (3.84) radiansin the normal (inverted) fit, with the least preferredparameter value near 0.8 radians disfavored by Δχ 2 ¼
2.7 (1.0). This preference is driven predominantly by the
sub-GeV e-like samples, via νμ → νe oscillations. Though
the effect of this parameter at these energies is a complicated function of both energy and the neutrino path length,
the point at 4.18 radians generally induces more electron
neutrino appearance in the sub-GeV e-like samples. At
higher energies the effect of δCP modulates the θ13 -driven
νμ → νe probability in the resonance region, but is secondary in size and induces more (less) appearance at 4.18
(0.8) radians. As there are fewer antineutrino events relative
to neutrino events in the atmospheric sample there isaccordingly more freedom to adjust θ13 to bring the MC
prediction into agreement with data in the inverted hierarchy fit. As a result a weaker constraint on δCP is obtained.
The consistency of these data with the presence of matter
effects is illustrated in Fig. 11. With sin2 θ13 set to0.0219  0.0012, the data prefer the normal hierarchy with
an electron density consistent with that of standard matter(α ¼ 1.0). Purely vacuum oscillations, represented by
α ¼ 0.0, are disfavored by the fit by χ 2α¼0 − χ 2min ¼ 5.2
after accounting for the hierarchy uncertainty. Based on toy
Monte Carlo studies, this corresponds to a significance of
excluding vacuum oscillations at 1.6σ.
V. ATMOSPHERIC NEUTRINOS WITH
EXTERNAL CONSTRAINTS
Though the atmospheric neutrino data are sensitive to thevalues of θ13 , θ23 , and jΔm232 j, the size of the mass
hierarchy signal is a function of these parameters. As such,larger uncertainties translate directly into reduced hierarchy
sensitivity. Indeed, toy MC data sets which were generated
with a particular hierarchy but were best fit to the alternative hierarchy often preferred values of the atmosphericmixing parameters different from the input values. For
example, a true normal hierarchy MC generated with θ23 in
the lower octant can be reasonably fit by the inverted
hierarchy hypothesis and the second octant of this parameter. Since there is relatively poor separation between
neutrino and antineutrino interactions, the expected
increase in the event rates in both scenarios is roughly
equal. Restricting the allowed regions of the atmosphericmixing parameters therefore provides increased hierarchy

PHYS. REV. D 97, 072001 (2018)sensitivity by effectively removing such degenerate combinations. The constraints adopted in the present analysis are
based exclusively on information available in the literature
and are described below.
A. Reactor constraint on θ13
Currently the most precise measurements of sin2 2θ13
come from the Daya Bay, RENO, and Double Chooz
experiments In the analysis described below the central
value of this parameter is taken to be sin2 θ13 ¼ 0.0219 
0.0012 based on the average of these measurements
presented in [29]. A systematic error representing the sizeof the uncertainty from this average is incorporated in theanalysis.
B. Constraints from T2K
The T2K (Tokai-to-Kamioka) long-baseline neutrino
experiment sends a beam composed primarily of νμ from
Tokai-village, Japan, 2.5° off-axis toward the SuperKamiokande detector 295 km away. A complex of detectors (the near detectors) located 280 m downstream of theneutrino production point and at the same off-axis angle isused to measure the unoscillated beam spectrum and tothereby constrain the expected spectrum at Super-K (the far
detector). A sharp beam profile peaking at 600 MeV is
expected at the far detector and provides for sensitive
measurements of θ23 and Δm232 . Currently T2K’s measurements [8,27] of these parameters are more constraining thanthe Super-K atmospheric neutrino measurement and provide a statistically-independent constraint. These together
with inherent correlations in some systematic error sources,
such as the detector response and cross section model,
make T2K a powerful input to the Super-K hierarchy
analysis. A more detailed description of the T2K experiment is presented elsewhere [31].
Since Super-K serves as the far detector for T2K many
aspects of the experiments are shared. Notably the detector
simulation as well as the neutrino interaction generator,
NEUT [32], and the event reconstruction tools at Super-K
are common between the two. From the standpoint ofSuper-K then, only the neutrino source and associatedsystematics differ between the beam and atmospheric
neutrino measurements. For this reason it is possible tocreate a reliable simulation of the T2K experiment using
software and methods specific to atmospheric neutrino
measurements, provided only information about the beam
flux and systematic errors. Accordingly, in addition to the
19 × 4 data samples presented in Sec. III, simulated T2K νeappearance and νμ disappearance samples are introduced
into the atmospheric analysis in order to directly incorporate T2K’s measurements. Monte Carlo corresponding tothese samples is constructed from reweighted atmospheric
neutrino MC and data are taken from the literature. This
scheme allows various oscillation hypotheses to be tested
against the published T2K data and in conjunction with the072001-13

K. ABE et al.

PHYS. REV. D 97, 072001 (2018)
0.4
0.2

Up - Down / Up + Down

Data
Normal Hierarchy0

Inverted Hierarchy-0.2
Multi-GeV e-like ν eMulti-GeV e-like ν e

104

104

0.2
0
-0.2
Multi-Ring e-like ν e-0.4

Multi-Ring e-like ν e

104Multi-Ring Other

104
Energy [MeV]104

FIG. 10. Upward- (cos θ < −0.4) to downward-going (cos θ > 0.4) event ratio as a function of energy. The error bars are statistical.
For the single-ring samples the energy is taken to be the visible energy assuming the light-producing particle was an electron. For themulti-ring samples the total energy is used after accounting for the particle type (electron or muon) of each reconstructed ring. The cyanline denotes the best fit from the normal hierarchy hypothesis, and the orange dashed line the best fit from the inverted hierarchy
hypothesis. The error on the prediction is dominated by the uncertainty in the ντ cross section and is not more than 3% (absolute) in any
bin of the figure.

Super-K data. Provided the model samples reproduce
T2K’s results when fit without the atmospheric neutrino
data, the results of a combined analysis can be taken as
reliable.
Neutrino MC samples at Super-K are generated according to the Honda 2011 flux calculation [17] and a sample
equivalent to a 500 year exposure of the SK-IV detector, therun period which contains the T2K beam data, is
10
9Inverted Hierarchy

8

Normal Hierarchy

Δ χ2

7
6
5
4

reweighted according to the beam flux prediction presented
in [33]. Detailed predictions assuming no oscillations areavailable for the νμ , ν̄μ , νe , and ν̄e components of both thebeam and atmospheric fluxes at Super-K. Atmospheric
neutrino interactions are reweighted according to neutrino
flavor, arrival direction, and energy to match the beam
spectrum. Though the T2K beam enters the Super-K tankfrom one direction and atmospheric neutrinos enter from all
directions, the uniformity of the detector’s response is suchthat this reweighting results in negligible biases in themodel samples. Both T2K analysis samples consideredhere are fully contained interactions based on the same
fiducial volume as the atmospheric neutrino sample. The
normalization of the reweighted MC (hereafter beam MC)
is computed based on the total neutrino interaction crosssection on 22.5 kton of water convolved with the beamflux. Table IV lists the interaction rate for 1.0 × 1021

3
2
1
0
0

0.2

0.4

0.6

0.8

1

α

1.2

1.4

1.6

1.8

FIG. 11. Constraints on the matter effect parameter α from theSuper-K
atmospheric
neutrino
data
fit
assuming
sin2 θ13 ¼ 0.0219  0.0012. Orange lines denote the inverted
hierarchy result, which has been offset from the normal hierarchy
result, shown in blue, by the difference in their minimum χ 2
values. Vacuum corresponds to α ¼ 0, while the standard matter
profile used in the rest of the analyses presented here corresponds
to α ¼ 1.

TABLE IV. Expected interaction rates within the SK 22.5 kton
fiducial volume for the T2K beam fluxes (Φ) and cross section
type (σ) presented in [33]. Rates correspond to the number ofinteractions per 1.0 × 1021 protons on target.
ΦσInt=22.5 kton

νμ
νμ
ν̄μ
νe
ν̄e

νe
νμ
ν̄μ
νe
ν̄e

1722.31643.3
53.3
29.3
4.3

072001-14

ATMOSPHERIC NEUTRINO OSCILLATION ANALYSIS …
protons on the T2K target for several combinations ofneutrino flux and cross section.
Separate T2K e-like and μ-like samples are constructedfrom the beam MC using the selection criteria presented in[6,34], respectively. Both samples are composed of fully
contained fiducial volume events with more than 30 MeVofvisible energy and a single reconstructed Cherenkov ring.
To be included in the e-like sample the PID of the ring isrequired to be e-like and must have more than 100 MeV ofvisible energy. Additionally, there must not be any activity

Erec
ν ¼

PHYS. REV. D 97, 072001 (2018)consistent with the electron from a decayed muon and thereconstructed neutrino energy (described below) must be
less than 1250 MeV. A final cut designed to reduce
backgrounds from NC π 0 interactions is applied according
to [6]. Events whose Cherenkov ring has μ-like PID with amomentum greater than 200 MeV=c and at most one decay
electron comprise the μ-like sample.
During the analysis, both samples are binned using
the reconstructed neutrino energy calculated assuming
charged-current quasi-elastic (CCQE) interactions in water:ðMn − V nuc Þ · El − m2l =2 þ Mn · V nuc − V 2nuc =2 þ ðM 2p − M2n Þ=2:
Mn − V nuc − El þ Pl cos θHere Mn ðM p Þ is the neutron (proton) mass and V nuc is theaverage nucleon binding energy in 16O, 27 MeV. The
charged lepton mass, ml , is assumed to be that of anelectron for the e-like sample and that of a muon for the μlike sample. Similarly, the total energy, El , is computed foreach sample using the corresponding ml and the reconstructed momentum, Pl . The cos θ term represents the
opening angle between the neutrino and lepton directions,which is computed using MC truth information for theparent neutrino and the reconstructed direction of thecharged lepton ring. Though the official T2K analyses
use maximum likelihood methods, without detailed information of each data event, reproducing the analyses exactly
using only published information is infeasible. Instead thedata are binned as specified in the T2K publications. The elike sample uses 50 MeV wide bins evenly spaced from 100
to 1250 MeV and the μ-like sample uses 50 MeV bins from
0.2 to 3.0 GeV, 100 MeV wide bins from 3.0 to 5.0 GeV,
and a single bin for more energetic events.
A critical component of the T2K analysis is the constraint coming from measurements of the unoscillated
neutrino flux and interactions at its near detector complex.
Measurements of the CC νμ interaction rate adjust thecentral values and uncertainties on parameters describing
the flux and cross section models underlying the simulation
at Super-K. Incorporation of these constraints alters theshape and composition of the expected spectrum at SuperK and is therefore essential for an accurate reproduction of
the T2K results. Energy dependent normalization parameters for the beam’s νμ , ν̄μ , νe , and ν̄e flux components from[27] are applied as additional weighting factors for the
beam MC. Constraints on the interaction model, such as thevalue of axial mass parameters for quasi-elastic processes
and pion production interactions via the Δ resonance, as
well as the CCQE, CC single pion, and NC π 0 cross section
normalizations are similarly incorporated as multiplicative
weighting factors. For example, the T2K-measured change
in the CCQE axial mass parameter, M QE
A from the default

ð14Þ

value of 1.21  0.45 to 1.33  0.20 is incorporated into thepresent analysis by computing the ratio of the CCQE crosssection for each MC event based on its generated lepton andhadron kinematics. Errors assigned to the flux and crosssection parameters in [27] are used in the construction ofsystematic error response coefficients discussed below. It
should be noted however, that the complete spectral
response of the T2K error model is not publicly available,
and the influence of systematic errors is often expressed as
the expected change in each sample’s event rate. In thesecases the error model used in the atmospheric neutrino
analysis is adapted to produce the same event rate change in
the T2K samples. In the combined analysis of atmospheric
data and the T2K model, detector and cross section
systematic errors are considered completely correlated
between the two data sets, while the flux errors areuncorrelated.
The model constructed here is based on 6.57 × 1020
protons on target taken with T2K’s neutrino-enhanced
beam. Though antineutrino data and contours are available
in the literature (c.f. [35]), the statistics are too low to
impact the sensitivity of the present analysis and are notincluded in the model. Figure 12 shows a comparison of themodel with T2K’s constraints on δCP and the mass
hierarchy after removing (profiling out) the effect of otheroscillation parameters. The expected impact of the T2K
model on the atmospheric neutrino sensitivity to the mass
hierarchy is illustrated in Fig. 13. For all assumed values ofsin2 θ23 , the T2K model’s constraint on the atmosphericmixing parameters strengthens the sensitivity.
It should be noted that other long-baseline neutrino
experiments have made precision measurements of atmospheric mixing parameters, which, when adapted as external
constraints in this analysis, could improve the expected
sensitivity in the same manner as T2K. For example, as
seen in Fig. 15, MINOS [36] constrains Δm232 roughly as
precisely as T2K, although T2K constrains sin2 θ23 better.
Moreover, the neutrino interactions in MINOS are on ironnuclei, not water, introducing an uncancelled systematic

072001-15

K. ABE et al.

PHYS. REV. D 97, 072001 (2018)7

C. Analysis
Inverted Hierarchy, Model
Inverted Hierarchy, T2K Official
Normal Hierarchy, Model
Normal Hierarchy, T2K Official

6

Δ χ2

5
4
3

90%

2
68%

1
0

-3

-2

-1

0

1

2

3

δCP
FIG. 12. Comparison of the T2K model (solid) with T2K
results digitized from Fig. 33 of Ref. [27] (dashed). Cyan andorange lines correspond to the normal hierarchy and inverted
hierarchy fits, respectively, with the offset in the two representing
the difference in their minimum χ 2 values.

uncertainty. Measurements by NOvA [37,38] of muon
neutrino disappearance and electron neutrino appearance
should benefit the present analysis; their inclusion is
anticipated in a future effort.

Incorrect Hierarchy Rejection Δ χ2

6Super-K with T2K νμ , νe (upper)
Super-K (lower)

5

4

3

2

1

0

0.4

0.45

0.5

0.55

0.6

sin2 θ23

FIG. 13. Expected sensitivity to the normal mass hierarchy as afunction of the true value of sin2 θ23 . Grey (lower) and cyan
(upper) bands show the sensitivity of the atmospheric neutrino
sample alone and when combined with the T2K model, respectively. The width of the bands corresponds to the uncertainty inδCP . The inverted hierarchy shows a qualitatively similar improvement in sensitivity and is not shown.

After the introduction of external constraints theatmospheric neutrino data are analyzed in two ways using
modified versions of the fitting scheme outlined in Sec. IV.
In the first analysis the same atmospheric neutrino data
samples and binning are fit over a restricted parameter
space, with sin2 θ13 constrained to 0.0219 as describedabove and other parameter ranges unchanged. An additional systematic error parameter representing the effect of
the uncertainty in external measurements of θ13 on the SK
analysis samples is included in the fit.
The second analysis imposes the same constraint but
introduces additional analysis bins and systematic errors toaccommodate the T2K analysis samples described above.
Using this model of the T2K samples the analysis isperformed over the same oscillation parameter grid anddoes not rely on knowledge of T2K’s published likelihood
surface. Systematic error parameters for the T2K samples
are fit simultaneously with those for the atmospheric
neutrino samples.
D. Results and discussion
Constraints on the atmospheric neutrino mixing parameters and δCP in the θ13 -constrained fit without the T2K
samples are shown in Fig. 14. As in the unconstrained fit
the data prefer the normal hierarchy over the inverted
hierarchy with Δχ 2 ≡ χ 2NH;min − χ 2IH;min ¼ −4.33. Whilethe best fit value of jΔm232 j has shifted slightly, it is withinerrors of the unconstrained fit and in good agreement withother measurements (cf. Fig. 15). Similarly, the preference
for the second octant of θ23 remains unchanged and no
significant change is seen in the width of the parameter’s
allowed region at 1σ. The best fit value of δCP is 4.18 forboth hierarchies, with a tighter constraint on other values
relative to the unconstrained fit. Parameter values and their
1σ errors are summarized in Table V.
In the second fit the addition of the T2K samples is
expected to improve the constraint on the atmosphericmixing parameters due to T2K’s more precise measurements. The left two panels of Fig. 16 show one-dimensional
constraints on these parameters and two-dimensional contours appear in Fig. 17. In the latter dotted lines denote theallowed region from the θ13 -constrained fit to the atmospheric neutrino data only and dashed lines show theallowed regions from the T2K model fit by itself. The
combination of the two data sets, depicted as the solid line,
shows that the fit to these parameters is dominated by theT2K model, with little improvement seen in the contour
when fit together with atmospheric neutrinos.
With less freedom to adjust the atmospheric mixing
parameters, the combination of atmospheric neutrinos with
the T2K model is expected to improve the mass hierarchy
sensitivity on average (see Fig. 13). By itself, the T2K
model favors the normal hierarchy by Δχ 2 ¼ −0.85 [27].
Though T2K has little mass hierarchy sensitivity on

072001-16

ATMOSPHERIC NEUTRINO OSCILLATION ANALYSIS …
20

PHYS. REV. D 97, 072001 (2018)20

20

Inverted Hierarchy

10

10

99%

5

10

99%

5

95%
90%

99%

5

95%
90%

68%

95%
90%

68%

0
0.001

0.002

0.003

0.004

0.005

Normal Hierarchy

15

Δ χ2

15

Δ χ2

Δ χ2

15

0
0.2

68%

0.4

0.6

0.8

0
0

2

4

6

δCP

sin2 θ23| Δ m232 | , | Δ m231 | eV2

FIG. 14. Constraints on neutrino oscillation parameters from the Super-K atmospheric neutrino data fit assuming
sin2 θ13 ¼ 0.0219  0.0012. Orange lines denote the inverted hierarchy result, which has been offset from the normal hierarchy
result, shown in cyan, by the difference in their minimum χ 2 values.

2[eV 2]
Δ m32

0.00350.003

Normal Hierarchy, 90% C.L.
Super-K
T2K
NOvA
IceCube
MINOS

0.0025

0.002
0.2

0.4

0.6

0.8

sin2 θ23

FIG. 15. Constraints on neutrino oscillation contours at the
90% C.L. from analyses assuming the normal mass hierarchy.
The Super-K contour (cyan) is taken from the analysis withsin2 θ13 assumed to be 0.0219  0.0012. Contours from the T2K(violet) [8], NOvA (dashed green) [7], MINOSþ (dashed blue)[36], and IceCube (red) [39] experiments are also shown.

average, Δχ 2 ¼ −0.4 at the Super-K best fit point, this result
is driven by an excess of observed events in its appearance
sample. When atmospheric neutrinos are combined with

T2K, the hierarchy preference strengthens to Δχ 2 ¼ −5.27,
with the majority of the expected sensitivity coming fromthe atmospheric samples appearing in Fig. 10.
Similar preferences in both samples for δCP near 3π=2
result in a stronger constraint on this parameter whenanalyzed together. The right panel of Fig. 16 shows theconstraint for both hierarchy assumptions, with the offset in
the two lines corresponding to the Δχ 2 between the two.
Naturally, this preference is consistent with an increased νe(as opposed to ν̄e ) rate in T2K relative to the expectationfrom the measured value of θ13 . Though the constraint fromthe normal hierarchy fit disfavors the region around π=2,
the contour includes the CP-conserving value δCP ¼ π atnearly 1σ.
VI. INTERPRETATION
It is known that the significance of a mass hierarchy
determination does not necessarily follow the expectationfrom a comparison of the χ 2 minima from each of thehierarchy hypotheses (cf. Ref. [40]). Indeed, the hierarchies
do not form a nested hypothesis and as a result Wilks’
theorem [41] is not applicable. To address the issue of the

TABLE V. Summary of 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 sin2 θ13 , respectively, as
described in the text. The expected absolute χ 2 value for the SK (SK þ T2K) fits is 559.9 (636.2). The p-value for obtaining a smaller χ 2than the data is 0.439 (0.482) in the NH θ13 -constrained fits.
Fit

Hierarchy

SK θ13 Free

NH
IH
NH
IH
NH
IH

SK θ13 Constrained
SK þ T2K θ13 Constrained

χ2

sin2 θ13

sin2 θ23

jΔm232;31 j [×10−3 eV2 ]δCP

571.29574.77571.33575.66639.43644.700.018þ0.029−0.0130.008þ0.017−0.007––––0.587þ0.036−0.0690.551þ0.044−0.0750.588þ0.031−0.0640.575þ0.036−0.0730.550þ0.039−0.0570.550þ0.035−0.0512.50þ0.13−0.31
2.20þ0.33−0.13
2.50þ0.13−0.20
2.50þ0.08−0.37
2.50þ0.05−0.12
2.40þ0.13
−0.05

þ1.45
4.18−1.66
þ2.38
3.84−2.12þ1.41
4.18−1.61
þ1.52
4.18−1.66
þ0.81
4.88−1.48
þ1.05
4.54−0.97072001-17

K. ABE et al.

PHYS. REV. D 97, 072001 (2018)

20

20

20

Inverted Hierarchy

10

10

99%

5

99%

5

95%
90%

68%

0
0.001

10

99%

5

95%
90%

95%
90%

68%

0.002

0.003

0.004

0.005

Normal Hierarchy

15

Δ χ2

15

Δ χ2

Δ χ2

15

0
0.2

68%

0.4

0.6

0.8

0
0

2

4

6

δCP

sin2 θ23| Δ m232 | , | Δ m231 | eV2

FIG. 16. Constraints on neutrino oscillation contours from a combined fit of Super-K atmospheric neutrino data and a model of theT2K experiment assuming sin2 θ13 ¼ 0.0219  0.0012. Orange lines denote the inverted hierarchy result, which has been offset fromthe normal hierarchy result, shown in cyan, by the difference in their minimum χ 2 values.

Normal Hierarchy

0.003

0.00250.00250.002

0.3

Inverted Hierarchy| Δ m231 | MeV2

Δ m232 MeV2

0.003

0.002

0.4

0.5

0.6

0.7

0.3

sin2 θ23

0.4

0.5

0.6

0.7

sin2 θ23

FIG. 17. Constraints on neutrino oscillation contours from a combined fit of Super-K atmospheric neutrino data and a model of theT2K experiment. The left figure shows 90% C.L. constraints from the atmospheric neutrino data (dotted), the T2K model (dashed), andtheir combination (solid) for the normal hierarchy. The right figure shows the same for the inverted hierarchy fit. In each contour sin2 θ13
is constrained to be 0.0219  0.0012. Normal and inverted hierarchy contours for each analysis are drawn relative the best fit among
the two.

hierarchy significance in the present analysis, ensembles ofpseudo data sets generated from the atmospheric neutrino
MC are used to estimate p-values for obtaining a differencein χ 2 between the hierarchy hypotheses more extreme thanthat observed in data. This condition is termed “rejecting”
the alternative hierarchy hypothesis for a given hierarchy
assumption in what follows.
For the Super-K analysis, two important issues need tobe considered. First, as shown in Fig. 13 the expected
sensitivity to the mass hierarchy is a strong function of theunderlying oscillation parameters and as such, p-value
calculations are expected to depend heavily on the parameters assumed in the generation of MC ensembles. Rather
than attempting a Bayesian-like treatment of the p-value
calculation and marginalizing over the effect of eachparameter, a range of p-values has been computed using
the 90% C.L. intervals obtained from the present analysis to
avoid ambiguities surrounding the choice of parameterpriors.
Second, it is also clear from the figure that at the current
level of statistics, Super-K has only modest sensitivity toreject either hypothesis, making the interpretation of the pvalue susceptible to fluctuations of the background. While
the p-value for rejecting the inverted hierarchy (IH)
hypothesis assuming the normal hierarchy (NH) may beunlikely, the p-value in the reverse scenario may be equallyunlikely, leading to an overestimation of the significancewhen stated in terms of the first p-value only. Following thelead of the LHC experiments, this issue is treated using theCLs method [42], where

072001-18

ATMOSPHERIC NEUTRINO OSCILLATION ANALYSIS …

Fraction of Pseudo Data Sets

10-1

10-2

normal hierarchy is driven primarily by upward-goingexcesses seen in hierarchy-sensitive e-like samples,
smaller values of p0 ðIHÞ and larger CLs are obtained
when assuming smaller values of sin2 θ23 or when δCP isnear π=2 since both of these regions predict the least
amount of electron neutrino appearance. For sin2 θ23 >0.60 both metrics decrease as there is sufficient electron
neutrino appearance to discriminate between the two
hierarchy hypotheses at the level seen in the data. In
contrast, both metrics are found to vary only slightly
with Δm232;31 .

Best fit, SK θ13 Constrained
SK True IH
SK True NH
p-value: 0.027
p-value: 0.613

10-3

10-4

-15

-10

-5

0

5

PHYS. REV. D 97, 072001 (2018)10

χ2 - χ2NH

IH

FIG. 18. 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 CLs value for the SK θ13 constrained
analysis. In the cyan (orange) histogram the pseudo data havebeen generated assuming the normal (inverted) hierarchy at theanalysis best fit shown in Table V. Shaded portions of thehistograms denote the fraction of pseudo data sets with more
extreme values than that observed in the data, Δχ 2data ¼ −4.33.

CLs ¼

p0 ðIHÞ
:1 − p0 ðNHÞ

ð15Þ

Here p0 ðIHÞ (p0 ðNHÞ) represents the p-value forobtaining a difference in the minimum χ 2 values between
both hierarchy hypotheses, Δχ 2 ≡ χ 2NH − χ 2IH smaller(larger) than that from the data, Δχ 2data , assuming the truehierarchy is the IH (NH). While CLs does not behave as afully frequentist p-value, it is a conservative method of
preventing erroneous rejection of the null hypothesis when
the overall sensitivity is limited.
MC ensembles were generated assuming statisticalfluctuations of the pseudo data sets according to thecurrent detector exposure, and Gaussian fluctuations of
the systematic errors. Figure 18 shows the distribution of
MC ensembles used in the calculation of the CLs value
for the SK θ13 -constrained fit. Table VI shows the rangeof p-values and CLs values based on ensembles generated with true oscillation parameters taken from the
90% C.L. bounds on θ23 and δCP and best fits fromthe analyses above. Since the data’s preference for theVII. CONCLUSION
Analysis of Super-Kamiokande atmospheric neutrino
data over a 328 kton-year exposure of the detector indicatesa weak preference for the normal mass hierarchy, disfavoring the inverted mass hierarchy at 93.0% assuming oscillation parameters at the analysis best fit point and preferring
matter over vacuum oscillations by 1.6σ. Assuming thenormal mass hierarchy the constraints on the atmospheric2
mixing parameters are sin2 θ23 ¼ 0.588þ0.031−0.064 and Δm32 ¼
þ1.41
2.50þ0.13−0.20 , with δCP ¼ 4.18−1.61 . Fitting in conjunction witha model of the T2K experiment generally enhances these
constraints and the preference for the normal mass hierarchy. Over the range of parameters allowed at 90% C.L.
the inverted mass hierarchy is disfavored by between
81.9% and 96.7% for SK by itself and by between
91.9% and 94.5% when SK is combined with T2K forthe θ13 -constrained fits.
ACKNOWLEDGMENTS
We gratefully acknowledge the cooperation of the
Kamioka Mining and Smelting Company. The SuperKamiokande 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 NRF2009-0083526 (KNRC) funded by the Ministry of Science,
ICT, and Future Planning, the European Union H2020
RISE-GA641540-SKPLUS, the Japan Society for the

TABLE VI. Normal hierarchy significance summarized in terms of the probability of observing a χ 2 preference for the NH more
extreme than that observed in data assuming an IH, p0 ðIHÞ, and CLs values for a range of assumed parameters. The best fit column
reports results assuming MC ensembles generated with oscillation parameters taken from the best fit point parameters from the NH fit ineach analysis. Similarly, the Lower 90% (Upper 90%) column reports values assuming true parameters generated on the lower (upper)
90% C.L. obtained from the fits.
p0 ðIHÞ
Fit
SK θ13 Constrained
SK þ T2K θ13 Constrained

CLs

Lower 90% C.L.

Best Fit

Upper 90% C.L.

Lower 90% C.L.

Best Fit

Upper 90% C.L.

0.012
0.004

0.027
0.023

0.020
0.024

0.181
0.081

0.070
0.075

0.033
0.056

072001-19

K. ABE et al.

PHYS. REV. D 97, 072001 (2018)

Promotion of Science, the National Natural Science Foundation of China under Grants No. 11235006, the National Scienceand 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).
APPENDIX: SYSTEMATIC UNCERTAINTIESTABLE VII. Flux-related systematic errors that are common to all SK run periods. The second column shows the best fit value of thesystematic error parameter, ϵj , in percent and the third column shows the estimated 1σ error size in percent.
Systematic error
Flux normalizationðνμ þ ν̄μ Þ=ðνe þ ν̄e Þ
ν̄e =νeν̄μ =νμUp/down ratioEν < 1 GeV
Eν > 1 GeVb
Eν < 1 GeV
1 < Eν < 10 GeV
Eν > 10 GeVc
Eν < 1 GeV
1 < Eν < 10 GeV
Eν > 10 GeVd
Eν < 1 GeV
1 < Eν < 10 GeV
Eν > 10 GeVe<400 MeVa>400 MeV
Multi-GeV
Multi-ring Sub-GeV
Multi-ring Multi-GeV

Horizontal/vertical ratio

PC<400 MeV>400 MeV
Multi-GeV
Multi-ring Sub-GeV
Multi-ring Multi-GeV
PC

K=π ratio in flux calculationf
Neutrino path length
Sample-by-sample

FC Multi-GeV
PC þ Stopping UP-μ

Matter effects

e-likeμ-like0-decay μ-likee-likeμ-like0-decay μ-likee-likeμ-likee-likeμ-likee-likeμ-likee-likeμ-like0-decay μ-likee-likeμ-like0-decay μ-likee-likeμ-likee-likeμ-likee-likeμ-likeFit value (%)

σ (%)

14.3
7.8
0.08
−1.1
1.6
1.6
3.3
−1.6
0.24
2.9
−2.9
−0.026−0.078−0.286−0.208−0.130−0.442−0.182−0.052−0.104−0.052−0.078−0.052−0.052
0.018
0.018
0.054
0.252
0.341
0.252
0.576
0.414
0.252
0.234
0.504
0.270
0.306
−9.3
−2.13
−6.6
0.22
0.52

25
15
2
3
5
5
5
8
2
6
15
0.1
0.3
1.1
0.8
0.5
1.7
0.7
0.2
0.4
0.2
0.3
0.2
0.2
0.1
0.1
0.3
1.4
1.9
1.4
3.2
2.3
1.4
1.3
2.8
1.5
1.7
10
10
5
5
6.8

Uncertainty decreases linearly with log Eν from 25%(0.1 GeV) to 7%(1 GeV).
Uncertainty is 7% up to 10 GeV, linearly increases with log Eν from 7%(10 GeV) to 12%(100 GeV) and then to 20%(1 TeV)c
Uncertainty linearly increases with log Eν from 5%(30 GeV) to 30%(1 TeV).
dUncertainty linearly increases with log Eν from 8%(100 GeV) to 20%(1 TeV).
e
Uncertainty linearly increases with log Eν from 6%(50 GeV) to 40%(1 TeV).
f
Uncertainty increases linearly from 5% to 20% between 100 GeV and 1 TeV.
ab072001-20

ATMOSPHERIC NEUTRINO OSCILLATION ANALYSIS …

PHYS. REV. D 97, 072001 (2018)

TABLE VIII. Neutrino interaction, particle production, and PMNS oscillation parameter systematic errors that arecommon to all SK run periods. The second column shows the best fit value of the systematic error parameter, ϵj , inpercent and the third column shows the estimated 1σ error size in percent.
Systematic error

Fit value (%)σ (%)

M A in QE
Single π Production, Axial Coupling
Single π Production, CA5
Single π Production, BKG
CCQE cross sectiona
CCQE ν̄=ν ratioa
CCQE μ=e ratioa
DIS cross section
DIS model comparisonsb
DIS Q2 distribution (high W)cDIS Q2 distribution (low W)c
Coherent π production
NC/CCντ cross section
Single π production, π 0 =π Single π production, ν̄i =νi (i ¼ e; μ)dNC fraction from hadron simulationπ þ decay uncertainty Sub-GeV 1-ring

−0.69
−4.4
−3.1
−8.7
6.7
9.2
0.67
−4.4
3.0
8.2
−5.8
−10.0
12.1
−13.8
−20.3
−11.0
−0.47
−0.17
−0.22
1.1
0.25
1.60
−0.2
−1.8
0.022
0.32
0.11

10
10
10
10
10
10
10
5
10
10
10
100
20
25
40
10
10
0.6
0.8
4.1
0.9
5.7
10
10
2.4
4.6
5.4

Final state and secondary interactionse
Meson exchange currentf
Δm2 21 [29]sin2 ðθ12 Þ [29]sin2 ðθ13 Þ [29]e-like 0-decay
μ-like 0-decay
e-like 1-decay
μ-like 1-decay
μ-like 2-decayaDifference from the Nieves [24] model is set to 1.0
Difference from CKMT [43] parameterization is set to 1.0cDifference from GRV98 [44] is set to 1.0
dDifference from the Hernandez [45] model is set to 1.0eError is set by the result of a fit to global data as presented in Ref. [27].
fDifference from NEUT without model from [24] is set to 1.0.
b

TABLE IX. Systematic errors that are independent in SK-I, SK-II, SK-III, and SK-IV. Columns labeled “fit” show the best fit value ofthe systematic error parameter, ϵj , in percent and columns labeled σ shows the estimated 1σ error size in percent.
SK-I

SK-II

SK-III

SK-IV

Systematic Error

Fit Valueσ

Fit Valueσ

Fit Valueσ

Fit Valueσ

FC reduction
PC reduction
FC/PC separation
PC stopping/through-going separation (bottom)
PC stopping/through-going separation (barrel)
PC stopping/through-going separation (top)Non-ν background
Sub-GeV μ-like
Multi-GeV μ-like
Sub-GeV 1-ring
0-decay μ-likePC−0.009
0.016
−0.10
−15.8
3.8
8.5
0.010
0.040
0.010

0.2
2.4
0.6
23
7
46
0.1
0.4
0.1

0.005
−3.43
0.077
−2.4
−5.7
−3.0
0.065
0.065
0.049

0.2
4.8
0.5
13
9.4
19
0.4
0.4
0.3

0.066
−0.012
−0.13
−0.32
−13.9
−12.6
0.105
0.105
0.084

0.8
0.5
0.9
12
29
87
0.5
0.5
0.4

0.68
−0.78
0.0004
−1.5
−0.40
−24.1
−0.011−0.011−0.052

1.3
1
0.02
6.8
8.5
40
0.02
0.02
0.09

0.020

0.2

0.115

0.7

0.381

1.8

−0.282 0.49
(Table continued)

072001-21

K. ABE et al.

PHYS. REV. D 97, 072001 (2018)

TABLE IX. (Continued)
SK-I

SK-II

SK-IIISK-IV

Fit Valueσ

Fit Valueσ

Fit Valueσ

Fit Value

0.068

0.5

0.000

0.2

−0.0040.2−0.000 0.02

0.014

0.1

0.000

0.3

−0.0140.7−0.000 0.08

3.6

13

−5.2

38

−1.0

27

2.6

18

3.7

12

3.8

11

0.75

11

0.34

12

−0.85
2
−0.11
2
0.22
2
−1.5
e-like 0.45
2.3 −1.07 1.3
0.80
2.3
0.96
μ-like 0.14
0.7 −1.91 2.3
1.04
3
1.79
> 400 MeV
e-like 0.078 0.4 −1.40 1.7
0.45
1.3 −0.60
μ-like 0.14
0.7 −0.576 0.7 0.208 0.6 −0.36
Multi-GeV
e-like 0.72
3.7 −2.14 2.6
0.45
1.3 −0.60
μ-like 0.33
1.7 −1.41 1.7
0.35
1
0.72
Multi-ring Sub-GeV
e-like −0.68 3.5
3.13
3.8
0.45
1.3
1.14
μ-like −0.88 4.5
6.75
8.2 −0.90 2.6
1.37
Multi-ring Multi-GeV
e-like −0.61 3.1
1.56
1.9 −0.38 1.1
0.54
μ-like −0.80 4.1 0.658 0.8 −0.73 2.1 −1.43
Particle identification (1 ring)
Sub-GeV
e-like 0.039 0.23 0.227 0.66 0.053 0.26 −0.123
μ-like −0.030 0.18 −0.172 0.5 −0.038 0.19 0.097
Multi-GeV
e-like 0.032 0.19 0.082 0.24 0.062 0.31 −0.154
μ-like −0.032 0.19 −0.089 0.26 −0.060 0.3 0.154
Particle identification (multi-ring) Sub-GeV
e-like −0.23 3.1 −3.446
3.49
9.5 −2.24
μ-like 0.049 0.66 1.38
2.5 −1.91 5.2
0.85
Multi-GeV
e-like 0.48
6.5
5.57
9.7 −1.80 4.9 −1.76
μ-like −0.21 2.9 −2.24 3.9
0.99
2.7
0.85
Multi-ring likelihood selection
Multi-ring e-likeνe , ν̄e −6.5
6.0
−1.3
3.8
−5.3
5.3
−2.3
Multi-ring Other
6.2
5.7
1.4
4.1
4.7
4.9
2.7
Energy calibration−0.75 3.3 −0.90 2.8
0.06
2.4
0.08
Up/down asymmetry energy calibration
0.26
0.6
0.24
0.6
0.74
1.3 −0.15
UP-μ reduction
Stopping−0.091 0.7 −0.090 0.7 0.162 0.7 0.087
Through-going−0.065 0.5 −0.064 0.5 0.115 0.5 0.052
UP-μ stopping/through-going separation
0.003 0.4 −0.004 0.6 0.030 0.4 −0.102
Energy cut for stopping UP-μ−0.043 0.9 −0.122 1.3 0.957
2 −0.122
Path length cut for through-going UP-μ−0.416 1.5 −0.826 2.3 0.993 2.8
1.47
Through-going UP-μ showering separation
7.53
3.4 −4.68 4.4
2.90
2.4 −3.30
Background subtraction for UP-μ Stoppinga
10.0
16
−3.1
21
−4.9
20
−6.7
Non-showeringa
−3.6
18
−3.6
14
1.4
24
2.1
Showeringa−12.3 18 −15.7 14
0.1
24
−0.9
νe =ν̄e Separation−0.98 7.2
6.96
7.9
0.45
7.7
2.46
Sub-GeV 1-ring π 0 selection
100 < Pe < 250 MeV=c
1.7
9
7.0
10
0.98
6.3
5.2
250 < Pe < 400 MeV=c
1.7
9.2
9.8
14
0.76
4.9
3.4
400 < Pe < 630 MeV=c
3.0
16
7.7
11
3.7
24
14.8
630 < Pe < 1000 MeV=c
2.6
14
11.2
16
1.3
8.2
19.4
1000 < Pe < 1330 MeV=c
2.2
12
6.8
9.8
1.7
11
27.4
Sub-GeV 2-ring π 0
1.3
5.6
−2.7
4.4
1.6
5.9 −0.72
Decay-e tagging
−3.2
10
−1.0
10
0.9
10
1.3
Solar Activity
−1.8
20
20.0
50
2.7
20
0.6

2
1.6
3
1
0.6
1
1.2
1.9
2.3
0.9
2.4
0.28
0.22
0.35
0.35
4.2
1.6
3.3
1.6
3.0
3.4
2.1
0.4
0.5
0.3
0.6
1.7
1.5
3
17
17
24
6.8
4.6
3
13
17
24
5.6
10
10

Systematic Error
Sub-GeV e-like(flasher event)Multi-GeV e-like(flasher event)Multi-GeV
1-ring e-like
Multi-GeV
Multi-ring e-like

Fiducial Volume
Ring separation< 400 MeV

The uncertainties in BG subtraction for upward-going muons are only for the most horizontal bin, −0.1 < cos θ < 0.

a072001-22σATMOSPHERIC NEUTRINO OSCILLATION ANALYSIS …
[1] B. Pontecorvo, Sov. Phys. JETP 26, 984 (1968).
[2] Z. Maki, M. Nakagawa, and S. Sakata, Prog. Theor. Phys.
28, 870 (1962).
[3] F. An et al. (Daya Bay Collaboration), Chin. Phys. C 37,
011001 (2013).
[4] J. Ahn et al. (RENO Collaboration), Phys. Rev. Lett. 108,
191802 (2012).
[5] Y. Abe et al. (DOUBLE-CHOOZ Collaboration), Phys.
Rev. Lett. 108, 131801 (2012).
[6] K. Abe et al. (T2K Collaboration), Phys. Rev. Lett. 112,061802 (2014).
[7] P. Adamson et al. (NOvA Collaboration), Phys. Rev. Lett.
118, 151802 (2017).
[8] K. Abe et al. (T2K Collaboration), Phys. Rev. Lett. 118,
151801 (2017).
[9] P. Adamson et al. (NOvA Collaboration), Phys. Rev. Lett.
118, 231801 (2017).
[10] C. Patrignani et al. (Particle Data Group), Chin. Phys. C 40,
100001 (2016).
[11] V. Barger, K. Whisnant, S. Pakvasa, and R. J. N. Phillips,
Phys.
p
ﬃﬃﬃ Rev. D 22, 2718 (1980); A sign error affecting the
2Gf N e term in Eq. 21 of Barger et al. was pointed out in P.
Langacker, J. P. Leveille, and J. Sheiman, Phys. Rev. D 27,
1228 (1983).
[12] A. M. Dziewonski and D. L. Anderson, Phys. Earth Planet.
Inter. 25, 297 (1981).
[13] H. Nishino, K. Awai, Y. Hayato, S. Nakayama, K. Okumura,
M. Shiozawa, A. Takeda, K. Ishikawa, A. Minegishi, and Y.
Arai, Nucl. Instrum. Methods Phys. Res., Sect. A 610, 710(2009).
[14] Y. Fukuda et al. (Super-Kamiokande Collaboration), Nucl.
Instrum. Methods Phys. Res., Sect. A 501, 418 (2003).
[15] K. Abe, Y. Hayato, T. Iida, K. Iyogi, J. Kameda et al.,Nucl. Instrum. Methods Phys. Res., Sect. A 737, 253(2014).
[16] R. Wendell et al. (Super-Kamiokande Collaboration), Phys.
Rev. D 81, 092004 (2010).
[17] M. Honda, T. Kajita, K. Kasahara, and S. Midorikawa, Phys.
Rev. D 83, 123001 (2011).
[18] Y. Hayato, Nucl. Phys. B, Proc. Suppl. 112, 171 (2002).
[19] R. Brun, F. Bruyant, F. Carminati, S. Giani, M. Maire, A.
McPherson, G. Patrick, and L. Urban, Report No. CERNW5013, 1994.
[20] K. Abe et al. (Super-Kamiokande Collaboration), Phys.
Rev. D 91, 052019 (2015).

PHYS. REV. D 97, 072001 (2018)[21] C. H. Llewellyn Smith, Phys. Rep. 3, 261 (1972).
[22] R. A. Smith and E. J. Moniz, Nucl. Phys. B43, 605 (1972);
B101, 547(E) (1975).
[23] R. Bradford, A. Bodek, H. S. Budd, and J. Arrington, Nucl.
Phys. B, Proc. Suppl. 159, 127 (2006).
[24] J. Nieves, J. E. Amaro, and M. Valverde, Phys. Rev. C 70,
055503 (2004); 72, 019902(E) (2005).
[25] D. Rein and L. M. Sehgal, Ann. Phys. (N.Y.) 133, 79(1981).
[26] K. M. Graczyk and J. T. Sobczyk, Phys. Rev. D 77, 053001(2008); 79, 079903(E) (2009).
[27] K. Abe et al. (T2K Collaboration), Phys. Rev. D 91, 072010(2015).
[28] A. Kayis-Topaksu et al. (CHORUS Collaboration), Eur.
Phys. J. C 51, 775 (2007).
[29] K. A. Olive et al. (Particle Data Group), Chin. Phys. C 38,
090001 (2014).
[30] G. L. Fogli, E. Lisi, A. Marrone, D. Montanino, and A.
Palazzo, Phys. Rev. D 66, 053010 (2002).
[31] K. Abe et al. (T2K Collaboration), Nucl. Instrum. Methods
Phys. Res., Sect. A 659, 106 (2011).
[32] Y. Hayato, Acta Phys. Pol. B 40, 2477 (2009).
[33] K. Abe et al. (T2K Collaboration), Phys. Rev. D 87, 012001(2013).
[34] K. Abe et al. (T2K Collaboration), Phys. Rev. Lett. 112,
181801 (2014).
[35] K. Abe et al. (T2K Collaboration), Phys. Rev. D 96, 011102(2017).
[36] P. Adamson et al. (MINOS Collaboration), Phys. Rev. Lett.
112, 191801 (2014).
[37] P. Adamson et al. (NOvA Collaboration), Phys. Rev. Lett.
116, 151806 (2016).
[38] P. Adamson et al. (NOvA Collaboration), Phys. Rev. D 93,
051104 (2016).
[39] M. G. Aartsen et al., Nucl. Phys. B908, 161 (2016).
[40] X. Qian, A. Tan, W. Wang, J. J. Ling, R. D. McKeown, andC. Zhang, Phys. Rev. D 86, 113011 (2012).
[41] S. S. Wilks, Ann. Math. Stat. 9, 60 (1938).
[42] A. L. Read, J. Phys. G 28, 2693 (2002).
[43] A. Capella, A. Kaidalov, C. Merino, and J. Tran Thanh Van,
Phys. Lett. B 337, 358 (1994).
[44] M. Gluck, E. Reya, and A. Vogt, Eur. Phys. J. C 5, 461(1998).
[45] E. Hernandez, J. Nieves, and M. Valverde, Phys. Rev. D 76,
033005 (2007).

072001-23