PHYSICAL REVIEW D 91, 052019 (2015)

Limits on sterile neutrino mixing using atmospheric neutrinos
in Super-Kamiokande
K. Abe,1,29 Y. Haga,1 Y. Hayato,1,29 M. Ikeda,1 K. Iyogi,1 J. Kameda,1,29 Y. Kishimoto,1,29 M. Miura,1,29 S. Moriyama,1,29
M. Nakahata,1,29 Y. Nakano,1 S. Nakayama,1,29 H. Sekiya,1,29 M. Shiozawa,1,29 Y. Suzuki,1,29 A. Takeda,1,29 H. Tanaka,1
T. Tomura,1,29 K. Ueno,1 R. A. Wendell,1,29 T. Yokozawa,1 T. Irvine,2 T. Kajita,2,29 I. Kametani,2 K. Kaneyuki,2,29,*
K. P. Lee,2 T. McLachlan,2 Y. Nishimura,2 E. Richard,2 K. Okumura,2,29 L. Labarga,3 P. Fernandez,3 J. Gustafson,4
E. Kearns,4,29 J. L. Raaf,4 J. L. Stone,4,29 L. R. Sulak,4 S. Berkman,5 H. A. Tanaka,5 S. Tobayama,5 M. Goldhaber,6,*
G. Carminati,7 W. R. Kropp,7 S. Mine,7 P. Weatherly,7 A. Renshaw,7 M. B. Smy,7,29 H. W. Sobel,7,29 V. Takhistov,7K. S. Ganezer,8 B. L. Hartfiel,8 J. Hill,8 W. E. Keig,8 N. Hong,9 J. Y. Kim,9 I. T. Lim,9 T. Akiri,10 A. Himmel,10
K. Scholberg,10,29 C. W. Walter,10,29 T. Wongjirad,10 T. Ishizuka,11 S. Tasaka,12 J. S. Jang,13 J. G. Learned,14 S. Matsuno,14
S. N. Smith,14 T. Hasegawa,15 T. Ishida,15 T. Ishii,15 T. Kobayashi,15 T. Nakadaira,15 K. Nakamura,15,29 Y. Oyama,15
K. Sakashita,15 T. Sekiguchi,15 T. Tsukamoto,15 A. T. Suzuki,16 Y. Takeuchi,16 C. Bronner,17 S. Hirota,17 K. Huang,17
K. Ieki,17 T. Kikawa,17 A. Minamino,17 A. Murakami,17 T. Nakaya,17,29 K. Suzuki,17 S. Takahashi,17 K. Tateishi,17
Y. Fukuda,18 K. Choi,19 Y. Itow,19 G. Mitsuka,19 P. Mijakowski,34 J. Hignight,20 J. Imber,20 C. K. Jung,20 C. Yanagisawa,20
H. Ishino,21 A. Kibayashi,21 Y. Koshio,21 T. Mori,21 M. Sakuda,21 R. Yamaguchi,21 T. Yano,21 Y. Kuno,22 R. Tacik,23,31
S. B. Kim,24 H. Okazawa,25 Y. Choi,26 K. Nishijima,27 M. Koshiba,28 Y. Suda,28 Y. Totsuka,28,* M. Yokoyama,28,29
K. Martens,29 Ll. Marti,29 M. R. Vagins,29,7 J. F. Martin,30 P. de Perio,30 A. Konaka,31 M. J. Wilking,31 S. Chen,32 Y. Zhang,32
K. Connolly,33 and R. J. Wilkes33(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
15High Energy Accelerator Research Organization (KEK), Tsukuba, Ibaraki 305-0801, Japan
16
Department of Physics, Kobe University, Kobe, Hyogo 657-8501, Japan17
Department of Physics, Kyoto University, Kyoto, Kyoto 606-8502, Japan18
Department of Physics, Miyagi University of Education, Sendai, Miyagi 980-0845, Japan
19
Solar Terrestrial Environment Laboratory, Nagoya University, Nagoya, Aichi 464-8602, Japan20
Department of Physics and Astronomy, State University of New York at Stony Brook,
New York 11794-3800, USA
21
Department of Physics, Okayama University, Okayama, Okayama 700-8530, Japan
22
Department of Physics, Osaka University, Toyonaka, Osaka 560-0043, Japan
23
Department of Physics, University of Regina,
3737 Wascana Parkway, Regina, Saskatchewan S4SOA2, Canada
24
Department of Physics, Seoul National University, Seoul 151-742, Korea
25
Department of Informatics in Social Welfare, Shizuoka University of Welfare,
Yaizu, Shizuoka 425-8611, Japan
26
Department of Physics, Sungkyunkwan University, Suwon 440-746, Korea
27
Department of Physics, Tokai University, Hiratsuka, Kanagawa 259-1292, Japan
28
The University of Tokyo, Bunkyo, Tokyo 113-0033, Japan

1550-7998=2015=91(5)=052019(22)052019-1© 2015 American Physical Society

K. ABE et al.

PHYSICAL REVIEW D 91, 052019 (2015)
29

Kavli Institute for the Physics and Mathematics of the Universe (WPI),
Todai Institutes for Advanced Study, University of Tokyo, Kashiwa, Chiba 277-8582, Japan30
Department of Physics, University of Toronto, 60 St. George Street, Toronto, Ontario M5S1A7, Canada
31
TRIUMF, 4004 Wesbrook Mall, Vancouver, British Columbia V6T2A3, Canada
32
Department of Engineering Physics, Tsinghua University, Beijing 100084, China
33
Department of Physics, University of Washington, Seattle, Washington 98195-1560, USA
34
National Centre For Nuclear Research, 00-681 Warsaw, Poland(Received 8 October 2014; published 24 March 2015)
We present limits on sterile neutrino mixing using 4,438 live-days of atmospheric neutrino data from theSuper-Kamiokande experiment. We search for fast oscillations driven by an eV2 -scale mass splitting andfor oscillations into sterile neutrinos instead of tau neutrinos at the atmospheric mass splitting. Whenperforming both of these searches we assume that the sterile mass splitting is large, allowing
sin2 ðΔm2 L=4EÞ to be approximated as 0.5, and we assume that there is no mixing between electron
neutrinos and sterile neutrinos (jU e4 j2 ¼ 0). No evidence of sterile oscillations is seen and we limit jU μ4 j2
to less than 0.041 and jU τ4 j2 to less than 0.18 for Δm2 > 0.1 eV2 at the 90% C.L. in a 3 þ 1 framework.
The approximations that can be made with atmospheric neutrinos allow these limits to be easily applied to
3 þ N models, and we provide our results in a generic format to allow comparisons with other sterile
neutrino models.
DOI: 10.1103/PhysRevD.91.052019

PACS numbers: 14.60.St, 14.60.Pq

I. INTRODUCTION
The flavor oscillations of massive neutrinos have been
well established by a wide range of experiments looking at
the disappearance of neutrinos produced in the atmosphere[1,2], in the Sun [3–8], in nuclear reactors [9–12], and atparticle accelerators [13–15] where recently the appearanceof electron and tau neutrinos were observed in primarily
muon neutrino samples [16–18]. The evidence from theseexperiments suggests two independent neutrino massdifferences, an “atmospheric” Δm2 ≈ 3 × 10−3 eV2 and a“solar” Δm2 ≈ 7 × 10−5 eV2 , requiring that three neutrinos
are participating in oscillations. Experiments at the Large
Electron-Positron collider (LEP) also probed the number ofneutrinos using the width of the Z0 mass peak, which
depends on the number of neutrino flavors into which a Z0
can decay. A combined analysis of all the LEP datameasured 2.980  0.0082 light neutrino families [19].
However, not all neutrino experiments are consistent
with this three-flavor picture; several hints of another,larger mass splitting have appeared. The LSND experimentobserved ν̄e appearance in a ν̄μ beam consistent with twoflavor oscillations with Δm2 ≈ 1 eV2 [20]. The later
MiniBooNE experiment saw some possible signs of ν̄μ →ν̄e as well as νμ → νe oscillations at a similar Δm2 [21].
Additional anomalies appear in experiments looking at
intense ν̄e and νe sources at distances too short for standard
oscillations: a lower rate of ν̄e ’s than predicted was seen atseveral reactor experiments [22,23], and the rate of νe ’s
from 51 Cr and 37 Ar calibration sources at gallium-basedsolar neutrino experiments was 3σ lower than the expected*

Deceased.

rate [24–28]. Both of these hints are consistent withoscillations driven by a Δm2 > 1 eV2 [29]. In order for
the interpretation of these measurements to coexist withthe well-established solar and atmospheric mass splittings,
at least one additional neutrino must be introduced. The
LEP measurements further require that any additional
neutrinos either be heavier than half the Z0 mass, which
would make it difficult for them to participate in oscillations,or not couple to the Z0 boson, and hence not participate
in weak interactions. These noninteracting neutrinos arecalled “sterile.”
Cosmological measurements are also sensitive to thenumber of neutrinos, albeit in a model-dependent way, byidentifying the neutrinos as the relativistic species present
in the early universe. Recent measurements are generally
consistent with an effective number of neutrinos a little
above three, but not excluding four [30–33].
If all of the hints and anomalies are interpreted asconsistent evidence of a single additional sterile neutrino(called the “3 þ 1” model), they require a νμ disappearance
signal with a similar Δm2, which has not been seen in shortbaseline νμ disappearance experiments like CCFR [34] or
MiniBooNE and SciBooNE [35], or in the long-baseline
experiment, MINOS [36]. Consequently, 3 þ 1 models fit
the combined global oscillation data poorly. Theories withadditional sterile neutrinos (3 þ 2, 3 þ 3, 1 þ 3 þ 1) havebeen investigated without a clear consensus interpretationof the experimental data [29,37,38].
The Super-Kamiokande (Super-K, SK) atmospheric data
sample can provide a useful constraint on sterile neutrinos
across a wide variety of proposed sterile neutrino models.
The atmospheric neutrino sample covers a wide range inboth energy, E, and distance traveled, L. The signatures of

052019-2

LIMITS ON STERILE NEUTRINO MIXING USING …

PHYSICAL REVIEW D 91, 052019 (2015)

sterile neutrino oscillations in SK data are valid over arange of mass splittings relevant to previous hints and thelimits set in the 3 þ 1 framework can be readily extended tomodels with more than one sterile neutrino.
II. THE SUPER-KAMIOKANDE EXPERIMENTSuper-Kamiokande is a cylindrical, underground, waterCherenkov detector, 41.4 m in height and 39.3 m indiameter. It is arranged into two optically separated regions.
The inner detector (ID) is instrumented with 11,129
twenty-inch photomultiplier tubes (PMTs) [39] and an
active-veto outer detector (OD) instrumented with 1,885
eight-inch PMTs, both filled with ultra-pure water. A
fiducial volume is defined 2 m from the walls of the ID
and has a mass of 22.5 kton.
Neutrinos are detected by observing the Cherenkov
radiation from the highly relativistic charged particles
produced in neutrino-nucleus interactions. The charged
particles must have a velocity greater than the speed oflight in water, introducing a total energy threshold which
depends on particle mass: 780 keV for electrons, 160 MeV
for muons, and 212 MeV for charged pions. The particles
radiate Cherenkov photons in a cone (42° in water forparticles with velocity close to c) as long as the particle isabove threshold, producing a circular pattern of light which
is projected onto the wall of the detector. Particles whichstop inside the detector produce a ring while those thatexit produce a filled circle. The timing of the Cherenkov
light allows the vertex to be reconstructed, and the
direction of travel of the particle is estimated from thevertex and the Cherenkov ring pattern. More energetic
particles typically produce more total light. Particle types
are identified based on the pattern of the hits making up
the ring. Electrons and photons produce electromagnetic
showers which create many overlapping rings and appear
as a single ring with a fuzzy edge. Nonshowering
particles (muons, pions, protons) produce concentric light
cones as they travel and appear as a single ring with asharp outer boundary.
The neutrino oscillation probability depends on theinitial neutrino flavor, the distance the neutrino travels,L, and the neutrino energy, E. We separate our data into
samples with enhanced νμ or νe flavor content and bin it
using observables correlated with L and E. Instead of
distance, we bin the data in zenith angle, cos θz , defined as
the angle between the event direction and the downward
vertical direction. The neutrinos with the shortest path
lengths are downward-going (cos θz near 1) and theneutrinos with the longest path lengths are upward-going(cos θz near −1). The simulation which predicts the numberof neutrino events in each bin includes a distribution ofneutrino production heights based on a model of theatmosphere described in more detail in [40]. This rangeof production heights introduces a smearing of theoscillation probability for a given zenith angle fordownward-going and horizontal events but is negligible
for upward-going events which cross most of the Earth. For
events with one visible ring, we bin in momentum and for
multiring events we bin in visible energy, defined as theenergy of an electron that would produce the total amountof light observed in the detector.
III. DATA SAMPLE
Super-K has had four run periods, summarized inTable I, with a total exposure 4,438 live-days which areeach considered separately in the simulation and analysis.
The previous atmospheric neutrino oscillation paper[41] included only the first three run periods. The current
period, SK-IV, began with the installation of new front-end
electronics (QTC Based Electronics with Ethernet, QBEE)
whose key component is a new high-speed charge-to-time
converter (QTC) ASIC [42]. The SK-IV data continues tobe accumulated, but this analysis includes only data taken
until September, 2013.
There are three basic event topologies used in theatmospheric neutrino analysis which cover different neutrino energies (plotted in Fig. 1). The fully contained (FC)sample includes events with vertices inside the fiducial
volume and which stop before leaving the inner detector. It
is the lowest-energy sample ranging from a few hundred
MeV up to about 10 GeV. These events have the best energy
resolution since all of the energy is contained within thedetector. However, they also have the worst direction
resolution (from 12° to 100°, depending on energy[43,44]) since the outgoing lepton direction is lesscorrelated with the incoming neutrino direction. In theoscillation analysis, the FC sample is divided into 13
subsamples, categorized based on visible energy into
sub-GeV, below 1.33 GeV, and multi-GeV, above
1.33 GeV. The FC subsamples are then binned in energy
and cos θz , though a few sub-GeV subsamples withparticularly poor direction resolution have only a single

TABLE I. Summary of the four SK data-taking periods. The
photocoverage was reduced during SK-II due to an accident in2001. SK-IV data taking is continuing, but this analysis includes
only data taken until September, 2013. The difference of livetimes between FC/PC and UP-μ is due to the insensitivity of theUP-μ reduction to noise such as “flasher” PMTs. Unlike the UP-μ
reduction, the FC and PC reductions exclude data close in time to
known flashing PMTs to avoid including fake events, reducing
the total livetime for those samples.
Live-days
SK-I
SK-II
SK-III
SK-IV

052019-3

1996–2001
2002–2005
2006–2008
2008–2013

FC/PC

UP-μ

Photocoverage (%)

1,489
799
518
1,632

1,646
828
635
1,632

40
19
40
40

K. ABE et al.

PHYSICAL REVIEW D 91, 052019 (2015)
FC Sub-GeV
FC Multi-GeV
PC Stop
PC Through3000

Events / 4600 Days

2000

1000

250

Up-μ Stop
Up-μ Through200

Up-μ Showering

150
100
50
0 -1
10

1

10

102

103

104

105

Neutrino Energy (GeV)

FIG. 1 (color online). The true energy distribution fromsimulation without oscillations of the fully contained (subGeV and multi-GeV), partially contained (stopping andthrough-going), and up-going muon (stopping, through-goingnonshowering, and through-going showering) samples.

cos θz bin. Details of which bins are used in which
subsample are shown in Table II. The sub-GeV events are
categorized into μ-like, e-like, and neutral-current π 0 -likesamples. The μ- and e-like subsamples are further divided
by number of decay electrons, which can signify thepresence of a charged pion produced below Cherenkov
threshold and thus help isolate neutral-current (NC)backgrounds. The multi-GeV subsamples are split into
μ-like and e-like, with the e-like events divided into νe -likeand ν̄e -like. The FC sample selection techniques aredescribed in greater detail in [41].
The partially contained (PC) sample contains events thathave vertices in the fiducial volume, but produce leptons
that leave the inner detector. They have long tracks and soare almost exclusively from νμ interactions and range inenergy from a few GeV up to tens of GeV. These events
have better direction resolution (9°–16° [43]) than FC
events due to their higher energy, but worse energy
resolution since the exiting muon carries some energy
out of the detector. They are divided into two subsamples
based on their energy deposition in the OD: stopping,which stop in the outer detector, and through-going, whichpass through the outer detector out into the rock [2]. They
are binned in both visible energy (based on light observedin the ID) and cos θz .
Up-going muon (UP-μ) events contain muons that startin the surrounding rock and then enter and pass through theouter detector into the inner detector. This subsample also
starts at a few GeV but extends up to hundreds of TeV.
These events are only included if they are up-going, where
the bulk of the Earth has shielded the detector from theotherwise overwhelming cosmic-ray muon background.
They are split into the lower-energy stopping (stops inthe inner detector) and the higher-energy through-going(exits out the far side of the detector) subsamples. The
through-going events are further subdivided into nonshowering (minimum-ionizing) and showering subsamples
based on the method described in [46]. The critical energy
at which the muon’s energy loss by radiative processes(primarily pair production and bremsstrahlung) equals
energy loss by ionization is 900 GeV [45] so evidenceof showering allows us to select a sample with higher
average energy despite an unknown fraction of the muon’senergy being deposited in the rock before reaching andafter leaving the detector. The UP-μ through-going subsamples are binned only in cos θz since the measured
energy is only a rough lower bound on the initial neutrinoenergy.
A summary of all the event samples used in this analysis,
including the binning used, the number of observed events,
and the number of events predicted by simulation, is shownin Table II.
Several improvements to the simulation have been
included since the last atmospheric neutrino oscillation
publication [41]. The neutrino interaction generator, NEUT[47], includes an updated tau-neutrino cross section and amore accurate calculation of the NC elastic scattering crosssection [48]. This version also includes an improved modelof photon emission from excited nuclei based on recent
experimental data [49,50] and improved spectroscopic
factor simulation [51]. The pion interaction model was
also improved: interaction probabilities were tuned toexisting pion scattering data [52], particularly at low
momentum, < 500 GeV=n, while at higher momenta themodel, including both interaction probabilities and kinematics, was updated to the SAID partial wave analysis ofworld data [53–55]. In all energy regimes, nucleon ejection
after pion absorption in the nucleus was implemented withmultiplicity determined by data in [56] and the kinematicsof the two-body ejection modeled with the data-based
parametrization in [57].
The detector simulation includes a model of the newelectronics and software triggers as well as an updated
tuning of the PMT response in the ID. Improved modelsof the PMT geometry and reflective Tyvek surfaces, aswell as tube-by-tube dark noise rates and saturation curves
based on in situ measurements have been implemented
into the OD simulation. Additionally, low momentum
pion interactions in the water are now simulated using
the pion interaction model from NEUT. The atmospheric
neutrino flux model is taken from [58]. The momentum
reconstruction algorithms have also been updated withsome minor improvements. More details on the event
generator, Monte Carlo simulation (MC) and reconstruction
can be found in [40] and more details on the recent
improvements can be found in [59].

052019-4

LIMITS ON STERILE NEUTRINO MIXING USING …

PHYSICAL REVIEW D 91, 052019 (2015)

TABLE II. Summary of the atmospheric neutrino data and simulated event samples. The oscillated MC has been calculated assuming
three-flavor mixing with Δm232 ¼ 2.51 × 10−3 eV2 , Δm221 ¼ 7.46 × 10−5 eV2 , sin2 ðθ12 Þ ¼ 0.305, sin2 ð2θ13 Þ ¼ 0.095, sin2 ðθ23 Þ ¼
0.514 [3,15,45]. Visible energy is defined as the energy of an electron required to produce all the Cherenkov light seen in the event. The
distribution of 0-, 1-, and 2-decay electron μ-like subsamples changes significantly in SK-IV compared to earlier periods due to theimproved decay-e tagging efficiency of the upgraded electronics. The fraction of UP-μ events classified as showering in the SK-IV datais large relative to SK-I due to the slow increase in the gain of the PMTs over time.
SK-I
Energy bins
Fully contained (FC) sub-GeV
e-like, single ring0 decay-e
5e momentum
1 decay-e
5e momentum
μ-like, single ring0 decay-e
5μ momentum
1 decay-e
5μ momentum
2 decay-e
5μ momentum0
π -likeSingle ring5e momentum
Two-ring
5π 0 momentum

SK-II

SK-III

SK-IV

cos θz bins

Data

MC

Data

MC

Data

MC

Data

MC

10 in ½−1; 1

2987
301

2975.2
310.5

1573
172

1549.1
170.3

1091
118

1052.2
108.8

3074
402

3126.0333.8

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

1025
2012
147

974.1
2042.1
145.4

561
1037
86

534.5
1068.4
76.7

336
742
61

338.1
735.0
60.7

583
2767
245

592.8
2741.2
255.0

181
493

183.6
492.4

111
251

109.1
265.8

59
171

60.7
175.3

194
548

167.7
546.3

Fully contained (FC) multi-GeV
Single ringνe -likeν̄e -like
μ-likeMultiringνe -likeν̄e -like
μ-like4e momentum
4e momentum
2μ momentum

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

191
665
712

170.3
664.4
730.5

79
317
400

82.4
338.2
384.1

68
206
238

59.8
230.3
250.5

238
626
788

221.3
641.3
794.4

3 visible energy
3 visible energy
4 visible energy

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

216
227
603

222.2
224.3
596.4

143
134
337

138.3
132.4
328.7

65
80
228

77.3
76.9
219.6

269
275
694

267.5
264.8
705.3

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

143
759

144.4
777.3

77
350

73.2
370.1

54
290

55.4
306.0

188
919

187.9
948.4

10 in ½−1; 0

432

444.7

206

216.2

194

172.1

416

417.1

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

1564
272

1532.4
325.0

726
110

741.4
117.1

613
110

569.5
142.7

1467
446

1435.8393.1

Partially contained (PC)
Stopping
Through-going

2 visible energy
4 visible energy

Upward-going muons (UP-μ)Stopping
Through-going
Nonshowering
Showering

3 visible energy

H ¼ UMð3þNÞ U † þ V e þ V s :IV. STERILE NEUTRINO PHENOMENOLOGYThe neutrino oscillation probabilities in this analysis are
based on the framework developed in [60]. With N additional sterile neutrinos, the Pontecorvo–Maki–Nakagawa–
Sakata (PMNS) mixing matrix [61,62] must be expanded toa ð3 þ NÞ × ð3 þ NÞ matrix:
0

U e1

BU
B μ1
B
B
U ¼ B Uτ1
B
B Us1
@
..
.

Ue2

U e3

Ue4

Uμ2

U μ3

Uμ4

Uτ2

Uτ3

Uτ4

U s2
..
.

Us3
..
.

U s4
..
.



1

C
C
C
C
C:
C
C
A
..
.

ð4:2Þ

The matrix Mð3þNÞ is the neutrino mass matrix,
Mð3þNÞ ¼

1diagð0; Δm2 21 ; …; Δm2 ð3þNÞ1 Þ;
2E

ð4:3Þ

which also depends on the neutrino energy E. V e and V s
are the potentials experienced by the electron and sterile
neutrinos respectively,
ð4:1Þ

This larger mixing matrix then appears in the completely
generic 3 þ N Hamiltonian,pﬃﬃﬃ
V e ¼ ðGF = 2Þdiagð2N e ; 0; …Þ

ð4:4Þ

pﬃﬃﬃ
V s ¼ ðGF = 2Þdiagð0; 0; 0; N n ; N n ; …Þ

ð4:5Þwhich depend on Fermi’s constant, GF , and the average
electron and neutron densities along the neutrino path, N e

052019-5

K. ABE et al.

PHYSICAL REVIEW D 91, 052019 (2015)

and N n , respectively. Depending on their type, neutrinos
experience one of three different potentials: νe ’s havecharged-current (CC) interactions with electrons andneutral-current (NC) interactions with electrons and nucleons, νμ’s and ντ ’s have only NC interactions, and any νs ’s
have no interactions. The NC interactions depend only on
the neutron density because the Z0 couplings to electronsand protons are equal and opposite and their densities are
identical in neutral matter. The factor of two between N e
and N n comes from the difference between the two currents
in the standard model. We have taken advantage of thefreedom to arbitrarily set the zero of the potential energy todefine lack of NC interactions as a potential experienced by
the sterile neutrinos.
In order to simplify the calculation in the analysis of
atmospheric neutrino data, we introduce a few assumptions.
We assume the sterile mass splittings are sufficiently largethat oscillations in all samples are “fast” and the L=E term,sin2 ðΔm2 L=4EÞ, can be approximated as hsin2 i ¼ 0.5. For
the SK data, this assumption is good for Δm2 > 10−1 eV2.
The complex phases introduced by the additional neutrinos
are also neglected in this treatment since they were shownin [60] to have a negligible impact on the atmospheric
neutrino sample. We also assume that there are no νe − νsoscillations. While hints of these jUe4 j2 -driven oscillations
have been seen in short-baseline νe =ν̄e disappearance [22],SK is not very sensitive to this parameter. We estimate thatallowing nonzero jU e4 j2 values of the size allowed by theseother experiments reduces SK’s sensitivity to jUμ4 j2 by
between 3% and 40%, depending on which experiment isused and what Δm2 value is assumed. In order to avoid
introducing a complex multiexperiment fit, we assume
jUe4 j2 ¼ 0. All of these assumptions are discussed further
in Appendix B.
The assumptions are chosen to eliminate features in theoscillation probability to which the atmospheric neutrinos
are not sensitive and focus only on the parameters that canbe measured. The assumptions do not generally limit theapplicability of the results (e.g. results can be comparedagainst experiments where different assumptions are madein the theory) except in certain specific cases like the valid
range of Δm2 ’s discussed in Appendix B.
We then definedμ ¼

X
jU μi j2

ð4:6Þ

i≥4

and divide the mixing matrix U into a standard neutrino
model part with only the 3 × 3 U PMNS surrounded by zeros~
and an ð3 þ NÞ × ð3 þ NÞ sterile part, U:

U¼

UPMNS
0


0 ~
U:
0

ð4:7Þ

With these assumptions and definitions, we can calculatethe νμ =νe oscillation probabilities following the methodof [60],Pee ¼ P~ ee ;

ð4:8Þ

Peμ ¼ ð1 − dμ ÞP~ eμ ;

ð4:9Þ

Pμe ¼ ð1 − dμ ÞP~ μe ;

ð4:10Þ

Pμμ ¼ ð1 − dμ Þ2 P~ μμ þX
jU μi j4 ;

ð4:11Þ

i≥4where P~ αβ is the probability derived from a three-neutrino
Hamiltonian,~ ¼ U PMNS Mð3Þ U†PMNS þ V e
H000
GF N n X B~ pﬃﬃﬃ@ 0 jU α2 j22 α¼sterile~ α2 U
~ α3
0 U

1
0
~ α2 U~ α3 CUA;

ð4:12Þ

~ α2 jjU2where the first term is the standard neutrino Hamiltonian in
vacuum, the second term is the matter potential in the Earthfrom νe CC interactions, and the third term gives thecomponent of the sterile matter potential which is rotated~into the three active flavors by the sterile mixing matrix U.
The scale of the sterile potential is set by Fermi’s constant
GF and the average neutron density along the neutrino’s
path N n , calculated using the four-layer PREM model ofthe density profile of the Earth [63]. Equations (4.8)–(4.11)
and Eq. (4.12) show that there are two dominant signatures
introduced by sterile neutrino mixing. The first is the
reduction of the νμ survival probability at all lengths andenergies from the ð1 − dμ Þ2 term in Eq. (4.11). The second
signature is the distortion of the oscillation probabilities
when passing through significant amounts of matter due tothe matter effects proportional to N n in Eq. (4.12).
It is not feasible to calculate the oscillation probabilities
generated by generic 3 þ N models since there are too~ by the sum overmany free parameters introduced into H
several α’s: as many as 2N magnitudes and N phases. So,
following the technique of [60], we reduce the parameter
space by introducing further approximations. These
approximations will allow us to perform the fit in thesimpler 3 þ 1 case, described below, and then extend thoseresults into more generic 3 þ N models in Sec. VI.
We examine two approximations, appropriate in different circumstances: the no-νe approximation which assumes
electron neutrinos are fully decoupled from μ − τ − soscillations, and the sterile vacuum approximation which
includes νe appearance via standard three-neutrino oscillations but assumes no sterile matter effects by setting
the neutron density in the Earth to be zero. The former

052019-6

LIMITS ON STERILE NEUTRINO MIXING USING …

PHYSICAL REVIEW D 91, 052019 (2015)

approximation includes both sterile oscillation signatures
but produces a biased estimate of dμ while the latter is only
sensitive to the dμ signature, but produces an unbiased
estimate of it.
Note that in both of the 3 þ 1 approximations, one
explicit parameter
been eliminated because
P has alreadydμ ¼ jU μ4 j2 andjUμi j4 ¼ d2μ ¼ jUμ4 j4 . In the following
sections we will use jUμ4 j2 , but will return to using dμ
in Sec. VI.

where Em are the eigenvalues of the new system, θm is the
new mixing angle, and A32 ¼ Δm232 =4E is the magnitudeof the eigenvalue of the Hamiltonian for two-flavor
oscillations in the atmospheric sector without any sterile
neutrinos.
Pulling together these pieces, the oscillation probabilities
in the no-νe approximation are

A. No-νe oscillation probabilities
The νe ’s are fully decoupled from oscillations by setting
θ13 ¼ θ12 ¼ 0, which allows Eq. (4.12) to be reduced to a
two-level system:
2~ ¼ Δm32
H
4E

− cos 2θ23

sin 2θ23
 ~ 2
GF N n jU s2 jﬃﬃﬃ
 p
~ s2 U~ s3
2
Usin 2θ23ð4:15Þ

As ¼

cos 2θs ¼;

ð4:16Þ

jUτ4 j2 − jUμ4 j2 ð1 − jUμ4 j2 − jU τ4 j2 Þ:ð1 − jUμ4 j2 ÞðjU μ4 j2 þ jUτ4 j2 Þ

ð4:17Þ

The complete system, which is itself Hermitian as the sumof two Hermitian matrices, can also be diagonalized to
produce new effective two-neutrino oscillation probabilities which are a function of the atmospheric mixing
parameters and the sterile parameters above:
E2m ¼ A232 þ A2s þ 2A32 As cosð2θ23 − 2θs Þ;

ð4:18Þ

A32 sinð2θ23 Þ þ As sinð2θs Þ;Emð4:19Þ

sin 2θm ¼
cos 2θm ¼

A32 cosð2θ23 Þ þ As cosð2θs Þ;Emð4:22Þ

Pμμ ¼ ð1 − jUμ4 j2 Þ2 ð1 − sin2 ð2θm Þsin2 ðEm LÞÞþ jUμ4 j4 ;

ð4:23Þ

− ð1 − jU μ4 j2 ÞAs sinð2θs Þ sinð4θm Þsin2 ðEm LÞ

ð4:13Þ

which in the 3 þ 1 model can be expressed in terms of the
only two independent sterile matrix elements, jU μ4 j2 andjUτ4 j2 :ð1 − jU μ4 j2 ÞðjU μ4 j2 þ jUτ4 j2 ÞPeμ ¼ Pμe ¼ 0;

× ð1 − sin2 ð2θm Þsin2 ðEm LÞÞ

Noting that the second matrix is Hermitian, it can be
diagonalized and then parametrized by one real eigenvalue,
As , and one angle, θs ,

GF N n− cos 2θs sin 2θspﬃﬃﬃ As;
ð4:14Þ
sin 2θscos 2θs2

sin 2θs ¼

ð4:21Þ

Pμτ ¼ ð2As þ 2As jU μ4 j2 − jUμ4 j4 − 1Þcos 2θ23~ s2 U~ s3 
U:~ s3 j2
jUðjUμ4 j2 þ jUτ4 j2 Þ;2
qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
2 jU μ4 j2 jUτ4 j2 ð1 − jUμ4 j2 − jU τ4 j2 ÞPee ¼ 1;

ð4:20Þ

þ ð1 − jUμ4 j2 Þð1 þ jUμ4 j2 − 2As Þ;

ð4:24Þ

PNCα ¼ Pαe þ Pαμ þ Pατ ;ð4:25Þ

where PNCα is the probability for a να to remain any active
species and is applied to the NC events in our simulation.
Note that these probabilities depend on jU μ4 j2 , As , and θs ,which in turn depend only on jU μ4 j2 and jU τ4 j2 (plus theatmospheric oscillation parameters). The main signature ofjUμ4 j2 is the reduction in the survival probability of νμ ’s
due to fast oscillations introduced by the ð1 − jUμ4 j2 Þ2coefficient in Eq. (4.23). The primary effect of jUτ4 j2
comes through As which scales the size of the sterile matter
effects in the matter Hamiltonian, as can be seen in
Eq. (4.14). A nonzero jU τ4 j2 also makes θs nonzero, which
enhances the matter effects further. While jU μ4 j2 also
contributes to As , for densities available on Earth the fast
oscillation effect is always much stronger than the matter
effect, so any measurement of jUμ4 j2 will come primarily
from fast oscillations. If, however, sterile matter effects areseen without accompanying fast oscillations, then that must
be caused by jUτ4 j2. The effects of both jU μ4 j2 and jUτ4 j2
on the νμ survival probability vs zenith angle and energy inSK are shown in Fig. 2.
Since the signature of nonzero jU μ4 j2 is a lower νμsurvival probability independent of distance and energy, it
manifests itself in the atmospheric neutrino data as a
reduction in the normalization of all the μ-like samples.
Since there are significant systematic uncertainties on theabsolute neutrino flux but much smaller uncertainties on
the relative flux of νμ ’s to νe ’s, the constraint on the μ-likenormalization depends on the normalization of the e-likesamples. While the νe appearance signal is not very large(approximately 7% of the multi-GeV νe samples),

052019-7

K. ABE et al.

PHYSICAL REVIEW D 91, 052019 (2015)
1.0

1.0

0.8
0.6

0.0
0.4
-0.5

1

10

10

2

0.6
0.0
0.4
-0.5

0.2

-1.0
-1
10

0.8

0.5

cosθz

cosθz

0.5

1.0

1.0

0.0

0.2

-1.0

10-1

1

0.0

102

10

Energy (GeV)

Energy (GeV)

(a)

(b)

FIG. 2 (color online). (a) The νμ survival probability without sterile neutrinos, plotted versus zenith angle and neutrino energy. This
includes the standard νe CC matter effect, which creates the distortion around a few GeV in the most upward-going zenith angles (cos θznear −1), which correspond to neutrinos that pass through the core of the Earth. (b) The νμ survival probability calculated using the no-νeapproximation with jU μ4 j2 ¼ 0.0018 and jU τ4 j2 ¼ 0.33. The distortion due to the νe CC matter effects is gone, but there is now a more
pronounced distortion introduced by the sterile matter effects which reduces the amount of νμ disappearance for the most upward-goingbins in the 10’s of GeV region. There is also a small amount of extra disappearance away from the standard oscillations introduced by thenonzero jU μ4 j2 which is most visible in the slight darkening of the upper-right part of the plot corresponding to the higher energy
downward-going events (cos θz near 1).
ð3Þ

Pμμ ¼ ð1 − jU μ4 j2 Þ2 Pμμ þ jU μ4 j4 ;completely ignoring it does introduce a bias towards lower
measured jU μ4 j2 .
The sterile matter effect signature, on the other hand,changes the shape of the zenith distribution in the PC and
UP-μ samples. Consequently, it is not dependent on the νesamples to control systematic uncertainties and so is notbiased by the no-νe assumption. The sterile matter effects
alter the zenith distribution since the sterile term in
Eq. (4.13) is enhanced by the high average N n experienced
by the most upward-going neutrinos that pass through thecore of the Earth. The distortion is most pronounced in thehigher-energy samples because the large neutrino energy E~
suppresses the standard model part of H.

1

P(νμ → νμ)0.8

Under the alternative sterile vacuum assumption, N n in
~ → HSM and the sterile neutrinos
Eq. (4.12) goes to 0, so H
experience only vacuum oscillations. (This assumption is
also called the “hydrogen-Earth” approximation in [29].)
Then, the P~ αβ terms in Eqs. (4.8)–(4.11) become the
ð3Þ
standard, three-flavor oscillation probabilities, Pαβ , which
are calculated following [64], consistent with previous SK
analyses. Then, the oscillation probabilities can be recalculated as
ð3Þð4:26Þ

Peμ ¼ ð1 − jUμ4 j2 ÞPeμ ;ð3Þð4:27Þ

ð3Þ

ð4:28Þ

Pμe ¼ ð1 − jUμ4 j2 ÞPμe ;

0.6

0.4

0.2

B. Sterile vacuum oscillation probabilities

Pee ¼ Pee ;ð4:29Þ

0
1

10

2

10

3

10

4

10

L/E (km/GeV)

FIG. 3 (color online). The νμ survival probability without sterile
neutrinos (black) and with jU μ4 j2 ¼ 0.058 (red), calculated using
the sterile vacuum oscillation probability, plotted versus L=E.
The oscillation probability is not unique for a given L=E since theνe CC matter effect dependence on L and E is more complicated.
So, many points corresponding to simulated neutrino events areplotted versus L=E, but with oscillation probabilities calculatedusing the individual simulated L and E values, to show the band of
possible oscillation probabilities. While the standard atmosphericoscillation pattern and the smaller variation due to νe CC matter
effects are persistent, introducing a sterile neutrino reduces the
maximum survival probability at all values of L=E. The effect ismost visible in regions without standard atmospheric oscillations.

052019-8

LIMITS ON STERILE NEUTRINO MIXING USING …

PHYSICAL REVIEW D 91, 052019 (2015)Pμτ ¼ ð1 − jUμ4 j2 Þð1 − Pμμ Þ;ð3Þ

ð4:30Þ

PNCα ¼ Pαe þ Pαμ þ Pατ ;ð4:31Þ

ð3Þ

where Pαβ is the standard three-flavor oscillation probability and PNCα gives the survival probability for NC
events. Figure 3 shows the effect of a nonzero jUμ4 j2 on
the νμ survival probability as a function of L=E in theatmospheric sample.
Since the νe appearance is included in this approximation, there is no bias introduced in the estimation of jU μ4 j2 .
However, without the sterile matter effects, there is nosensitivity to jUτ4 j2 .
V. OSCILLATION ANALYSES WITH ONE
STERILE NEUTRINO
The data samples described in Sec. III are fit simultaneously to search for evidence of sterile neutrinos using thesame technique as in [41] with some updates, including
adding the SK-IV data and updating some systematic
uncertainties. Each run period, SK-I, SK-II, SK-III, and
SK-IV, has its own 500 years-equivalent sample of MC toreflect the different physical and operational conditions
during the four run periods.
The oscillation fit minimizes a “pulled” χ 2 [65] which
compares the MC expectation for a particular set ofoscillation parameters with the data based on a Poisson
probability distribution:
XX
X~ ~ϵÞ −
χ2 ¼ 2ðθ;OSKn
E~ SKniiinn

P SKn XOSKn
þ χ 2penalty ð~ϵÞ;þ
Oi ln P n i
SKn
~
~
~
ð
θ;
ϵ
Þ
E
n
n i

ð5:1Þ

where n indexes the four SK run periods, i indexes theanalysis bins, OSKn
is the number of observed events in bini~ ~ϵÞ is the MC expectation in bin i~i during SKn, and ESKn
ðθ;
i
~ and
in SKn with the oscillation parameters being tested, θ,systematic parameters, ~ϵ. The data and expectation are
divided into 480 bins of cos θz and/or energy, depending on
sample, as detailed in Table II. The binning has beenchosen to ensure enough events are in each bin to have astable χ 2 calculation. While the expectation in each bin is
calculated separately for each run period, the four run
periods are summed together for the comparison between
data and MC.
The effects of the systematic errors on the expectation
are approximated as linear:


X
ϵ
SKn j
;
1þ
f i;j
σj
j

~ ~ϵÞ ¼ ESKn ðθÞ~E~ SKn
ðθ;
i
i

ð5:2Þ

~ is the MC
where j indexes the systematic errors, ESKnðθÞi
expectation in bin i in SKn without systematic shifts, andf SKni;j is the fractional change in bin i in SKn due to σ j , the1-sigma change in systematic j. The constraints on theseparameters are included as a penalty term in Eq. (5.1):
χ 2penalty ð~ϵÞ ¼

X ϵj 2
j

σj

:

ð5:3Þ

The two analyses consider 155 systematic error parameters; some of them are common to all four SK run periods
and some are calculated separately and treated as independent for each period. The common errors include
uncertainties in the atmospheric neutrino flux, neutrino
interaction cross sections, particle production within nuclei,
and the standard PMNS oscillation parameters. They comefrom the Honda flux calculation [58], external neutrino
interaction measurements as well as model comparisons,and other oscillation measurements, respectively. For these
uncertainties, f SKni;j is the same in SK-I, SK-II, SK-III, and
SK-IV. The period-specific errors are generally related todetector performance: uncertainties on reconstruction, particle identification, energy scale, and fiducial volume differ
between run periods since they depend on the specificgeometry and hardware of the detector, which are determined using control samples and simulation studies. For
these uncertainties, f SKni;j will be nonzero in one run period
and zero in all the others. All the systematic uncertainties
and their sizes are listed in Appendix D.
Equation (5.1) is minimized with respect to ~ϵ for eachchoice of θ~ in a fit’s parameter space. A set of linear
equations in the ϵj ’s are derived from Eq. (5.1) using thefact that the derivative ∂χ 2 =∂ϵi is zero at the minimum [65].
These equations can then be solved iteratively to find theminimum profile likelihood for that set of oscillation~ The best fit
parameters, building up a map of χ 2 vs θ.
point is defined as the global minimum of this map. Tests
performed with high-statistics simulation, both without andwith simulated sterile neutrino signals, showed no significant biases in the extracted best fit points.
In order to focus the analysis on the sterile neutrino
parameters, the standard oscillation parameter values wereconstrained to external measurements and their uncertainties
taken as systematic uncertainties. The T2K measurementof νμ disappearance, jΔm232 j ¼ ð2.51  0.10Þ × 10−3 eV2and sin2 ðθ23 Þ ¼ 0.514  0.055 [15], is used because its
narrow-band beam makes it less sensitive to the sterile
effects being measured in this analysis. The mixing angle
sin2 ð2θ13 Þ ¼ 0.095  0.01 is taken from the PDG world
average [45], the solar terms are taken from the global fit
performed by the SK solar þ KamLAND analysis, Δm221 ¼
ð7.46  0.19Þ × 10−5 eV2 , sin2 ðθ12 Þ ¼ 0.305  0.021 [3],and we assume Δm232 > 0 and δcp ¼ 0, though the precise

052019-9

K. ABE et al.

PHYSICAL REVIEW D 91, 052019 (2015)
PC Through-going

1

400

0.8
300

0.6
τ4

|U |2

200
100
0
-1.0

0.4
-0.5

0.0

0.5

1.0

0.2

cos θzUp-μ Stopping

0
-3
10

250

-1

10

|U |

1

2

μ4

200

FIG. 5. The 90% and 99% upper limits on jU τ4 j2 vs jU μ4 j2 fromthe no-νe fit are shown by the solid and dashed lines, respectively.
The best fit point is marked by a black dot. The light gray region
is excluded at 90% and the dark gray region is disallowed byunitarity.

150
100
50
0
-1.0

-2

10

-0.8

-0.6

-0.4

-0.2

0.0

cos θzFIG. 4 (color online). Zenith angle distributions summed across
SK-I through SK-IV of the PC through-going and UP-μ stoppingsubsamples shown for the data (black points with statisticalerror bars), the MC prediction without sterile neutrinos (blacksolid line), and the MC prediction with a large (approximately 5σ
sensitivity) sterile signal of jU τ4 j2 ¼ 0.31. Both MC predictions
are shown after fitting the systematic uncertainties to the data.
These two subsamples are shown because they contain the 10’s of
GeV neutrinos most sensitive to the sterile matter effect. The
prediction with a sterile component shows an up-turn for the mostup-going events which corresponds to the distortion of theoscillogram shown in Fig. 2(b).

value of these choices have negligibly small effects on thisanalysis.
A. No-νe analysis
As described in Sec. IVA, the analysis with the no-νeapproximation fits both jUμ4 j2 and jU τ4 j2 . Since it does notinclude normal νe matter effects it is systematically biased
towards smaller jUμ4 j2 values than the CC matter effect fit.
The fit is performed on a two-dimensional grid of 200
points, 50 points in jUμ4 j2 distributed logarithmicallybetween 10−3 and 10−1 and 40 points in jUτ4 j2 distributed
linearly between 0 and 0.4.
The best fit is at jUμ4 j2 ¼ 0.012 and jUτ4 j2 ¼ 0.021 with2
χ min ¼ 531.1 over 480 bins (goodness-of-fit p ¼ 0.05).
Figure 4 show the zenith angle distributions for the
subsamples most sensitive to the jU τ4 j2 parameter andan example of what a large sterile contribution would look

like. The Δχ 2 to the no-sterile prediction is 1.1, consistent
with no sterile neutrinos at the 1σ level with two degrees
of freedom. We limit jU τ4 j2 to less than 0.18 at 90% andless than 0.23 at 99%. These limits are independent of the
new Δm2 above 0.1 eV2 (see Appendix B 3). The contours
in jUτ4 j2 vs jUμ4 j2 can be seen in Fig. 5. The jUμ4 j2 best fit
point and limit are discussed in the next section in theanalysis which focuses on that parameter.
B. Sterile vacuum analysis
The analysis with the sterile vacuum approximation fits
only jU μ4 j2, the term which drives fast oscillations, creating
extra disappearance at all energies and zenith angles in all
μ-like samples. The fit is performed on a one-dimensional
grid of 200 points distributed logarithmically between 10−3
and 10−1 . The best fit is at jUμ4 j2 ¼ 0.016 with χ 2min ¼
532.1 over 480 bins (goodness-of-fit p ¼ 0.05). No sterile
oscillations is slightly disfavored by Δχ 2 ¼ 1.1.
Figure 6 shows the best fit from this analysis in several
μ-like samples which closely matches the prediction without sterile neutrinos (the ratio is approximately unity acrossall bins). In fact, there is no net difference in χ 2 between thebest fit point and the prediction without sterile neutrinos
looking just at the difference between the data and theprediction in each bin. All of the difference in χ 2 at the best
fit come from the reduction in the systematic penalty termfrom introducing a nonzero jU μ4 j2 . The dashed line in
Fig. 6 shows the prediction with the best fit sterile
parameter, but without separately minimizing the systematic uncertainties. It shows the effect of a nonzero jUμ4 j2 inisolation: it lowers the normalization in the μ-like samples

052019-10

LIMITS ON STERILE NEUTRINO MIXING USING …

PHYSICAL REVIEW D 91, 052019 (2015)

Sub-GeV μ-like 1 decay-e

1.0

TABLE III. The best fit pull values, shown for both no sterile
neutrinos and the best fit point from the sterile vacuum analysis,of the systematics which change the most between those twopoints. The values at the best sterile fit are all significantly smaller
than the values assuming no sterile neutrinos, reducing the χ 2
penalty term.

0.9

Systematic uncertainty

0.8
-1.0

ðνμ þ ν̄μ Þ=ðνe þ ν̄e Þ, < 1 GeVðνμ þ ν̄μ Þ=ðνe þ ν̄e Þ, 1–10 GeV
CCQE νμ =νe

Ratio to No Steriles

1.2
1.1

-0.5

0.0

0.5

cos θz1.0

Multi-GeV μ-like
Ratio to No Steriles

1.1
1.0
0.9

-0.5

0.0

cos θz0.5

1.0Up-μ Showering
1.2

Ratio to No Steriles

Best fit (σ)

−0.49
−0.50
0.36

−0.13
−0.09
0.01

is concentrated in three systematic errors: the ðνμ þ ν̄μ Þ=ðνe þ ν̄e Þ ratio in the atmospheric flux below 1 GeVand from
1–10 GeV as well as the charged-current quasielastic(CCQE) νμ =νe cross-section ratio, summarized in Table III.
All three of these systematic errors relate to the relative
normalization between the μ-like subsamples, which havesterile oscillations, and the e-like subsamples, which donot. These two flux systematics specifically affect the lowenergy subsamples and the CCQE interaction mode isdominant at lower energies, so it affects the same subsamples. The flux uncertainty is calculated as part of theneutrino flux model, which uses direct muon flux measurements plus simulations of hadronic interactions in theatmosphere constrained by hadron production experiments[58]. The uncertainty is between 2% and 3% in size at theseenergies. The CCQE cross-section uncertainty comes from

1.2

0.8
-1.0

No steriles (σ)

1.1
1.0

10

2

0.9
0.8
-1.0

-0.8

-0.6

-0.4

-0.2

0.0

FIG. 6 (color online). Ratios to the MC prediction withoutsterile neutrinos, binned in zenith angle and summed across SK-I
through SK-IV, for three μ-like subsamples at low (FC sub-GeV),
medium (FC multi-GeV), and high energies (though-going UPμ). The prediction without sterile neutrinos has been fit to the datausing the systematic uncertainties. The black points represent thedata with statistical error bars and the solid red line shows the MC
prediction with the best fit for sterile neutrinos (jU μ4 j2 ¼ 0.016),
including the best fit systematic uncertainties. In all the samples itlines up close to unity, meaning the prediction is nearly identicalto the prediction without sterile neutrinos. The dashed red line
shows the MC prediction with the same sterile component(jU μ4 j2 ¼ 0.016), but now with the same systematic uncertainty
parameters as the denominator, showing the effect of just thesterile oscillations: the normalization is shifted downward by
approximately 3% in every μ-like sample.

by approximately 3%. By introducing this normalization
change with the sterile oscillation parameter, several systematic error parameters can be moved closer to theirnominal values, reducing the χ 2 penalty term. The reduction

Δm241 (eV2)cos θz

10

1

-1

10

10

-3

10

-2

10

-1

|U |2
μ4

FIG. 7. The 90% and 99% upper limits on jU μ4 j2 from thesterile vacuum fit to Super-K is shown in the solid and dashed andvertical lines, respectively. The gray filled region is excluded at
90%. This analysis is not sensitive to Δm2 , but the experiments
who also measure jU μ4 j2 are, so here the one-dimensional SuperK result is shown in two dimensions. The dotted line is the 90%
limit placed by the joint analysis of MiniBooNE and SciBooNE[35] and the dot-dashed line is the 90% limit placed by the CCFR
experiment [34].

052019-11

K. ABE et al.

PHYSICAL REVIEW D 91, 052019 (2015)

the difference between the default model in NEUT [66] anda local Fermi gas model [67] and is 1% to 1.5% in size.
While the sterile oscillations create effects in basicallyevery μ-like sample, these low-energy samples are the mostimportant in this analysis since they have the highest
statistics and thus the smallest statistical uncertainties.
We limit jU μ4 j2 to less than 0.041 at 90% and less than
0.054 at 99%. These limits are independent of the new Δm2
above 0.1 eV2 (see Appendix B 3) and can be compared toother limits on sterile-driven νμ disappearance from shortbaseline experiments in Fig. 7. The limits on this parameterare dominated by the systematic uncertainties on the lowenergy normalization and the sensitivity improvement withincreased statistics will be relatively small unless bettersystematic constraints are included. The expected sensitivity to this parameter is a limit at 0.024 at 90%, somewhat
tighter than the observed limit since it assumes a best fit
with no sterile neutrino component.
VI. EXTENDING THE ANALYSES TO
ADDITIONAL STERILE NEUTRINOS
The oscillation probabilities from Sec. IV, both the no-νeand sterile vacuum approximations, were developed toallow extensions to multiple sterile neutrinos.
A. Extending the sterile vacuum analysis
The most straightforward extension is with the sterilevacuum analysis. Starting again with the oscillation probabilities from Eqs. (4.8)–(4.11), we perform the substitution
ð3Þ
P~ αβ → Pαβ , but leave the probabilities in terms of dμ ,Precalling that dμ ¼ jUμi j2 for i ≥ 4:
ð3Þ

Pee ¼ Pee ;ð6:1Þ

Peμ ¼ ð1 − dμ ÞPeμ ;

ð3Þ

ð6:2Þ

ð3Þ

ð6:3Þ

Pμe ¼ ð1 − dμ ÞPμe ;ð3Þ

Pμμ ¼ ð1 − dμ Þ2 Pμμ þX
jU μi j4 ;

ð6:4Þ

i≥4
ð3Þ

Pμτ ¼ ð1 − dμ Þð1 − Pμμ Þ;ð6:5Þ

PNCα ¼ Pαe þ Pαμ þ Pατ :ð6:6Þ

In [60], the authors note that these expressions are almost
equivalent to Eqs. (4.26)–(4.29) with jU μ4 j2 → dμ , except
P
for the constant term
jUμi j4 from Eq. (6.4), which does
not equal dμ 2 due to potential cross terms. Following their
method, we can write the νμ survival probability asð3Þ

Pμμ ¼ ð1 − dμ Þ2 Pμμ þ d2μ ð1 þ ξ2μ Þ=2;

ð6:7Þwhere ξμ parametrizes the second order deviation from d2μ
in the constant term introduced by additional sterile
neutrinos. Their studies show the effect of ξμ on the dμ
limit from atmospheric neutrinos in the context of a 5νmodel and show that it has no significant effect on the limit
placed on dμ [68].
The independence from ξμ derives from how the scalingP
term ð1 − dμ Þ2 and constant termjU μi j4 affect theoscillation probability. The primary effect we observe inatmospheric neutrinos comes from the scaling term which
creates extra disappearance that is independent of baseline
and energy since it scales the entire νμ survival probability.
This effect is visible almost everywhere in the atmosphericð3Þdata, except where Pμμ → 0 (see Fig. 3). The constant term
creates an opposing, but smaller effect which reduces
disappearance but it is usually overwhelmed by the scaling
term. The effect of the constant term can only be seen in the
ð3Þ“valleys” of the oscillation probability where Pμμ → 0. In
atmospheric neutrinos, the bottoms of these valleys are notclearly resolved, so this effect is vanishingly small (asopposed to in long-baseline experiments which precisely
measure the first oscillation minimum). Since the value ofξμ can be neglected when performing this fit, the oscillation
probabilities in Eqs. (6.1)–(6.6) are in fact equivalent tothose in Eqs. (4.26)–(4.29).
Due to this equivalence, the limit on jU μ4 j2 from the3 þ 1 fit shown in Fig. 7 can be taken as the limit on dμ in
general 3 þ N models.
B. Extending the no-νe analysis
The results from the no-νe analysis can also be extended,at least in an approximate way, to theories with additional
sterile neutrinos. Recall that the oscillation probabilities in
Eqs. (4.21)–(4.25) depend on the solutions to the two-level
Hamiltonian,
F Nn~ ¼ HSM  Gp
ﬃﬃﬃ Hs :
H
2

ð6:8Þ

With additional sterile neutrinos, the dependence canbecome quite complicated since there is a sum overmultiple sterile species, α, from Eq. (4.12):~ α2 j2 U~ α2 U~ α3 
X  jU
Hs ¼:
ð6:9Þ
~ α2 j2~ ~jU
α¼sterile U α2 U α3
In the most general case, Hs depends on 3N free parameters(two magnitudes and a phase difference for each sterile
species α). However, this matrix is 2 × 2 and Hermitian, so
no matter how many independent terms go into the matrix,
there can only be two free parameters after diagonalization,an eigenvalue we labeled As and a mixing angle θs .
~ first in terms of a genericIn Sec. IVA, we rewrote H
diagonalized Hermitian matrix parametrized by As and θs ,052019-12

LIMITS ON STERILE NEUTRINO MIXING USING …

PHYSICAL REVIEW D 91, 052019 (2015)

9

1

8
0.8

7

sin(2θs)

6
0.6

5
4

0.4
3
2

0.2

1
0
10

-3

10

-2

10

-1

0

As
FIG. 8 (color online). The Δχ 2 from the fit to the atmospheric
neutrino data in the no-νe approximation, plotted versus the two
effective parameters, As and sinð2θs Þ, with the third freeparameter, dμ , profiled out.

and then calculated those parameters by explicitly diagonalizing Eq. (6.9) with only one sterile neutrino species.
~ depend only on the two free
After that, the solutions of H
sterile parameters in the mixing matrix, jU μ4 j2 and jU τ4 j2 ,and thus the oscillation probabilities in Eqs. (4.21)–(4.25)
depend only on those parameters as well.
To constrain models with additional sterile neutrinos, we
perform a fit using the same oscillation probabilities, but
we do not solve explicitly for Eqs. (4.15)–(4.17), meaning
the oscillation probabilities, and hence the χ 2 surface
produced by the fit, depend on the two generic parameters,
As and θs , plus dμ on which the oscillation probabilities in
Eqs. (4.21)–(4.25) have an explicit dependence (we havesubstituted dμ for jUμ4 j2 as described in the previous
section). The values of these parameters can be calculated
easily from the sterile part of the mixing matrix U for any
sterile neutrino theory, and they can then be used to look up
the Δχ 2 from this atmospheric fit, allowing constraints to be
put on the parameters in that theory. See the supplemental
material for a table containing the full three-dimensionaldelta log likelihood surface [69]. As a demonstration, theΔχ 2 surface for As vs sinð2θs Þ (dμ has been profiled out) isshown in Fig. 8.
The μ → τ and hence μ → s probabilities are approximations in this case since they neglect some potential
cross-terms introduced in the sum over α, but the fit isdominated by the νμ disappearance signal, so this approximation in the NC and τ oscillation probabilities will havelittle effect on the results.

oscillations with an additional sterile neutrino. The fit was
performed with two different approximations appropriate
for setting limits on the two new matrix elements in the 3 þ 1
framework to which Super-K is sensitive: jUμ4 j2 and jU τ4 j2 .
No significant evidence for fast oscillations driven by a newlarge Δm2 or of the matter effect associated with oscillations
from νμ to νs are seen. We limit the 3 þ 1 parameters jU μ4 j2
to less than 0.041 and jUτ4 j2 to less than 0.18 at 90%. While
the measurement of jU μ4 j2 is limited by systematic uncertainties on the neutrino flux and cross section around 1 GeV,
the constraint on jUτ4 j2 can potentially improve with additional atmospheric data. Assuming only a single sterile
neutrino, these new limits increase the known tension
between νμ disappearance measurements and the hints seenin the νe appearance and disappearance channels. Since
these limits are independent of the size of the new Δm2, theyexclude some new regions of parameter space at low mass
splittings where beam experiments are not sensitive. They
can also be extended readily to 3 þ N models which might
resolve the tensions between the three channels, and theresults are provided in a format to allow tests of more generalmodels in the supplemental materials [69].
ACKNOWLEDGMENTS
The authors would like to thank M. Maltoni for his
help calculating and implementing the sterile oscillation
probabilities. The authors gratefully acknowledge thecooperation of the Kamioka Mining and Smelting
Company. Super-K has been built and operated from funds
provided by the Japanese Ministry of Education, Culture,
Sports, Science and Technology, the U.S. Department ofEnergy, and the U.S. National Science Foundation. This
work was partially supported by the Research Foundationof Korea (BK21 and KNRC), the Korean Ministry of
Science and Technology, the National Science Foundationof China, the European Union FP7 (DS laguna-lbno PN284518 and ITN invisibles GA-2011-289442) the National
Science and Engineering Research Council (NSERC) of
Canada, and the Scinet and Westgrid consortia of Compute
Canada.
APPENDIX A: RESULTS IN OTHER
PARAMETRIZATIONSThere are several mostly equivalent parametrizations thatcan be used for the sterile oscillation parameters in 3 þ 1models. While we have chosen to present our results interms of the magnitude of the matrix elements, we present
in Tables IV and V the limits in some other choices ofparameters:
sin2 θ24 ¼ jUμ4 j2 ;ðA1ÞVII. CONCLUSION

sin2 θ34 ¼ jUτ4 j2 =ð1 − jU μ4 j2 Þ;ðA2ÞThe atmospheric neutrino data from all four periods ofSuper-Kamiokande have been fit to look for evidence ofsin2 2θμμ ¼ 4jU μ4 j2 ð1 − jU μ4 j2 Þ;ðA3Þ

052019-13

K. ABE et al.

PHYSICAL REVIEW D 91, 052019 (2015)

jU s4 j2 ¼ 1 − jUμ4 j2 − jUτ4 j2 :ðA4Þ

TABLE IV.

90% C.L.’s from the sterile vacuum fit.

jU μ4 j2

sin2 θ24

θ24

sin2 2θμμ

0.041

0.041

7.7°

0.071

APPENDIX B: ASSUMPTIONS IN THE
OSCILLATION MODEL
A number of assumptions and approximations are madein order to make the 3 þ 1 calculations easier and to allow
those results to be extended to more general 3 þ N models.
This Appendix presents the justification for the validity of
three of the major assumptions.

TABLE V. 90% C.L.’s from the no-νe fit. The profiled value ofjU μ4 j2 ¼ 0.010 for this point.
jU τ4 j2

jU s4 j2

sin2 θ34

θ34

0.18

0.81

0.18

25°

1. No sterile-electron neutrino mixing
Following the method in Appendix C 2 of [60], we canapproximate the primary effect of a nonzero jUe4 j2 byconsidering only its effect on the νe survival probability
Pee , taken as analogous to Pμμ :ð3Þ

Pee ¼ ð1 − jU e4 j2 Þ2 Pee þ jU e4 j4 ;

ðB1Þ

ð3Þ

where Pee is the standard three-flavor νe survival probability. When this extra free parameter is introduced, thelimit on jUμ4 j2 turns out to be correlated with the limit on
jUe4 j2 , as shown in the sensitivity contours in Fig. 9. With
jUe4 j2 unconstrained, the expected 90% limit on jU μ4 j2
becomes 0.067, 180% larger than the 0.024 90% sensitivity
limit with the assumption of jU e4 j2 ¼ 0. However, once
constraints from other experiments are introduced the effectis significantly reduced. The [60] paper introduces aconstraint of jU e4 j2 < 0.012 at the 1σ level based on a
0.1

0.08

e4

|U |2

0.06

0.04

0.02

0

0

0.02

0.04

0.06

0.08

0.1

2

|U |
μ4

FIG. 9. The 90% sensitivity contour for the sterile vacuum fit
with the effect Pee from Eq. (B1) included. Allowing the freedom
in the electron sample normalization reduces the sensitivity tojU μ4 j2 as can be seen from the bowing outward on the right sideof the contour. Note that on this plot jU μ4 j2 is shown in linear
scale so the correlation with jU e4 j2 is clear.

value from the Bugey [70] limit around Δm2 ¼ 1 eV2 .
Applying this constraint to this analysis leads to a 17%change in our sensitivity. In the low-Δm2 region, where ourresults are most competitive, the change is only 3%, whileat the highest Δm2 ’s the change can be as large as 30%. If
instead the nonzero hints from global fits are used as
constraints, the change in our limit ranges from 10% to
40%, with the larger effects again occurring at higher Δm2.
A proper accounting of these constraints would require a
global fit to multiple experiments introducing Δm2 andjUe4 j2 as fit parameters, which is beyond the scope of thisanalysis; instead we take the approach used in [60] and theatmospheric section of [29] and assume jUe4 j2 ¼ 0.
2. No three-flavor matter effects in the no-νe fit
The main effect of setting θ13 to zero in the no-νe fit,eliminating multi-GeV νe appearance, was already discussed in Sec. VA. However, this assumption has a second
effect: it eliminates the distortion in the νμ survival
probability from matter effects in the Earth. These distortion can be seen in the few-GeV region for the mostupward going events (cos θz ≊ − 1) in Fig. 2(a).
Neglecting this matter effect turns out to have little effecton the jUτ4 j2 limit. A sensitivity fit using the no-νe model toa MC prediction made using the full three-flavor oscillation
probability which includes these distortions finds a best fit
at jUτ4 j2 ¼ 0 and jUμ4 j2 equal to its minimum value (it isbinned in log scale and so does not go to zero). The threeflavor distortions in the νμ survival probability turn out tobe relatively small (at most a few percent in the PC throughgoing and stopping UP-μ samples) and to not affect thethrough-going UP-μ samples which are distorted significantly by the sterile matter effects.
3. Sterile-induced fast oscillations
This assumption posits that the oscillations driven byΔm2 are so fast that individual oscillation periods cannot be
resolved in the experiment and that functions of Δm2 canbe replaced with their average values:

052019-14

LIMITS ON STERILE NEUTRINO MIXING USING …

PHYSICAL REVIEW D 91, 052019 (2015)

0.6

0.6

FC Sub-GeV μ-like
0.55

0.55

5% error on <sin2><sin2>

Δm2=0.045 eV2<sin2>FC Multi-GeV μ-like

0.5

5% error on <sin2>0.5

Δm2=0.055 eV2
0.45

0.45

0.4
-2
10

10

-1

1

10

Δ m (eV )
(a)
2

2

10

0.4
-2
10

2

10

-1

0.6

Δ m (eV )
(b)
2

10

10

2

0.6

Up-μ

PC
0.55

0.55

5% error on <sin2><sin2><sin2>

1

2

0.5

5% error on <sin2>0.5

Δm2=0.13 eV2

Δm2=0.1 eV2

0.45

0.45

0.4
-2
10

10

-1

1

Δ m (eV )
(c)
2

2

10

10

0.4
-2
10

2

10

-1

1

Δ m (eV )
(d)
2

2

10

10

2

FIG. 10 (color online). The average of sin2 ðΔm2 L=4EÞ calculated event-by-event for a range of values of Δm2 in four SK samples:(a) FC sub-GeV, (b) FC multi-GeV, (c) PC, and (d) UP-μ. When the event-by-event average deviates significantly (here defined as 5%)
from hsin2 i ¼ 0.5, the “fast-oscillation” assumption is no longer valid.

 2 
Δm L
→ hsini ¼ 0;
sin
4E

ðB2Þ

 2 
Δm L
1
→ hsin2 i ¼ :
4E
2

ðB3Þ

sin2

However, since the phase in these terms depends on L andE as well as Δm2 , the ranges over which they are valid
could vary for the different samples used in the analysis.
For a sufficiently small Δm2 , this fast oscillation
assumption will break down, and the higher the energy
and shorter the path length, the larger of a value of Δm2 thatis invalid. We can estimate this lower limit by calculating
the value of sin2 ðΔm2 L=4EÞ for many MC events in thevarious SK samples (FC sub- and multi-GeV, PC, and UPμ) at a range of possible values of Δm2 . The average is thencalculated from the event-by-event values at each Δm2 andthe point where the actual average deviates significantly
from one half can be found. These averages vs Δm2 for thefour samples can be seen in Fig. 10.
Setting a threshold of 5% error on the value of hsin2 i,
we find that the fast oscillation sample is valid until

approximately 10−1 in all four samples. The highest
limit is 0.13 in the PC sample where there are both highenergies and the very short track lengths from down-going
events.
Meeting this assumption only sets the bottom of thevalid Δm2 range. The upper limit on the mass forwhich the limits are valid is set by the requirement thatthe mass splitting is sufficiently small that the neutrinos remain coherent. A sufficiently heavy neutrino,
approximately 1 keV or so, will separate from theother light neutrinos and thus not be able to participate
in oscillations.
APPENDIX C: ZENITH ANGLE AND
MOMENTUM DISTRIBUTION
Below are shown the zenith angle or energy distributions, summed across SK-I through SK-IV, for all thesamples in the analysis. Figure 11 shows the μ-like FC, PC,
and UP-μ subsamples while Fig. 12 shows the e-like andNCπ 0 -like samples. For subsamples binned in both zenith
angle and energy, the projection into only zenith angle isshown. The plots show the data represented by points with

052019-15

K. ABE et al.

PHYSICAL REVIEW D 91, 052019 (2015)
Sub-GeV μ-like no decay-e

Sub-GeV μ-like 1 decay-e
800

300

Sub-GeV μ-like 2 decay-e
150

600
200

100

400
100

0
-1.0

50

200

-0.5

0.0

0.5

1.0

cos θz

0
-1.0

0

-0.5

0.0

0.5

1.0

500

1000

1500

μ momentum (MeV/c)cos θzMulti-Ring μ-like

Multi-GeV μ-like
250

300

SKI-IV Data

200
No Sterile

200

150

100

0
-1.0

-0.5

0.0

cos θz

0.5

1.0

100

No-νe Best Fit

50

Hydrogen-Earth Best Fit

0
-1.0

400

60

300

40

200

20

100

0.0

cos θz0.5

1.0

PC Through-going

80

-0.5

0.0

cos θzPC Stopping

0
-1.0

-0.5

0.5

1.0

0
-1.0

-0.5

0.0

cos θz0.5

1.0Up-μ Non-Showering

Up-μ Stopping

Up-μ Showering

800

250
200
150

600

150

400

100

200

50

100
50
0
-1.0 -0.8 -0.6 -0.4 -0.2

cos θz0.0

0-1.0 -0.8 -0.6 -0.4 -0.2

cos θz0.0

0-1.0 -0.8 -0.6 -0.4 -0.2

0.0

cos θzFIG. 11 (color online). Distributions of zenith angle or energy, summed across SK-I through SK-IV, for the μ-like FC, PC, andUP-μ subsamples. They are projected into zenith angle when binned in both angle and energy and the sub-GeV 2 decay-e sample isbinned only in momentum. The black points represent the data with statistical error bars, while the solid blue line represents the noνe best fit, the dashed red line represents the sterile vacuum best fit, and the solid black line represents the MC prediction withoutsterile neutrinos.

052019-16

LIMITS ON STERILE NEUTRINO MIXING USING …

PHYSICAL REVIEW D 91, 052019 (2015)
Multi-Ring νe-like

Multi-GeV νe-like

Sub-GeV e-like no decay-e
1000

100

120

800

80

100

600

60

400

40

200

20

80
60

0
-1.0

-0.5

0.0

0.5

1.0

40
20

0
-1.0

-0.5

0.0

0.5

1.0

0
-1.0

-0.5

Multi-GeV νe-like

Sub-GeV e-like 1 decay-e
300

300

0.0

0.5

1.0

cos θzcos θzcos θz140

Multi-Ring νe-like

120
100

200

200

80
60

100

100

40
20

0

500

1000

1500

e momentum (MeV/c)

0
-1.0

-0.5

0.5

1.0

0
-1.0

-0.5

0.0

0.5

1.0

cos θzcos θz2-ring π01-ring π0
200

500

150

400

SKI-IV Data
No Sterile

300

100
50
0

0.0

500

1000

1500

e momentum (MeV/c)

200

No-νe Best Fit

100

Hydrogen-Earth Best Fit

0

200

400

600

800 1000π0 momentum (MeV/c)

FIG. 12 (color online). Distributions of zenith angle or energy, summed across SK-I through SK-IV, of the e- and NCπ 0 -like FC
subsamples. They are projected into zenith angle when binned in both. In the fits for sterile neutrinos, these samples serve primarily tocontrol the normalization of the atmospheric neutrino flux. As in Fig. 11, the black points represent the data with statistical error bars,
while the solid blue line represents the no-νe best fit, the dashed red line represents the sterile vacuum best fit, and the solid black line
represents the MC prediction without sterile neutrinos. The small deviation of the no-νe best fit in the multi-GeV and multiring e-likesamples is from setting θ13 to zero.

statistical error bars as well as the best fits from the twoanalyses (no-νe as solid blue and sterile vacuum as dashedred) as well as the MC prediction without sterile neutrinos(represented by a black line), with systematic uncertainties
still fit to the data.
The best fits generally agree quite well with the prediction without sterile neutrinos, though both fits favor anonzero sterile oscillation component because it allows forsome systematic uncertainties to fit closer to their nominal
values, as discussed in Sec. IV B.

APPENDIX D: SYSTEMATIC UNCERTAINTIES
Tables VI–VIII summarize the best fit systematic error
parameters for the best fit point from the sterile vacuum
analysis. The pull values for any given systematic error are
expected to be normally distributed across many experiments. For this particular data set, the pull values approximately follow a Gaussian distribution, though with a width
narrower than one due to the interaction between the lowenergy normalization uncertainties and the sterile parameter jU μ4 j2 described in Sec. V B.

052019-17

K. ABE et al.

PHYSICAL REVIEW D 91, 052019 (2015)

TABLE VI. Flux-related systematic errors that are common to all SK run periods. The flux uncertainties come from the Honda flux
calculation [58] and are themselves based on the external data sets used as inputs to the calculation. The second column shows the best
fit value of the systematic 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 ratio

Fit value (%)
Eν < 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 ratioPC< 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 effectsae-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μ-like

21
1.7
−0.25
−0.26
6.7
2.5
2.6
2.6
0.021
1.9
4.2
−0.0037−0.011−0.041−0.029−0.018−0.063−0.026−0.0074−0.015−0.0074−0.011−0.0074−0.0074
0.011
0.011
0.033
0.15
0.21
0.15
0.35
0.25
0.15
0.14
0.31
0.17
0.19
1.3
0.094
−5.8
0.79
1.8

σ (%)
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.
b052019-18

LIMITS ON STERILE NEUTRINO MIXING USING …

PHYSICAL REVIEW D 91, 052019 (2015)

TABLE VII. Neutrino interaction, particle production, and PMNS oscillation parameter systematic errors that are common to all SK
run periods. These uncertainties come primarily from comparisons between different cross section models and external neutrino
interaction measurements. The neutrino oscillation parameter errors come from the cited measurements. The second column shows thebest fit value of the systematic error parameter, ϵj , in percent and the third column shows the estimated 1-σ error size in percent.
Systematic error
M A in QE and single π
CCQE cross sectiona
CCQE ν̄=ν ratioa
CCQE μ=e ratioa
Single meson production cross section
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Δm2 32 [15]sin2 ðθ23 Þ [15]Δm2 21 [3]sin2 ðθ12 Þ [3]sin2 ð2θ13 Þ [45]e-like 0-decay
μ-like 0-decay
e-like 1-decay
μ-like 1-decay
μ-like 2-decay

Fit value (%)

σ (%)

−6.4
1.8
18
0.12
14
2.2
−1.5
0.003
−3.1
1.8
9.8
−4.6
−35
−11
−3
−0.48
−0.64
3.3
0.71
4.5
2
2.8
0.079
0.42
−0.55

10
10
10
10
20
5
10
10
10
100
20
25
40
10
10
0.6
0.8
4.1
0.9
5.7
3.98
10.9
2.55
6.89
10.5

a

Difference from the Nieves [67] model is set to 1.0.
Difference from CKMT [71] parametrization is set to 1.0.
c
Difference from GRV98 [72] is set to 1.0.
dDifference from the Hernandez [73] model is set to 1.0.
b

TABLE VIII. Systematic errors that are independent in SK-I, SK-II, SK-III, and SK-IV. The detector uncertainties are determinedusing control samples like cosmic ray muons and 2-ring π 0 ’s, and simulation studies. Columns labeled “fit” show the best fit value of thesystematic error parameter, ϵj , in percent and columns labeled σ shows the estimated 1-σ error size in percent.
SK-I

SK-II
σ

Fit value σ

Systematic error

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 μ-like
PC
Sub-GeV e-like
Multi-GeV e-like
Multi-GeV 1-ring e-like
Multi-GeV multiring e-like
Fiducial volume

0.006 0.2
0.007 0.2
−0.99
2.4 −3.47
4.8−0.027 0.6
0.081 0.5
−22.4
23
0.2
13
1.88
7
−5.54
9.4
8.3
46
−3.3
19
0.009 0.1
0.009 0.1
0.036 0.4
0.009 0.1
0.009 0.1
0.009 0.1
0.018 0.2
0.062 0.7
0.016 0.5
0.003 0.2
0.003 0.1
0.002 0.1
3.3
13 −15.0
38
1.1
12
2.5
11
−0.04
2
0.08
2

SK-IIIFit value

SK-IV
σ0.038 0.8
−0.041 0.5
0.003 0.9
−0.2 12−9.0 2916.0
87−0.009 0.1−0.009 0.1−0.018 0.2
−0.16
1.8
−0.003 0.1−0.013 0.4
5.1
27
−6.1
11
−0.42
2

Fit value

σ

0.030
0.3
−0.24 10.0001 0.02
−1.06 6.8
−0.65 8.5
−3.3 40−0.026 0.1−0.026 0.1−0.211 0.8−0.129 0.49
−0.000 0.1−0.000 0.1
1.1 183.1 120.40
2(Table continued)

052019-19

K. ABE et al.
TABLE VIII.

PHYSICAL REVIEW D 91, 052019 (2015)(Continued)
SK-I

Systematic error
Ring separation

Fit value< 400 MeVσ

Fit value σ

SK-IIIFit value

SK-IVσ

Fit valueσe-like 1.07
2.3 −1.09
1.3
0.79
2.3
0.05
1.6
μ-like 0.324 0.7 −1.93
2.3
1.03
3
0.09
3
e-like 0.185 0.4 −1.43
1.7
0.44
1.3 −0.03 1
μ-like 0.324 0.7 −0.588 0.7
0.205 0.6 −0.018 0.6
e-like 1.71
3.7 −2.18
2.6
0.44
1.3 −0.03 1
μ-like 0.79
1.7 −1.43
1.7
0.34
1
0.04
1.2
e-like −1.62
3.5
3.19
3.8
0.44
1.3
0.06
1.9
μ-like −2.08
4.5
6.88
8.2 −0.89
2.6
0.07
2.3
e-like −1.44
3.1
1.59
1.9 −0.38
1.1 0.027 0.9
μ-like −1.90
4.1
0.671 0.8 −0.72
2.1 −0.07 2.4
e-like 0.016 0.23 0.099 0.66 0.023 0.26 −0.025 0.28
μ-like −0.013 0.18 −0.075 0.5 −0.016 0.19 0.020 0.22
e-like 0.013 0.19 0.036 0.24 0.027 0.31 −0.031 0.35
μ-like −0.013 0.19 −0.039 0.26 −0.026 0.3 0.031 0.35
e-like −0.31
3.1 −3.39
6
5.09
9.5
2.15
4.2
μ-like 0.066 0.66 1.45
2.5 −2.79
5.2 −0.80 1.6
e-like 0.64
6.5
5.54
9.7 −2.63
4.9
1.71
3.3
μ-like −0.29
2.9 −2.24
3.9
1.43
2.7 −0.80 1.6
0.00
1.1 −0.20
1.7
0.65
2.7 −0.36 2.3
0.293 0.6 −0.070 0.6
0.36
1.3 −0.109 0.3−0.185 0.7 −0.131 0.7
0.111 0.7 0.126 0.5−0.132 0.5 −0.094 0.5
0.080 0.5 0.075 0.3
0.007 0.4
0.016 0.6
0.034 0.4 −0.109 0.6
0.085 0.9
0.11
1.3
0.87
2
0.01
1.7
0.86
1.5
1.50
2.3 −0.12
2.8 −1.87 1.5
3.59
3.4 −2.84
4.4
2.35
2.4 −4.88 3
10.2
16
−4.0
21
−2.2
20
−6.7 17
−4.0
18
0.8
14
0.6
24
1.8 17−7.5
18 −12.914
2.624
9.6 24
−2.67
7.2
0.08
7.9 −9.19
7.7 −4.07 6.8
3.47
9
2.9
10
2.23
6.3
1.92
4.6
3.55
9.2
4.1
14
1.73
4.9
1.25
3
6.1
16
3.3
11
8.4
24
5.6 13
5.2
14
4.8
16
2.90
8.2
7.0 17
4.5
12
2.87
9.8
3.9
11
9.9 240.31
5.6 −2.42
4.4 −1.17
5.9
1.78
5.6
−5.5
10
−2.7 101.5
10
1.1 10
0.1
20
17.2
50
2.0
20
0.3 10> 400 MeV
Multi-GeV
Multiring Sub-GeV
Multiring multi-GeV
Particle identification (1 ring)

SK-II

Sub-GeV
Multi-GeV

Particle identification (multiring) Sub-GeV
Multi-GeV
Energy calibration
Up/down asymmetry energy calibration
UP-μ reduction
Stopping
Through-going
UP-μ stopping/through-going separation
Energy cut for stopping UP-μ
Path length cut for through-going UP-μThrough-going UP-μ showering separation
Background subtraction for UP-μ Stoppinga
Nonshoweringa
Showeringa
νe =ν̄e separation
100 < Pe < 250 MeV=c
Sub-GeV 1-ring π 0 selection
250 < Pe < 400 MeV=c
400 < Pe < 630 MeV=c
630 < Pe < 1000 MeV=c1000 < Pe < 1330 MeV=c
Sub-GeV 2-ring π 0
Decay-e tagging
Solar activity

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

a[1] Y. Fukuda et al. (Super-Kamiokande Collaboration), Phys.
Rev. Lett. 81, 1562 (1998).
[2] Y. Ashie et al. (Super-Kamiokande Collaboration), Phys.
Rev. Lett. 93, 101801 (2004).
[3] K. Abe et al. (Super-Kamiokande Collaboration), Phys.
Rev. D 83, 052010 (2011).
[4] B. T. Cleveland, T. Daily, R. Davis, Jr., J. R. Distel, K. Lande,
C. K. Lee, P. S. Wildenhain, and J. Ullman, Astrophys. J. 496,
505 (1998).
[5] J. N. Abdurashitov et al. (SAGE Collaboration), Phys. Rev.
C 80, 015807 (2009).
[6] M. Altmann et al. (GNO Collaboration), Phys. Lett. B 616,
174 (2005).

[7] W. Hampel et al. (GALLEX Collaboration), Phys. Lett. B
447, 127 (1999).
[8] B. Aharmim et al. (SNO Collaboration), Phys. Rev. C 88,
025501 (2013).
[9] S. Abe et al. (KamLAND Collaboration), Phys. Rev. Lett.
100, 221803 (2008).
[10] F. P. An et al. (DAYA-BAY Collaboration), Phys. Rev. Lett.
108, 171803 (2012).
[11] Y. Abe et al. (Double Chooz Collaboration), Phys. Lett. B
723, 66 (2013).
[12] J. K. Ahn et al. (RENO Collaboration), Phys. Rev. Lett. 108,
191802 (2012).
[13] M. H. Ahn et al. (K2K Collaboration), Phys. Rev. D 74,
072003 (2006).

052019-20

LIMITS ON STERILE NEUTRINO MIXING USING …

PHYSICAL REVIEW D 91, 052019 (2015)[14] P. Adamson et al. (MINOS Collaboration), Phys. Rev. Lett.
110, 251801 (2013).
[15] K. Abe et al. (T2K Collaboration), Phys. Rev. Lett. 112,
181801 (2014).
[16] K. Abe et al. (Super-Kamiokande Collaboration), Phys.
Rev. Lett. 110, 181802 (2013).
[17] K. Abe et al. (T2K Collaboration), Phys. Rev. D 88, 032002(2013).
[18] N. Agafonova et al. (OPERA Collaboration), J. HighEnergy Phys. 11 (2013) 036.
[19] S. Schael et al. (ALEPH Collaboration, DELPHI Collaboration, L3 Collaboration, OPAL Collaboration, SLD Collaboration, LEP Electroweak Working Group, SLD
Electroweak Group, SLD Heavy Flavour Group), Phys.
Rep. 427, 257 (2006).
[20] A. Aguilar et al. (LSND Collaboration), Phys. Rev. D 64,
112007 (2001).
[21] A. A. Aguilar-Arevalo et al. (MiniBooNE Collaboration),
Phys. Rev. Lett. 110, 161801 (2013).
[22] P. Huber, Phys. Rev. C 84, 024617 (2011).
[23] G. Mention, M. Fechner, Th. Lasserre, Th. A. Mueller, D.
Lhuillier, M. Cribier, and A. Letourneau, Phys. Rev. D 83,
073006 (2011).
[24] P. Anselmann et al. (GALLEX Collaboration), Phys. Lett. B
342, 440 (1995).
[25] W. Hampel et al. (GALLEX Collaboration), Phys. Lett. B
420, 114 (1998).
[26] J. N. Abdurashitov, V. N. Gavrin, S. V. Girin, V. V.
Gorbachev, T. V. Ibragimova et al., Phys. Rev. Lett. 77,
4708 (1996).
[27] J. N. Abdurashitov et al. (SAGE Collaboration), Phys. Rev.
C 59, 2246 (1999).
[28] J. N. Abdurashitov, V. N. Gavrin, S. V. Girin, V. V.
Gorbachev, P. P. Gurkina et al., Phys. Rev. C 73, 045805(2006).
[29] J. Kopp, P. A. N. Machado, M. Maltoni, and T. Schwetz, J.
High Energy Phys. 05 (2013) 050.
[30] P. A. R. Ade et al. (Planck Collaboration), Astron. Astrophys. 571, A16 (2014).
[31] S. Das, T. Louis, M. R. Nolta, G. E. Addison,
E. S. Battistelli et al., J. Cosmol. Astropart. Phys. 04(2014) 014.
[32] C. L. Bennett et al. (WMAP), Astrophys. J. Suppl. Ser. 208,
20 (2013).
[33] K. T. Story, C. L. Reichardt, Z. Hou, R. Keisler, K. A. Airdet al., Astrophys. J. 779, 86 (2013).
[34] I. E. Stockdale, A. Bodek, F. Borcherding, N. Giokaris, K.
Lang et al., Phys. Rev. Lett. 52, 1384 (1984).
[35] G. Cheng et al. (MiniBooNE Collaboration, SciBooNE
Collaboration), Phys. Rev. D 86, 052009 (2012).
[36] P. Adamson et al. (MINOS Collaboration), Phys. Rev. Lett.
107, 011802 (2011).
[37] J. M. Conrad, C. M. Ignarra, G. Karagiorgi, M. H. Shaevitz,
and J. Spitz, Adv. High Energy Phys. 2013, 1 (2013).
[38] C. Giunti, Sterile neutrino status, arXiv:1311.1335.
[39] There were 11,129 PMTs during SK-III and SK-IV. During
SK-I there were 11,146 PMTs and during SK-II there were
5,182 PMTs.
[40] Y. Ashie et al. (Super-Kamiokande Collaboration), Phys.
Rev. D 71, 112005 (2005).

[41] R. Wendell et al. (Super-Kamiokande Collaboration), Phys.
Rev. D 81, 092004 (2010).
[42] 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).
[43] M. Ishitsuka, L=E analysis of the atmospheric neutrino datafrom Super-Kamiokande, Ph.D. thesis, University of Tokyo,
2004.
[44] F. Dufour, Precise study of the atmospheric neutrino
oscillation pattern using Super-Kamiokande I and II,
Ph.D. thesis, Boston University, 2009.
[45] J. Beringer et al. (Particle Data Group), Phys. Rev. D 86,
010001 (2012).
[46] S. Desai et al. (Super-Kamiokande Collaboration), Astropart. Phys. 29, 42 (2008).
[47] Y. Hayato, Acta Phys. Pol. B 40, 2477 (2009).
[48] R. Bradford, A. Bodek, H. S. Budd, and J. Arrington, Nucl.
Phys. B, Proc. Suppl. 159, 127 (2006).
[49] K. Kobayashi, H. Akimune, H. Ejiri, H. Fujimura, M.
Fujiwara et al., De-excitation gamma-rays from the s-holestate in N-15 associated with proton decay in O-16, arXiv:
nucl-ex/0604006.
[50] M. Yosoi, Structures and fragmentations of the deephole states in 11B and 15N, Ph.D. thesis, Kyoto University,
2003.
[51] T. Mori, R. Yamaguchi, M. Sakuda, A. M. Ankowski, andO. Benhar, AIP Conf. Proc. 1405, 121 (2011).
[52] T. S. H. Lee and R. P. Redwine, Annu. Rev. Nucl. Part. Sci.
52, 23 (2002).
[53] C.-H. Oh, R. A. Arndt, I. I. Strakovsky, and R. L. Workman,
Phys. Rev. C 56, 635 (1997).
[54] R. Workman, SAID partial wave analysis program.
Center for Nuclear Studies, George Washington University,
2009.
[55] R. A. Arndt, W. J. Briscoe, I. I. Strakovsky, R. L.
Workman, and M. M. Pavan, Phys. Rev. C 69, 035213(2004).
[56] D. Rowntree et al. (LADS Collaboration), Phys. Rev. C 60,
054610 (1999).
[57] B. G. Ritchie, Phys. Rev. C 44, 533 (1991).
[58] M. Honda, T. Kajita, K. Kasahara, and S. Midorikawa, Phys.
Rev. D 83, 123001 (2011).
[59] K. Abe, Y. Hayato, T. Iida, K. Iyogi, J. Kameda et al.,Nucl. Instrum. Methods Phys. Res., Sect. A 737, 253(2014).
[60] M. Maltoni and T. Schwetz, Phys. Rev. D 76, 093005(2007).
[61] B. Pontecorvo, Sov. Phys. JETP 26, 984 (1968).
[62] Z. Maki, M. Nakagawa, and S. Sakata, Prog. Theor. Phys.
28, 870 (1962).
[63] A. M. Dziewonski and D. L. Anderson, Phys. Earth Planet.
Inter. 25, 297 (1981).
[64] V. Barger, K. Whisnant, S. Pakvasa, and R. J. N. Phillips,
Phys. Rev. D 22, 2718 (1980).
[65] G. L. Fogli, E. Lisi, A. Marrone, D. Montanino, and A.
Palazzo, Phys. Rev. D 66, 053010 (2002).
[66] C. H. Llewellyn Smith, Phys. Rep. 3, 261 (1972).
[67] J. Nieves, J. E. Amaro, and M. Valverde, Phys. Rev. C 70,
055503 (2004).

052019-21

K. ABE et al.

PHYSICAL REVIEW D 91, 052019 (2015)[68] See particularly Figs. 12(a) and 12(b) of Appendix C of [60].
[69] Likelihood surfaces are available for download as tables:
http://www‑sk.icrr.u‑tokyo.ac.jp/indico/conferenceDisplay.py?confId=2047.
[70] Y. Declais, J. Favier, A. Metref, H. Pessard, B. Achkar et al.,Nucl. Phys. B434, 503 (1995).

[71] A. Capella, A. Kaidalov, C. Merino, and J. Tran Thanh Van,
Phys. Lett. B 337, 358 (1994).
[72] M. Gluck, E. Reya, and A. Vogt, Eur. Phys. J. C 5, 461(1998).
[73] E. Hernandez, J. Nieves, and M. Valverde, Phys. Rev. D 76,
033005 (2007).

052019-22