PHYSICAL REVIEW D 94, 052001 (2016)

Measurements of the atmospheric neutrino flux by Super-Kamiokande:
Energy spectra, geomagnetic effects, and solar modulation
E. Richard,2 K. Okumura,2,30 K. Abe,1,30 Y. Haga,1 Y. Hayato,1,30 M. Ikeda,1 K. Iyogi,1 J. Kameda,1,30 Y. Kishimoto,1,30
M. Miura,1,30 S. Moriyama,1,30 M. Nakahata,1,30 T. Nakajima,1 Y. Nakano,1 S. Nakayama,1,30 A. Orii,1 H. Sekiya,1,30
M. Shiozawa,1,30 A. Takeda,1,30 H. Tanaka,1 T. Tomura,1,30 R. A. Wendell,1,30 R. Akutsu,2 T. Irvine,2 T. Kajita,2,30
K. Kaneyuki,2,30,* Y. Nishimura,2 L. Labarga,3 P. Fernandez,3 J. Gustafson,4 C. Kachulis,4 E. Kearns,4,30 J. L. Raaf,4
J. L. Stone,4,30 L. R. Sulak,4 S. Berkman,5 C. M. Nantais,5 H. A. Tanaka,5 S. Tobayama,5 M. Goldhaber,6,* W. R. Kropp,7
S. Mine,7 P. Weatherly,7 M. B. Smy,7,30 H. W. Sobel,7,30 V. Takhistov,7 K. S. Ganezer,8 B. L. Hartfiel,8 J. Hill,8 N. Hong,9J. Y. Kim,9 I. T. Lim,9 R. G. Park,9 A. Himmel,10 Z. Li,10 E. O’Sullivan,10 K. Scholberg,10,30 C. W. Walter,10,30
T. Wongjirad,10 T. Ishizuka,11 S. Tasaka,12 J. S. Jang,13 J. G. Learned,14 S. Matsuno,14 S. N. Smith,14 M. Friend,15
T. Hasegawa,15 T. Ishida,15 T. Ishii,15 T. Kobayashi,15 T. Nakadaira,15 K. Nakamura,15,30 Y. Oyama,15 K. Sakashita,15
T. Sekiguchi,15 T. Tsukamoto,15 A. T. Suzuki,16 Y. Takeuchi,16,30 T. Yano,16 S. V. Cao,17 T. Hiraki,17 S. Hirota,17 K. Huang,17T. Kikawa,17 A. Minamino,17 T. Nakaya,17,30 K. Suzuki,17 Y. Fukuda,18 K. Choi,19 Y. Itow,19 T. Suzuki,19 P. Mijakowski,20
K. Frankiewicz,20 J. Hignight,21 J. Imber,21 C. K. Jung,21 X. Li,21 J. L. Palomino,21 M. J. Wilking,21 C. Yanagisawa,21,†
D. Fukuda,22 H. Ishino,22 T. Kayano,22 A. Kibayashi,22 Y. Koshio,22 T. Mori,22 M. Sakuda,22 C. Xu,22 Y. Kuno,23
R. Tacik,24,32 S. B. Kim,25 H. Okazawa,26 Y. Choi,27 K. Nishijima,28 M. Koshiba,29 Y. Totsuka,29,* Y. Suda,35
M. Yokoyama,35,30 C. Bronner,30 M. Hartz,30 K. Martens,30 Ll. Marti,30 Y. Suzuki,30 M. R. Vagins,30,7 J. F. Martin,31
A. Konaka,32 S. Chen,33 Y. Zhang,33 and R. J. Wilkes34(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, Japan
20
National Centre For Nuclear Research, 00-681 Warsaw, Poland
21
Department of Physics and Astronomy, State University of New York at Stony Brook,
New York 11794-3800, USA
22
Department of Physics, Okayama University, Okayama, Okayama 700-8530, Japan
23
Department of Physics, Osaka University, Toyonaka, Osaka 560-0043, Japan
24
Department of Physics, University of Regina, 3737 Wascana Parkway,
Regina, Saskatchewan S4SOA2, Canada
25
Department of Physics, Seoul National University, Seoul 151-742, Korea
26
Department of Informatics in Social Welfare, Shizuoka University of Welfare,
Yaizu, Shizuoka 425-8611, Japan
27
Department of Physics, Sungkyunkwan University, Suwon 440-746, Korea
28
Department of Physics, Tokai University, Hiratsuka, Kanagawa 259-1292, Japan
29
The University of Tokyo, Bunkyo, Tokyo 113-0033, Japan

2470-0010=2016=94(5)=052001(31)052001-1© 2016 American Physical Society

E. RICHARD et al.

PHYSICAL REVIEW D 94, 052001 (2016)30

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-8582, Japan
31
Department of Physics, University of Toronto, 60 St. George Street, Toronto, Ontario M5S1A7, Canada
32
TRIUMF, 4004 Wesbrook Mall, Vancouver, British Columbia V6T2A3, Canada
33
Department of Engineering Physics, Tsinghua University, Beijing 100084, China
34
Department of Physics, University of Washington, Seattle, Washington 98195-1560, USA
35
Department of Physics, University of Tokyo, Bunkyo, Tokyo 113-0033, Japan(Received 29 October 2015; published 2 September 2016)
A comprehensive study of the atmospheric neutrino flux in the energy region from sub-GeV up to several
TeV using the Super-Kamiokande (SK) water Cherenkov detector is presented in this paper. The energy andazimuthal spectra, and variation over time, of the atmospheric νe þ ν̄e and νμ þ ν̄μ fluxes are measured. Theenergy spectra are obtained using an iterative unfolding method by combining various event topologies
with differing energy responses. The azimuthal spectra depending on energy and zenith angle, and theirmodulation by geomagnetic effects, are also studied. A predicted east-west asymmetry is observed in boththe νe and νμ samples at 8.0σ and 6.0σ significance, respectively, and an indication that the asymmetry
dipole angle changes depending on the zenith angle was seen at the 2.2σ level. The measured energy andazimuthal spectra are consistent with the current flux models within the estimated systematic uncertainties.
A study of the long-term correlation between the atmospheric neutrino flux and the solar magnetic activity
cycle is performed, and a weak preference for a correlation was seen at the 1.1σ level, using SK-I–SK-IV
data spanning a 20-year period. For several particularly strong solar activity periods, corresponding to
Forbush decrease events, no theoretical prediction is available but a deviation below the typical neutrino
event rate is seen at the 2.4σ level. The seasonal modulation of the neutrino flux is also examined, but thechange in flux at the SK site is predicted to be negligible, and, as expected, no evidence for a seasonalcorrelation is seen.
DOI: 10.1103/PhysRevD.94.052001

I. INTRODUCTION
Atmospheric neutrinos are one of the main experimentally available neutrino sources, observed in a wide energy
region from 100 MeV up to the PeV scale [1]. They are
generated after the interaction of cosmic rays with the air
molecules in Earth’s atmosphere, from the decay of
secondary particles such as π and K mesons. The formof the energy spectrum is well approximated by a powerlaw, although it is relatively suppressed below the GeV
scale due to the rigidity cutoff effect on the primary cosmic
rays, caused by Earth’s magnetic field. The flavor ratio ofνμ þ ν̄μ to νe þ ν̄e below the GeV scale is approximately 2,
considering the dominance of the π  decay chains, but
increases towards higher energies. In neutrino detectors,
atmospheric neutrinos are observed coming from all
directions, as the Earth is mostly transparent for neutrinos
below the PeV scale; the flight length from the production
point to the detector thus varies from Oð10Þ to Oð104 Þ km,
depending on the zenith angle of the arrival direction. The
neutrino flux in the horizontal direction is generally higher
than that in the vertical direction, due to the longer path ofthe parent particles in the atmosphere; however, in the subGeV region there is an east-west asymmetry due to theazimuthal dependence in the rigidity cutoff energy of the*

Deceased.
Also at Science Department, BMCC/CUNY, New York,New York, USA.
†

primary cosmic rays. Towards the higher energies, fewer
neutrinos are produced as the π and K decay lengths
become longer than their path lengths in the atmosphereand the parent particles reach the ground before decaying.
Above approximately 100 TeV, so-called “prompt” neutrinos coming from the fast decay of charmed mesons are
expected to dominate, due to their much shorter decay
length.
Since the first detection of the atmospheric neutrino inunderground experiments in the 1960s [2,3], further measurements brought the discovery of neutrino oscillation in1998 [4]. The continuing series of independent neutrino
oscillation measurements in solar- [5,6], atmospheric[7,8], reactor- [9–12], and accelerator-sourced [13–16]
neutrinos are consistent with three neutrinos mixing as
described by the 3 × 3 Pontecorvo-Maki-Nakagawa-Sakata(PMNS) matrix [17,18]; this is generally accepted as thestandard framework of neutrino oscillation, while theCP-violation phase and the mass ordering between thesecond and third mass states are not yet known.
The study of the atmospheric neutrino is generally basedon predictions of the expected flux, coming fromMonte Carlo (MC) simulations. In this paper we will
discuss three such simulations performed by the HKKM
[19,20], Bartol [21], and FLUKA [22] research groups. The
HKKM model defines a primary cosmic-ray spectrum
based on measurements from the BESS [23,24] andAMS [25,26] Collaborations. For interactions in theatmosphere, the DPMJET-III [27] hadronic interaction

052001-2

MEASUREMENTS OF THE ATMOSPHERIC NEUTRINO FLUX …

model is used, with some customized tuning for better
agreement with the cosmic ray muon data. The JAM
nuclear interaction model [28], which has better agreementwith the π production measurements by the HARP experiment [29], is also introduced below 32 GeV in the more
recent HKKM11 model [20], resulting in a relative increaseof the neutrino flux below 1 GeV when compared to
HKKM07 [19]. The Bartol model adopts a primary proton
spectrum that is relatively lower below 50 GeV, and thehigh energy nucleon flux has a flatter energy dependence[30]. TARGET version 2.1 [31] is used for the decay andinteraction generator. The FLUKA flux model is so named
because it uses the FLUKA Monte Carlo simulation code[32], a widely used hadronic and electromagnetic interaction model. In all of these simulations, three-dimensional
particle tracking is performed for the primary and secondary particles.
The estimated uncertainties on the atmospheric neutrino
flux predictions are currently between ∼5%–25% in the
100 MeV → 100 GeV range. As the neutrino energy
increases above 10 GeV, the uncertainties in the π and K
fluxes become the dominant error sources [33,34]. The
neutrino flux predictions are consistent within ∼20% below
32 GeV [20], with larger differences in the higher-energy
region [19].
The current generation of neutrino detection experiments
have increased their statistics and reduced systematic errors,such that direct measurements can now be made withuncertainties comparable to those on the flux models. In thispaper, we discuss several types of direct measurements madeby the Super-Kamiokande (SK) detector. First we briefly
discuss these measurements in a conceptual sense, especiallytaking into account the fact that the Super-Kamiokande
atmospheric data was previously used to make an original
measurement of the neutrino oscillation parameters.
In a general sense, the measured event rate N of an
atmospheric neutrino detector is expressed by the convolution of atmospheric neutrino flux Φ, neutrino oscillation probability O, neutrino cross section σ, and detector
efficiency ϵ,
N ¼ Φ ⊗ O ⊗ σ ⊗ ϵ:

ð1:1ÞTo measure the oscillation probability O, the quantities Φ, σ,and ϵ must be determined by other independent measurements in advance. As described above, Φ is calculated byMonte Carlo simulations, σ is modeled based on the resultsof neutrino beam experiments, and ϵ is determined based on
detector calibrations and simulations. Combining these
inputs and the measured event rate N, the oscillation
probability O and, thus, the oscillation parameters havebeen previously measured by Super-Kamiokande, such as
the mass difference jΔm232 j. Conversely, it is possible tomeasure the atmospheric neutrino flux Φ if the neutrino
oscillation probability O is given. The measurements of the

PHYSICAL REVIEW D 94, 052001 (2016)

oscillation parameters by the Super-Kamiokande atmospheric data have since been confirmed by several independent experiments, and a consistent and accurate three-flavor
oscillation parameterization has emerged by combining datafrom each of the independent experiments. Thus, by taking
the Particle Data Group values of the oscillation parameters[35] as the input, the neutrino flux Φ can be measured at
Super-Kamiokande.
We separate our measurements of the atmospheric
neutrino flux into three main categories, as follows.
The first measurement is the energy spectra of the νeand νμ fluxes. Over the past decade, the simulations have
improved their statistics and calculation methods and now
provide accurate predictions across a wide energy range
from sub-GeV to 10 TeV [20–22]. These predictions areimportant for studies of the atmospheric neutrino itself as
well as serving as a background model for rare event
searches such as proton decay [36,37] or dark matter [38].
An initial measurement of the energy spectra was madeby the Frejus experiment in 1995 [39], before the existenceof neutrino oscillation was known. More recently,
measurements were made above TeV energies by the
cubic-kilometer detector IceCube [40–45]. At theseenergies, the understanding of the atmospheric flux isimportant with respect to searches of astronomical neutrinos, a flux of which was recently discovered by the
IceCube Collaboration [46–48].
A precise measurement of atmospheric neutrino flux
below 100 GeV has not yet been published, and the SuperKamiokande detector is able to make several significant
improvements with respect to the Frejus measurement inthis region. The larger detector size increases event statistics by around 2 orders of magnitude and extends themeasurement down to 100 MeV and up to around 10 TeV,
which overlaps with the low-energy end of the cubickilometer detector measurements. The high resolution of
the Cherenkov ring imaging technique used in SuperKamiokande leads to excellent identification of the neutrino flavor and background elimination. Combined withthe progress in understanding neutrino cross sections (see,e.g., Sec. 49 of [35] and references therein) and oscillations
by several independent experiments, an accurate measurement of the flux based on the number of observed events inthe detector is possible.
The second measurement is the azimuthal distribution of
the atmospheric neutrino flux. The geomagnetic field
deflects incoming primary cosmic rays, depending on their
momentum and nuclear composition; the geomagnetic
rigidity is defined as pc=Ze for momentum p, atomic
number Z, and the elementary charge e, and is often given
in units of GV. For a given arrival location and direction on
the Earth’s surface, only cosmic rays with rigidity above acertain threshold will have been able to traverse thegeomagnetic field to that point, excluding trajectories that
had intersected the Earth; this threshold is known as the052001-3

E. RICHARD et al.

PHYSICAL REVIEW D 94, 052001 (2016)

geomagnetic rigidity cutoff. The structure of the geomagnetic field causes variations in this cutoff, the primaryeffect of which is an east-west oriented anisotropy incosmic rays at the Earth’s surface; this was originallydetected in the 1930s [49] and used to infer that cosmic rays
are generally positively charged.
A similar east-west anisotropy is expected in the atmospheric neutrino flux; in a previous Super-Kamiokande
measurement [50], it was discovered in the νe sample(5σ) but seen only with low significance in the νμ sample(2σ). Since that measurement, the flux simulations have also
progressed from simple one-dimensional to three-dimensional calculations, and began to use complex geomagnetic
field models instead of dipole approximations. The HKKM
and Bartol simulations have further included bending of
secondary particles due to the geomagnetic field in theatmosphere. These changes have led to significant modifications in the predictions of the azimuthal anisotropies[51–53]. Measuring in detail the angular distributions, inaddition to the previously well-studied zenith distributions,
can thus be used as a further test of the flux simulations andtheir implementations of geomagnetic effects, and also to
confirm the discovery of such azimuthal anisotropies.
The third measurement is the modulation of the neutrino
flux over time. The solar cycle is an oscillatory change inthe solar activity, such as the level of plasma emissions,with an average period of approximately 11 yr. The cosmicray flux at Earth is well known to be anticorrelated with thesolar activity [54]. This is essentially because the plasma
flux (or “solar wind”) from the sun can scatter cosmic rays
entering the solar system, and, therefore, during periods ofhigh solar activity the cosmic-ray flux is relatively reduced.
Consequently, the atmospheric neutrino flux is predicted toalso be anticorrelated with the solar cycle, although this has
not previously been measured.
Historically, the solar activity was measured by its
correlation with the appearance of sunspots. However,since 1948, the use of neutron monitors (NMs) has
provided a method to accurately and continuously trackthe neutron flux at the Earth’s surface resulting from cosmic
ray impacts [55]; the NM counts are generally believed to
be well correlated with the primary cosmic-ray activity. In
this paper we test for the anticorrelation of the atmospheric
neutrino flux with the solar cycle by searching for acorrelation between the neutrino flux at SK and the neutron
detection rates at various NMs operated by other institutes.
This method is able to test on short time scales of Oð1 hÞ,as both neutrino and neutron observations take place on
Earth, and the propagation of the solar wind within the solar
system (with speed of the order 100 km s−1 ) need not be
considered. We assume that the effect of the solar wind is
uniform in the neighborhood of the Earth, such that we may
expect a good correlation between the neutrino flux at theSK site and the neutron flux monitored at the NMs invarious locations around the Earth.

Yearly changes in the atmospheric neutrino flux are also
expected, due to seasonal temperature variations. In thesummer months the atmospheric density is increased at
higher altitudes, and relatively more neutrinos are createdby secondary particles decaying in flight. Such changes are
predicted [56] to be strongest at the polar regions, with anormalization change of a few percent around the GeV to
TeV scale, but to become negligible moving towards to theequator where the seasonal variation in air density is less.
While the variation of the atmospheric neutrino flux at SK
is thus expected to be minimal, we also test for such a
correlation in this paper.
This paper continues by explaining the detector andatmospheric neutrino data set in Sec. II, the energy
spectrum measurement in Sec. III, the azimuthal spectrum
measurement in Sec. IV, the solar modulation measurementin Sec. V, and concludes with a summary in Sec. VI.
II. DETECTOR AND ATMOSPHERIC
NEUTRINO DATA SET
A. Super-Kamiokande detector
Super-Kamiokande (SK) is a water Cherenkov detector
located in the Mozumi mine of Gifu prefecture, Japan [57],at geographic coordinates 36°250 32.600 N 137°180 37.100 E
and altitude 370 m in the WGS-84 system [58]. About
11 000 20-in. photomultipliers (PMTs) are mounted on thewall of the detector facing inwards, and observe theCherenkov light emitted by charged particles producedby neutrino interactions in the ultrapure water. An opticallyseparated region on the outer side also contains about 1800
PMTs, which act as a veto for incoming cosmic ray events.
The detector has excellent particle identification (PID)
capability by using the Cherenkov ring’s pattern andopening angle, separating showering-type events fromtrack-type events, which are denoted as e-like or μ-likeevents, respectively.
SK has had four experimental phases so far. These aredesignated as SK-I (1996–2001), SK-II (2002–2005), SKIII (2006–2008), and SK-IV (2008–). The major detector
changes that distinguish these periods are as follows. The
SK-I period was ended when an implosion accident
destroyed about half of the PMTs; the remaining PMTs
were rearranged for even coverage during SK-II. For
SK-III, full PMT coverage was restored, and with SK-IV
came an improved front-end electronics system [59]. SKIV is an ongoing experiment; the data sets used in this paper
include data until September 2014 for the energy spectrum
and azimuthal analyses described in Secs. III and IV, anduntil April 2015 for the solar modulation analysis in Sec. V.
B. Atmospheric data types
The atmospheric neutrino data events are separated into
three main samples with different event topologies. In the
fully contained (FC) sample, the reconstructed neutrino

052001-4

MEASUREMENTS OF THE ATMOSPHERIC NEUTRINO FLUX …

interaction vertex is inside the 22.5-kton fiducial volume,
which is the inner region with a boundary 2 m inside theinner wall, and all visible secondary particles are contained
inside the inner detector. In the partially contained (PC)sample, the vertex is also inside this fiducial volume, but
outgoing particles are allowed to exit the inner detector. PC
events typically have longer charged particle tracks, whichare, therefore, mostly muons induced by νμ charged-current(CC) interactions. In the upwardgoing-muon (UPMU)
event sample, neutrinos interacting with the surroundingrock create muons which enter the detector from below(downgoing muons are ignored, as these are overwhelmingly produced by cosmic rays). The UPMU sample is alsoa predominantly νμ -induced sample. The efficiency inselecting true neutrino events is estimated for fully contained events as > 97%, and for partially contained
events as >80% (improving to > 95% for SK-III and
SK-IV) [60]. Non-neutrino background contaminationfrom cosmic-ray muons or light-flashing PMTs is lessthan a few percent for all samples.
In this analysis, these three main event samples arefurther divided as follows, according to their detailed
properties identified by the event reconstruction algorithm,
such as PID, number of reconstructed Cherenkov rings,
reconstructed energy, and presence or absence of electrons
from delayed muon decays. These subsample definitions
are similar to the standard ones used in the SK atmospheric
neutrino oscillation analysis [8], although compared to thatanalysis some subsamples are combined or discarded inthis paper. The detector cannot distinguish neutrinos fromantineutrinos by measuring the sign of the produced lepton,
as no magnetic field is applied in the detector; therefore,
when we refer to νμ or νe samples this includes also ν̄μ andν̄e , respectively (some separation techniques in fact exist[61–63] that take advantage of various kinematicdifferences, but they are not used in this paper).
Fully contained events are divided into e-like and μ-likesamples according to the PID of the most energetic Cherenkov
ring, corresponding generally to νe and νμ samples, respectively. They are further divided according to whether thenumber of the identified Cherenkov rings is one (single-ring)or more (multiring). Events are also classified into “sub-GeV”
or “multi-GeV” samples at a threshold of 1.33 GeV in total
visible energy, which is the sum of the visible energy among
all fitted Cherenkov rings. Single-ring events are of higher
purity in neutrino flavor, while multiring events tend to havehigher neutrino energy. Additional selection criteria to reducethe background of neutral-current (NC) and wrong-flavor CC
events are applied; for the sub-GeV single-ring e-like events,
NC-produced π 0 backgrounds are reduced by a π 0 finderalgorithm [64] which identifies events using a likelihood
method, based around the reconstructed invariant mass
assuming a π 0 → 2γ decay event (where the second γ ring
was originally missed due to low energy or ring overlap).
Sub-GeV multiring e-like events are also not used in these

PHYSICAL REVIEW D 94, 052001 (2016)

analyses in order to further avoid NC π 0 background. Anotherlikelihood-based cut is applied for multi-GeV multiring e-likeevents, which is the same method as described in [8] but
extended to include SK-IV, to reduce background contamination from hadronic events dominated by γ rings producedby π 0 created in CC νμ and NC-induced events.
Partially contained events are divided into two samples,“PC stopping” and “PC throughgoing,” based on whether
the muon particle stops in the outer detector, or leaves thedetector completely. The separation is made according to
the charge sum in the outer detector around the muon
exiting point, considering the detector geometry [65].
Upwardgoing-muon events are divided into threesamples: “UPMU stopping,” “UPMU nonshowering,”and “UPMU showering.” UPMU stopping events have
only an entering hit cluster in the outer detector, while theother two samples have both entering and exiting hit
clusters. UPMU nonshowering and showering events are
separated by the estimated deposited energy per unit muon
range, with the aim to separate muon events with radiativeemission, such as bremsstrahlung, which dominates whenthe muon’s energy is greater than around 1 TeV [66].
C. Data and Monte Carlo sets
Table I summarizes the detector live times and numbersof events in the data and MC sets used by this analysis,
where the MC simulations are generated corresponding to
500 yr of live time for each SK period. Neutrino interactions are generated from the HKKM07 [19] flux andNEUT [67] neutrino interaction models, and then passed toour detector simulation. The recent updates to NEUT and
the detector simulation are described in [68]. The data
reduction and event reconstruction processes, which are thestandard ones used for the SK atmospheric neutrino
analysis [8], are applied to both the data and MC samples.
Monte Carlo correction weights from the HKKM07
model to the HKKM11 [20] model are applied on an
event-by-event basis, considering the differences in theenergy, azimuthal, and zenith angle distributions. In somecases in Sec. III, weights will instead be applied to reweight
other flux models, according to their predicted energy
spectra. Weights due to neutrino oscillation are also applied
per event, using the three-flavor parameterization in [35],under the assumptions of normal hierarchy (Δm232 > 0) andCP symmetry (δCP ¼ 0). Finally, weights due to changes in
the average degree of solar activity were estimated based on
cosmic NM data [69] and applied depending on the SK
period, except for the solar modulation analysis in Sec. V,
where a more accurate per-event solar correction is
described and used.
The fractions of the interaction mode for each subsample(i.e., νe CC, νμ CC, or NC) are estimated from the MC data
sets, and are also shown in Table I. We see that all e-likeand μ-like subsamples have a high purity of νe and νμ CC
interactions, respectively. The systematic uncertainties on

052001-5

E. RICHARD et al.

PHYSICAL REVIEW D 94, 052001 (2016)

TABLE I. Live times and numbers of events for each SK period and each subsample. Numbers in brackets are the MC expectation
based on the HKKM11 [20] flux model, which are calculated separately for each SK period and then scaled to the data live time. The
oscillation and solar activity weights are included. UPMU data live times are longer due to less strict conditions for good run selection,and UPMU event numbers are fractional due to the subtraction of horizontally going cosmic muon backgrounds. The fractions of theneutrino interaction modes for each subsample, estimated from the MC data sets averaging over all SK periods, are shown in therightmost four columns.
Subsample

SK-I

SK-II

SK-IIIFC and PC
UPMU

1489
1646

799
518
828
636
Number of events

SK-IV

Total

Live time (days)
1993
1993

4799
5103
Interaction (%)νe CC νμ CC ντ CC NC

FC e-like
Sub-GeV single-ring 3288 (3104.7) 1745 (1632.8) 1209 (1100.7) 4251 (4072.8) 10493 (9911.0)
Multi-GeV single-ring 856 (842.8)
396 (443.7) 274 (299.5) 1060 (1080.0) 2586 (2666.0)
Multi-GeV multiring
449 (470.1)
267 (252.1) 140 (161.9)
634 (654.9)
1490 (1539.0)
FC μ-like
Sub-GeV single-ring 3184 (3235.6) 1684 (1731.8) 1139 (1152.0) 4379 (4394.7) 10386 (10514.0)
Multi-GeV single-ring 712 (795.4)
400 (423.9) 238 (273.9)989(1051.5)
2339 (2544.7)
Multi-GeV multiring
603 (656.5)
337 (343.8) 228 (237.9)
863 (927.8)2031 (2166.0)
PC
Stop
143 (145.3)
77 (73.2)
54 (53.3)
237 (229.0)
511 (500.8)
Through
759 (783.8)
350 (383.0) 290 (308.8) 1093 (1146.7) 2492 (2622.3)
UPMU
Stop
432.0 (433.7) 206.4 (215.7) 193.7 (168.3) 492.7 (504.1) 1324.8 (1321.8)
Nonshowering
1564.4 (1352.4) 726.3 (697.5) 612.9 (504.1) 1960.7 (1690.3) 4864.3 (4244.4)
Showering271.7 (291.6) 110.1 (107.0) 110.0 (126.0) 350.1 (274.4)
841.9 (799.0)

these purity estimations were calculated during this study,
and found to be relatively small; less than 4 and 1.5
percentage points for all e-like and μ-like subsamples,
respectively. Figure 1 shows the MC true neutrino energy
distributions of the subsamples, showing that each sample
has a different energy sensitivity. The FC μ-like sample
extends above 10 GeV, while the e-like sample extends up
to 100 GeV, indicating that μ-like events in the fiducial
volume beyond 10 GeV tend to exit the inner detector andbe classified as PC events. The FC e-like sample has alower neutrino energy threshold than the μ-like sample(100 MeV as opposed to 200 MeV), due to the difference in
Cherenkov threshold between the electrons/positrons andthe muons which are produced in the CC interactions inwater. The PC and UPMU samples approximately range
from 1 GeV to 10 TeV, and classification into the respective
subsamples is also seen to cover different energy regions.
III. ENERGY SPECTRUM ANALYSIS
A. Event classification
In order to measure the flux of neutrinos as a functionof their energy, we select all events shown in Table I. This
gives us a generally high purity selection of data induced byeither νe or νμ CC interactions, across a wide energy range.
Noting that PC and UPMU samples generally contain νμinteractions, and that the FC sample may be well separated
by its νe and νμ components, we define two data samples

94.1
86.3
73.0

1.5
3.2
7.6

<0.1 4.4
1.7 8.8
3.3 16.1

0.9
0.4
3.4

94.2
99.1
90.5

<0.1 4.9
0.3 0.2
0.6 5.5

12.7
0.8

81.7
98.2

0.9
0.4

1.0
0.2
0.1

97.7
99.4
99.8

1.0 0.3
0.3 0.1<0.1 0.1

4.6
0.5

for the energy spectrum measurement: a νe sample containing FC e-like events and a νμ sample from FC μ-like,PC, and UPMU events. For the νe (νμ ) sample, the signal
event is defined as a νe (νμ ) CC interaction, and thebackgrounds are CC interactions of other flavor neutrinos
and NC interactions of all flavor neutrinos.
In each subsample, the linear sum of the reconstructed
momentum of each visible particle (or Cherenkov ring) isused as a estimator of the neutrino energy, and denoted as thereconstructed energy Erec. The reconstructed momentum ofeach particle is itself calculated using the observed charge
pattern of the associated ring and the particle’s PID. Theenergy binning is defined as shown in Table II, in log10 Erec .
These bin widths are determined considering the energy
resolution of each subsample, which is typically ∼0.2 in
Δlog10 ðEν =GeVÞ for FC, but ∼0.2–0.3 for PC, ∼0.4 for
UPMU stopping, and ∼0.7 for other UPMU samples, since
the energy deposited outside the detector cannot be observed.
Thus, no energy binning is performed for UPMU showering
or nonshowering events. We refer hereafter to the number of
events in each data bin as Mj for bin index j. Figure 2 shows
the values of Mj for SK-I–SK-IV data.
B. Flux unfolding
To measure the energy spectrum, we employ an
unfolding method. This is a class of method in which atrue spectrum is deconvolved from an experimentally

052001-6

MEASUREMENTS OF THE ATMOSPHERIC NEUTRINO FLUX …
1600
1400
1200
1000
800
600
400
200

multi-GeV single-ring e-like
multi-ring e-like

1600
1400
1200
1000
800
600
400
200

ν μ sample

ν e sample

4000

sub-GeV single-ring e-like

Number of Events

SK I-IV dataν e CCν μ CCν τ CC
NC

3000
2000
1000

multi-GeV single-ring μ-likePC stopping
PC through-going200

UPμ stop

PC thru

PC stop

MR μ

MG 1R μSG 1R μMR e

jFIG. 2. Number of events for each reconstructed energy bin jfor each sample. Colored regions represent νe CC, νμ CC, ντ CC,
and NC contributions from the MC simulations, and red crosses
represent the data values for the entire SK-I–SK-IV period.
Monte Carlo distributions are normalized by the data live times.
Only statistical uncertainties are shown in the error bars of thedata points.

300
250

MG 1R e0multi-ring μ -likeUPμ non-sh
UP μ shower

sub-GeV single-ring μ -like

SG 1R-eNumber of Events / 0.1 Log (E ) / 500 yrs
10ν

PHYSICAL REVIEW D 94, 052001 (2016)

5000

150
100
50

300
250

UPMU stopping
UPMU non-showering
UPMU showering

200
150
100
50
0
-1
10

1

10

2

10

3

10

4

10

Neutrino Energy [GeV]

FIG. 1. Monte Carlo neutrino true-energy distributions for eachsubsample, corresponding to the total expected number of events
across the SK-I–SK-IV periods accounting for neutrino oscillation. From top, the FC e-like, FC μ-like, PC, and UPMU
samples are shown.

measured one, based on a knowledge of the experimental
reconstruction. This differs from forward-fitting methods,
where the experimentally measured spectrum is not deconvolved, but compared directly to model predictions thathave been passed through a simulated experimental
reconstruction; for example, in the case of measuring theenergy spectrum, a simple model could be defined such as a
power-law spectrum with normalization and spectral index
parameters. The reconstructed energy spectra depending
these parameters would then be predicted, and the values ofthe parameters measured by finding the best fit to theexperimentally measured spectra. The benefit of unfolding
methods, on the other hand, is that they allow direct

TABLE II. Binning definitions for the data bins M j (by the reconstructed energy Erec ) for the 11 subsamples.
There are 34 bins in total, among which 12 bins are assigned to the νe sample and 22 bins are assigned to the νμsample.
Subsample
FC e-like
Sub-GeV single-ring (SG 1R e)
Multi-GeV single-ring (MG 1R e)
Multi-GeV multi-ring (MR e)
FC μ-like
Sub-GeV single-ring (SG 1R μ)
Multi-GeV single-ring (MG 1R μ)
Multi-ring (MR μ)PC
Stopping (PC stop)
Throughgoing (PC through)
UPMU
Stopping (UP μ stop)
Nonshowering (UP μ nonshowering)
Showering (UP μ showering)

No. of binsj

Bin edges in log10 ðErec =GeVÞ

5
4
3

1–5
6–9
10–12

−1.0; −0.6; −0.4; −0.2, 0.0, 0.2
0.0,0.4,0.7,1.0,3.0
0.0,0.4,0.7,3.0

5
2
4

13–17
18, 1920–23−0.8; −0.6; −0.4; −0.2, 0.0, 0.2
0.0, 0.4, 2.0−1.0, 0.0, 0.4, 0.7, 2.0

2424, 2526–29−1.0, 0.4, 2.0−1.0, 0.0, 0.4, 0.7, 3.0

3
1
1

30–32
33
34

0.0, 0.4, 0.7, 3.0052001-7

E. RICHARD et al.

PHYSICAL REVIEW D 94, 052001 (2016)

comparisons of the unfolded spectra between experiments,
without restricting the measurement to a particular choiceof parameterization.
In our case, the reconstructed energy spectrum Mj will
be unfolded into a true neutrino energy spectrum, ofcorrect-flavor CC interactions only, which we denote
N CC(the binning of this unfolded spectra, and the coni
version to the actual neutrino flux values, will be explainedlater in this section). In a general sense, the relationshipbetween the true and reconstructed spectra is expressed by
the detector response matrix Aji as
X
Aji N CCMj ¼
i

UPμ stop

1

Aji
0.8

PC thru
PC stop
MR μ

0.6

MG 1R μSG 1R μ

0.4

ð3:1Þ
MR e

i

where Aji can be estimated by the detector MC simulations,and accounts for the inability to reconstruct perfectly thetrue neutrino energy. We can write the inverse relationshipusing the unfolding matrix Uij as
N CC
i ¼

UPμ shower
UPμ non-shjX
U ij M j :However, taking Uij as A−1
ji is in principle a poor approach,as the response matrix will not be estimated perfectly (andmay not even be invertible). Among the several moreadvanced algorithms available, we adopt an iterative
Bayesian method [70,71], using the ROOUNFOLD library[72] for the practical implementation. The method is known
as Bayesian due to its use of Bayes’s theorem in theconstruction of the unfolding matrix at each iterative step.
We chose this library as it is known to be reliable and easy
to implement.
The specific details of the iterative unfolding procedure
are as follows [70]. The first estimation of the detector
response matrix Aji is made using the SK MC data set,using correct-flavor CC events reconstructed into the final
samples. The contributions of background events (wrongflavor or NC events) to each reconstructed bin j are also
recorded, and will be accounted for in the final unfolded
spectrum.
We also apply the normalization condition
PA¼
1,
which means that at this stage detector ineffij jiciencies are not accounted for in the matrix. Figure 3 showsthe estimated response matrix for the SK-IV period, notincluding the background events.
To construct the unfolding matrix, we first define thenotation PðjjiÞ as the probability for an event in true energy
bin i to be detected in the reconstructed energy bin j, wherethese probabilities are exactly equivalent to the values of
the response matrix Aji . We also define a set of priorprobabilities P0 ðiÞ for a single event to fall in true energy
bin i as
ð3:3Þ

0.2

SG 1R-e

νμ

νe

ð3:2Þ

j

N CC
P0 ðiÞ ¼ P MC;i;CC
k N MC;k

MG 1R ei0FIG. 3. Response matrix Aji for the SK-IV period, constructedusing MC events. The vertical axis represents the binning by
subsample and reconstructed energy as defined in Table II, andthe horizontal axis represents the νe and νμ sample binnings
defined later in Table IV. The content of the matrix is normalized
so that the sum of each column becomes unity.

where N CC
MC;i is the default event spectrum taken from theMC simulation, and the denominator ensures the probabilities sum to unity. We can then state Bayes’s theoremPðjjiÞP0 ðiÞ;
P0 ðjÞ

ð3:4Þ

XPðjjiÞP0 ðiÞ:

ð3:5ÞPðijjÞ ¼
where
P0 ðjÞ ¼

i

Now that we have the estimated inverse probabilitiesPðijjÞ, we can use our data Mj to make a first estimationof the true number of events as
N̂ CC
i ¼

X
PðijjÞMj :

ð3:6Þ

j

Considering that this estimation takes inputs from both MC
and real data, it generally lies somewhere between thedefault MC spectrum and the true values. Thus, it is usefulto proceed in an iterative way, by using the normalized N̂ CCi
as new set of prior probabilities P0 ðiÞ to generate newPðijjÞ, and in turn iteratively update N̂ CC
i . The finalCC
iteration is denoted N i , and the operation of the above

052001-8

MEASUREMENTS OF THE ATMOSPHERIC NEUTRINO FLUX …

procedure on the data essentially takes the place of theunfolding matrix U ij in Eq. (3.2).
This unfolding method is seen to depend on the numberof iterations; a low number of iterations may be too close to
the statistically smooth prior values and not fully reflect the
information input by the data, while a high number of
iterations will lead statistical fluctuations in the data to
distort the shape of the unfolded spectrum. In general any
prior will thus cause some small bias in the unfolded
spectrum which cannot be perfectly corrected for, but
importantly can be accurately estimated and included asa systematic error; this will be shown in Sec. III C. In
practice, according to tests with our MC data, the iterative
Bayesian method is strongly data driven and produces goodresults compared to other unfolding methods after only a
few iterations. We set the number of iterations based on a
study using several sets of toy MC data generated from theHKKM spectra, by modifying the normalization andspectral index, and testing their reproducibility dependingon the number of iterations. Five iterations was chosenbecause the unfolded spectra were generally well reproduced and stable.
We next discuss our choice of binning for N CC
i . Theenergy range of the unfolded atmospheric neutrino flux
spectrum is defined for the νe flux as [−0.8, 2.0] in units oflog10 ðEν =GeVÞ, and divided into 11 bins of variable size;for the νμ flux, the range is [−0.6, 4.0], which is divided
into 12 bins. The binning definition is shown later, together
with the measurement results, in Table IV.
These energy ranges and bin widths are determinedconsidering the neutrino true energy coverage and resolution of the flux data samples. Smaller bins of width 0.2 inlog10 ðEν =GeVÞ are adopted below 10 GeV due to therelatively finer energy resolution of the FC sample, andwider bins are adopted for higher energies due to thedeterioration of the energy resolution in the PC and UPMU
samples. The global flux binning containing all νe and νμbins is denoted as i ¼ 1; …; 23.
We obtain the measured flux values Φi observed at thedetector position from the predicted να þ ν̄α flux valuesαΦνMCðĒi Þ (where α ¼ e, μ), the expected number of CC
interactions N CC
MC;i , and the number of CC events obtainedby unfolding the data N CC
i . Considering the difference indetector configuration by each operational period, we
perform an unfolding separately for each period denotedby the subscript SK, so
PSK CCNναΦi ¼ ΦMC ðĒi Þ · PSK CCi;SK
ð3:7Þ
N MC;i;SK
PSK
where
is a sum over the four SK periods. The
αpredicted flux values ΦνMCðĒi Þ are themselves calculated
as follows. First, the predicted unoscillated να and ν̄αα ;ν̄αdifferential fluxes ϕνMCðEν ; θz ; ϕA Þ as a function of neutrino energy Eν, and zenith and azimuthal angles of arrival

PHYSICAL REVIEW D 94, 052001 (2016)

direction θz and ϕA , are calculated by interpolating thetabulated HKKM11 flux in energy and direction as
described in [20]. The predicted oscillated flux at theαdetector position ΦνMCðEν Þ is then calculated by integrating
the differential flux over a 4π solid angle asαðEν Þ ¼
ΦνMC

Z
4π

dΩ

e;μ
X
νβ
fϕMC
ðEν ; θz ; ϕA ÞOνβ →να ðEν ; θz ; θO Þ
β

ν̄β
þ ϕMC
ðEν ; θz ; ϕA ÞOν̄β →ν̄α ðEν ; θz ; θO Þg

ð3:8Þ

where Oνβ →να ðEν ; θz ; θO Þ is the νβ → να oscillation probability from the production point to the detector, calculatedwith the standard three-flavor oscillation model withparameters θO as shown in Table III and including thematter effect when propagating inside the Earth. Finally, thepredicted flux values at the mean energy of the ith energy
bin Ēi are calculated according to [73], using and the lower
and upper edges of the ith energy bin Ei and Eiþ1 , asαðĒi Þ ¼
ΦνMC

Z E

1
Eiþ1 − Ei

iþ1

Ei

α
ΦνMC
ðEν ÞdEν ;ð3:9Þ

which are the values appearing in Eq. (3.7).
For comparison, the predicted neutrino fluxes
without oscillation are also calculated by replacingOνβ →να ðEν ; θz ; θO Þ ¼ 0ð1Þ for α ≠ β (α ¼ β) in Eq. (3.8).
Figure 4 shows the comparison between the predicted νμflux with and without oscillation, also showing the upwardgoing and downwardgoing flux components separately. It
can be seen that the oscillation effect is different between
the upward and downward fluxes, since the neutrino path
length is not a linear function of the zenith angle andchanges rapidly near the horizontal direction.
TABLE III. Nominal value and 1σ uncertainty (except for CPviolation phase and mass hierarchy) of the oscillation parameters
based on [35]. The unit of mass square difference is eV2 . The
probability density distributions of the error are assumed to be
Gaussian distributions while the uniform probability is considered for δCP and mass hierarchy. The correlations among theseparameters are not taken into account.
Parameter

NominalUncertainty

Δm221
jΔm232 jsin2 θ12
sin2 θ23
sin2 θ13

7.50 × 10−5

ð7.30–7.70Þ × 10−5
ð2.20–2.44Þ × 10−50.296–0.333
0.388–0.612
0.0217–0.0263
0–2π
Both > 0 and < 0δCP
Sign of Δm232

052001-9−3

2.32 × 10
0.31
0.50
0.024
0
>0

E. RICHARD et al.

PHYSICAL REVIEW D 94, 052001 (2016)

10-1

νμ

-2

E2Φ [GeV/cm2/sec/sr]

10

10-3
10-4

Total

10-5

upward (x0.1)

10-6

w/ oscillation
w/o oscillation

downward (x0.01)

-7

10

10-8
-1

-0.5

0

0.5

1

1.5

2

2.5

3

3.5

4

Log (Eν/GeV)
10

FIG. 4. Calculated νμ flux in Eq. (3.8) at the SK site, based on
the HKKM11 flux model with and without applying neutrino
oscillation. The calculations of only upwardgoing (cos θz < 0)and downwardgoing (cos θz > 0) are also shown separately with
their normalizations scaled by ×0.1 and ×0.01.

C. Systematic uncertainties
The systematic uncertainties on the number of events ineach data bin Mj , which we define to include the uncertainties on the MC prediction, are estimated as follows.
Using a systematic error database maintained by the SK
Collaboration, the effects of errors relating to the event
selection and reconstruction in the detector (116 errors),
and those relating to neutrino interaction (17 errors) can becalculated. For the detector-related errors, the same 29
systematic uncertainty sources described in [8] are considered and estimated for each SK period, giving 116 total
errors. Their effects are evaluated as systematic error
coefficients f jk , representing the fractional shift in eachdata bin M j resulting from a 1σ shift of the kth systematic
error source. The modified expectation of the number of~ j , is thus
events in each bin, summing over each error, M


N sys
X
~
f jk gk ;Mj ðgÞ ¼ MMC;j × 1 þð3:10Þ

kwhere MMC;j are the nominal MC expectations without any
systematic effect, N sys is the number of systematic error
sources, and g ¼ ðg1 ; …; gN sys Þ represents the appliedstrength of each systematic in units of σ. The effects of
the oscillation parameters are treated separately, being
directly calculated by shifting their values in Eq. (3.8).
The uncertainties of the oscillation parameters are shown in
Table III.
For propagation of both the systematic and statisticaluncertainties from Mj to the flux values N CC
i we employ a
toy MC method, since while ROOUNFOLD provides anaccurate theoretical calculation of the statistical errors[72], simultaneous treatment of the systematic uncertainties
within the software is not implemented. In this method, toydata sets are generated by randomly fluctuating MC dataaccording to their systematic and statistical error probability density functions (PDFs), as follows: first, system~ j ðgÞ are generated using a random set
atically shifted data M
of N sys systematic strengths g, where the PDF for each gk
parameter is a Gaussian distribution and there is nocorrelation between the gk parameters. The final toy data
Mtoy
j is then generated by a Poisson distribution with the~ j ðgÞ. Two thousand sets of toy data aremean values of M
generated and analyzed with the same analysis method as
described in the previous section. The variations among
these toy flux measurements are taken as the uncertainties,
including both statistical and systematic uncertainties; thecovariances of the unfolded flux among the energy bins arealso recorded into the covariance matrix, which will be
used for the χ 2 calculation in Sec. III D. The method is also
repeated for statistical errors only, to allow for the separate
calculation of the statistical and systematic components.
A final contribution to the uncertainty in the values N CC
icomes from the regularization error, i.e., the bias that may
be caused by inaccuracy in the initial estimation of the
response matrix and Bayesian prior. As mentioned inSec. III B, such a bias should be small, but it may be
noticeable in unfolded bins with low statistics. This bias
cannot be corrected for exactly, but we may estimate an
associated error by unfolding pseudodata sets with anenergy spectrum reasonably far from the MC prediction,and observing the resulting difference. Such pseudodata
sets were produced from the MC data by reweighting
events to define a modified energy spectrum Φ0MC;i asΦ0MC;i ¼ ð1 þ ΔαÞ

Ēi
1 GeV

Δγ

ΦMC;i

ð3:11Þ

where Δα and Δγ represent the deviation of the normalization and spectral index from each model, and are modified
in the ranges of 0.05 and 20%, respectively, and ΦMC;i isshorthand for the predicted flux values ΦMC ðĒi Þ. In Fig. 5
the comparisons of each input pseudodata set and its
unfolded flux spectra are shown, along with the fractional
deviation between the input and unfolded output in eachcase. The difference becomes largest at the highest energy
bin for both νe and νμ spectra, and the error sizes are 6%and 8% at most in the case of a 0.05 and 20% change.
These differences are taken as the estimated uncertainties
due to regularization in the unfolded energy spectrum, andare included in the systematic uncertainties.
Figure 6 shows the breakdown of the estimateduncertainties for each error source group. The absolute
errors of four groups (statistical, neutrino interaction,
detector response, and neutrino oscillation þ regularization), are individually calculated and compared. Currently
the neutrino interaction errors are the dominant source ofsystematic uncertainty.

052001-10

MEASUREMENTS OF THE ATMOSPHERIC NEUTRINO FLUX …
1.8

Φ / ΦMC

1.2

1
0.8

1
0.8
0.6

Δ α = +0.1, Δ γ = -0.025
Δ α = -0.1, Δ γ = +0.025
Δ α = +0.2, Δ γ = -0.05
Δ α = -0.2, Δ γ = +0.05

0.2

0.1 -1
0.05
0
-0.05
-0.1

-0.5

0

0.5

-1

-0.5

0

0.5

1

1

Δ α = -0.1, Δ γ = +0.025
Δ α = +0.1, Δ γ = -0.025
Δ α = -0.2, Δ γ = +0.05
Δ α = +0.2, Δ γ = -0.05

0.4
0.2

1.5

2

2.5

3

3.5

4

1.5

2

2.5

3

3.5

4

DeviationΦ / ΦMC

1.4

1.2

0.6

νμ

1.6

1.4

0.4

Deviation

1.8

νe

1.6

PHYSICAL REVIEW D 94, 052001 (2016)

0.1 -1
0.05
0
-0.05
-0.1

-0.5

0

0.5

-1

-0.5

0

0.5

Log (Eν /GeV)

1

1

1.5

2

2.5

3

3.5

4

1.5

2

2.5

3

3.5

4

Log (Eν /GeV)

10

10

FIG. 5. Pseudodata input and the unfolded energy spectra (upper) and their ratio (lower) for νe (left) and νμ (right), respectively.
Vertical axis in the upper panels is the ratio to HKKM11 [20] (ΦMC ) in the upper figure. Four sets of pseudodata inputs, where thedeviations of flux normalization and spectral index ðΔα; ΔγÞ ¼ ð10%; 0.025Þ, ð20%; 0.05Þ, are tested.

D. Results and discussionsΔΦ/ Φ [%]

25
20

E2Φ [GeV/cm2/sec/sr]

νe

10-5

1.6
1.4
1.2
1
0.8
0.6

νμ

-1

-0.5

νμ

νe

10-4

νe

15
10

νμ

10-3

1.6
1.4
1.2
1
0.8
0.6

total
statisticalneutrino interaction
detector
oscillation and regularization

30

10-2

Ratio

35

10-1

Ratio

Figure 7 shows the obtained νe and νμ energy spectrausing all SK-I–SK-IV data. The energy binning, meanenergy, measured flux, and error are also described inTable IV for νe and νμ . The measured energy spectrum
agrees with the oscillated HKKM11 flux within the estimated uncertainties, which, as mentioned above, are dominated by neutrino interaction uncertainties. The unoscillated
flux is also plotted, such that the deficit of νμ flux due toneutrino oscillation becomes apparent below 100 GeV.
The observed fluxes are compared to several flux
models, including HKKM11 [20], HKKM07 [19],FLUKA [22], and Bartol [21], in Fig. 8. In order to

0

0.5

1

1.5

2

2.5

3

3.5

4

Log (Eν / GeV)
10

5
0

0

2

4

6

8

10
12
i bins

14

16

18

20

22

FIG. 6. Breakdown of the uncertainties in measured flux
estimated for each error group, showing neutrino cross section,
detector response, and neutrino oscillation þ regularization.

FIG. 7. Measured flux spectra, using all SK-I–SK-IV data, forνe (blue points) and νμ (red points) as a function of the logarithmof the neutrino energy. Error bars correspond to the uncertainties,
including all statistical and systematic errors. The νe (blue line)and νμ (red line) flux curves show the HKKM11 [20] model with(solid) and without (dashed) neutrino oscillation. In the lower partof the figure, the data-to-HKKM11 ratio is shown.

052001-11

E. RICHARD et al.

PHYSICAL REVIEW D 94, 052001 (2016)

TABLE IV. Electron and muon neutrino fluxes measured by
SK-I-SK-IV data. Error written in percentage, including both
statistical and systematic uncertainties.
i

log10
ðE=GeVÞ

log10
ðĒi =GeVÞνe
1 −0.8 to −0.6
2 −0.6 to −0.4
3 −0.4 to −0.2
4 −0.2 to 0.0
5 0.0 to 0.2
6 0.2 to 0.4
7 0.4 to 0.6
8 0.6 to 0.8
9 0.8 to 1.0
10 1.0 to 1.5
11 1.5 to 2.0
νμ12 −0.6 to −0.4
13 −0.4 to −0.2
14 −0.2 to 0.0
15 0.0 to 0.2
16 0.2 to 0.4
17 0.4 to 0.6
18 0.6 to 0.8
19 0.8 to 1.0
20 1.0 to 1.5
21 1.5 to 2.0
22 2.0 to 3.0
23 3.0 to 4.0

ΔΦνi =Φνi(%)

Ē2i Φνi(GeV=cm2 = sec =sr)

−0.71
−0.51
−0.27
−0.09
0.10
0.30
0.50
0.70
0.90
1.22
1.72

1.21 × 10−2
1.46 × 10−2
1.50 × 10−2
1.37 × 10−2
1.16 × 10−2
8.55 × 10−3
6.09 × 10−3
3.73 × 10−3
2.32 × 10−3
9.42 × 10−4
2.03 × 10−4

18
17
16
15
17
17
18
19
18
15
18

−0.51
−0.32
−0.09
0.10
0.30
0.51
0.71
0.90
1.21
1.73
2.40
3.39

1.58 × 10−2
1.77 × 10−2
1.86 × 10−2
1.68 × 10−2
1.38 × 10−2
9.59 × 10−3
6.68 × 10−3
4.79 × 10−3
2.62 × 10−3
1.20 × 10−3
2.49 × 10−4
1.46 × 10−5

21
16
15
16
18
19
19
17
13
16
18
21

Ratio
Ratio

νe

1.6
1.4
1.2
1
0.8
0.6

νμ

-1

HKKM11

BartolHKKM07

FLUKA

1
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
0

2

4

6

8 10 12 14 16 18 20 22-1

FIG. 9. Error correlation matrix, Ĉij , of the observed spectrum
bins, used for the χ 2 calculation. The matrix contains both νe andνμ , with νe ranging from the 1st to the 11th bin, and νμ from the12th to the 23rd bin. Strong correlations are seen between νe andνμ as many systematics apply to both flavors in similar ways.

χ2 ¼

N X
N
X
i

perform a quantitative comparison, χ 2 values are calculated
between the observed energy spectrum and the flux model
predictions while taking into account the error correlations
between the energy bins,
1.6
1.4
1.2
1
0.8
0.6

22
20
18
16
14
12
10
8
6
4
2
0

ðΦi − ΦMC;i ÞT C−1
ij ðΦj − ΦMC;j Þ

ð3:12Þ

jwhere N is the number of the data bins, Φi and ΦMC;i are theobserved flux and the expectation of the chosen flux modelat the ith energy bin, respectively, and Cij is the covariancematrix of the observed spectrum Φi (which is calculated byobserving the correlation of the error fluctuations for N CCi
within the set of toy MC data described in Sec. III C).
Figure 9 shows the correlation matrix, where each elementĈij is obtained by normalizing the diagonal term of thecovariance matrix Cij to 1, i.e., Ĉij ¼ Cij =ðΔΦi ΔΦj Þ,
where ΔΦi is the error of the observed flux at the ithbin shown in Table IV. The value of each element of thecorrelation matrix is given in the Appendix. Both statisticaland systematic uncertainties are taken into account in theχ 2 calculation.
The calculation is performed for three cases: both νe andνμ (degrees of freedom N ¼ 23), νe only (N ¼ 11), and νμonly (N ¼ 12); the obtained χ 2 s for each flux model are
shown in Table V. The test including νe and νμ together
takes the systematic correlations between the two flavors
TABLE V. χ 2 values calculated by testing the measured flux
against each flux model prediction according to Eq. (3.12). The
number of degrees of freedom (DOF) in each test is also shown.

-0.5

0

0.5

1
1.5
2
Log (Eν / GeV)

2.5

3

3.5

4

χ2

10

FIG. 8. Measured energy spectra compared to the flux model
predictions (dotted lines) from the HKKM11 [20] (black),
HKKM07 [19] (light green), Bartol [21] (red), and FLUKA[22] (light blue) models. The vertical axis represents the flux ratio
from the data (or each model) to HKKM11. Error bars include
both statistical and systematic uncertainties.

Flux modelνe and νμνe onlyνμ only

HKKM11 [20]HKKM07 [19]Bartol [21]FLUKA [22]
DOF

22.2
22.5
30.7
25.6
23

5.3
6.8
6.6
5.4
11

12.2
12.1
17.0
15.2
12

052001-12

0.2

Ratio

Normalization ( Δ α)

MEASUREMENTS OF THE ATMOSPHERIC NEUTRINO FLUX …

0.1
0
-0.1

νμ

νe

-0.2

Ratio

Spectral index (Δ γ )

0.06
0.04
0.02
0

-0.02

PHYSICAL REVIEW D 94, 052001 (2016)
1.6
1.4
1.2
1
0.8
0.6

νe

1.6
1.4
1.2
1
0.8
0.6

νμ

upwarddownward

upwarddownward

-0.04

-1

-0.06

Ratio
Ratio

νμ

-1

-0.5

1

SK1
SK2

1.5

2

2.5

3

3.5

4

FIG. 12. Ratio of the measured flux to the HKKM11 model,
separated into upwardgoing (blue) and downwardgoing (red)data. Error bars show the statistical uncertainty only.

consistency of this structure across each SK period, theunfolded energy spectrum of each period is compared inFig. 11. Though there are larger statistical fluctuations(especially in the high-energy regions in νe due to thesmaller number of observed events), the spectral shapeseems consistent across all SK periods within the statisticalerrors. This shape is, however, consistent within the rangeof the estimated systematic effects, in particular theneutrino interaction uncertainties.

1.1
1
0.9
0.8
0.7
0.6

Residual from input

0.7
0.6
0.5
0.4
0.3
0.2
0.1

0.5

0

0.4

-0.1

1.2
1.1

νμ

0.7
0.6

1

0.5

0.9

0.4

0.8

0.3

0.7
0.6

SK3
SK4

Input
1st
2nd
3rd
4th
5th
6th
8th
10th

νe

0.2
Residual from input

0.1

0.5

0

0.4
-1

Residual

1.2

Ratio

1.6
1.4
1.2
1
0.8
0.6

0.5

10

into account. In this test the χ 2 values are not strongly
inconsistent for any of the flux models, though HKKM11
has the smallest values among these models; the p values
are p ¼ 0.51 for HKKM11, p ¼ 0.32 for FLUKA, and
p ¼ 0.13 for Bartol. As HKKM11 is preferred above
HKKM07, we see that the updates to the hadron simulation[20] in the HKKM11 model (which cause changes in theenergy region below 1 GeV) seem to bring better agreement. For the individual tests, both νe and νμ agree well.
We also evaluate each model against the data by fitting
an energy spectrum with variable parameters, Δα and Δγ,
which represent the deviations in normalization and spectral index from each model, respectively, in the same way as
was defined in Eq. (3.11). The χ 2 statistic is defined by
Eq. (3.12) with ΦMC;i replaced with Φ0MC;i , and the best-fit
values and 1σ errors of Δα and Δγ parameters are obtainedby minimizing the χ 2 . Figure 10 shows the result for each
flux model case. The normalization and spectral index
agrees within the 1σ error for every model, except from thefitted spectral index of FLUKA νμ which deviates by 2.7σ.
We notice that a similar wavy structure is seen for boththe νe and νμ unfolded fluxes. In order to verify theνe0Log (Eν / GeV)

FIG. 10. Best-fit and 1σ error of Δα and Δγ parameters for each
flux model by χ 2 calculation, representing the deviation in normalization and spectral index from the flux prediction, respectively.

1.6
1.4
1.2
1
0.8
0.6

-0.5

FLUKA

Residual

BartolRatio

HKKM11 HKKM07

-0.1
-0.5

0

0.5

1

1.5

2

2.5

3

3.5

4

Log (Eν /GeV)
10

0

0.5

1

1.5

2

2.5

3

3.5

4

Log (Eν / GeV)
10

FIG. 11. Ratio of the measured flux to the HKKM11 model foreach SK period: SK-I (green), SK-II (blue), SK-III (cyan), and
SK-IV (magenta). Error bars show the statistical uncertainty only.

FIG. 13. Validation of the unfolding procedure, showing theunfolded spectra at each iteration (dashed lines with different
colors) compared to the truth input (solid line with dots), which inthis case is set identically to our spectrum unfolded from the data.
The number of iterations was set to 5 in this paper (red dashedline with asterisks). The vertical axis on the left side shows theratio of the unfolded flux to the HKKM11 flux, and on the rightside shows the residuals compared to the input spectrum. The
shaded area shows the size of the previously estimated regularization error.

052001-13

E. RICHARD et al.

PHYSICAL REVIEW D 94, 052001 (2016)The νe and νμ fluxes are also separated into upwardgoingand downwardgoing data sets, using their reconstructed
direction, and measured in Fig. 12. This is in order to checkany possible bias in the flux calculation due to thedifferences of the data sets; for example, neutrino oscillation has a stronger effect in the upwardgoing data above
GeV energies, and the UPMU data is an upwardgoingsample only. As seen in the figure, no obvious difference
exists between the fluxes measured using these twodata sets.
Although we thoroughly validated the accuracy of theunfolding procedure using the HKKM-like pseudodata andincluded the estimated regularization bias as a systematic
uncertainty, there has still been some concern about theability of the unfolding procedure to accurately reproduce
more complicated spectral shapes, such as the wavy shapethat was eventually obtained, and, in particular, whether or
not such a shape would be more strongly affected by the

number of iterations. We therefore perform a further post
hoc check using the same validation method as before, butusing our actual unfolded spectra from the data as thepseudodata truth input. Figure 13 shows the unfolded
spectrum as a function of number of iterations, from1 up to 10, on top of the previously estimated regularization
error. It is seen that at around five iterations, which we had
adopted based on the pseudodata test, the unfolded spectra
are stable. Furthermore, the spectra are reproduced approximately within the estimated regularization uncertainties.
Therefore, we conclude that the shapes of our unfolded
spectra are not due to an unexpected additional bias fromthe unfolding procedure.
In Fig. 14 this flux measurement is shown with theresults from other experiments. Our measured data provide
significantly improved precision below 100 GeV. The
minimum of the observed energy range is extended
below 1 GeV, and at higher energies overlaps with νμ

10-1
10-2

E2 Φ [GeV cm-2 sec-1 sr-1]

10-3
10-4
-5

10

Super-Kamiokande I-IV ν μFrejus ν μ
IceCube ν μ unfolding
IceCube ν μ forward folding
AMANDA-II ν μ unfolding
AMANDA-II ν μ forward folding
ANTARES ν μ
HKKM11 ν μ+ν μ (w/ osc.)

10-6
-7

10

Super-Kamiokande I-IV ν eFrejus ν e
IceCube/DeepCore 2013 ν e
IceCube 2014 ν e
HKKM11 ν e+ν e (w/ osc.)

10-8
10-9

-1

0

1

2
3
Log (Eν /GeV)

4

5

10

FIG. 14. The measured energy spectra of the atmospheric νe and νμ fluxes by SK, shown with measurements by other experiments,
i.e., Frejus [39], AMANDA-II [40,41], IceCube [42–45], and ANTARES [75]. The phrase “forward folding” used by IceCube and
AMANDA-II is synonymous with forward fitting. The HKKM11 flux model predictions for the Kamioka site are also shown in solid(with oscillation) and dashed (without oscillation) lines. The error bars on the SK measurement include all statistical and systematic
uncertainties.

052001-14

MEASUREMENTS OF THE ATMOSPHERIC NEUTRINO FLUX …

PHYSICAL REVIEW D 94, 052001 (2016)

measurements by AMANDA-II, IceCube, and ANTARES,
which should allow a better constraint on the flux normalization at the energies beyond 100 GeV. The measured
fluxes are consistent with these measurements in theoverlap energy region greater than 100 GeV. There is somediscrepancy with the Frejus measurement below 1 GeV;
however, due to the difference in geomagnetic cutoff at the
Frejus site, the low-energy flux is expected to be a few tensof percent higher at that location [74].

B. Analysis method

A. Data sample
The azimuthal analysis uses only the FC sample
described in Sec. II B, and selects sub-GeV and multiGeV e-like and μ-like events with a single reconstructed
Cherenkov ring. For the sub-GeV e-like events, the NC π 0 like events identified by the π 0 finder algorithm areremoved. These samples have comparatively high statistics
among the SK data samples, as seen in Table I in
Sec. II C, and cover a wide energy range from 100 MeV
to ∼10 GeV, which includes the range where the east-west
asymmetry is expected to be visible. The single-ring events
have the highest accuracy for vertex and directional
reconstruction, and the neutrino flavor accuracy in thissample is estimated as 93.6% for νe and 98.5% for νμ, CC
and NC inclusive.
The azimuthal angle ϕ is defined by a particle’s forwardgoing direction, clockwise from true south at the SK
detector position (i.e., ϕ ¼ 0° for southgoing, ϕ ¼ 90°
for westgoing, ϕ ¼ 180° for northgoing, and ϕ ¼ 270°for eastgoing). The zenith angle θ is given in terms of cos θ,
where −1 represents an upwardgoing particle. Each neutrino event is binned using the reconstructed properties ofthe produced lepton, including azimuthal angle ϕrec (binned
in 12 evenly sized bins), zenith angle θrec (binned in 5
evenly sized bins in cos θrec ), and energy Erec. The
reconstructed energy Erec was described in Sec. III A,
and is binned as 4 energy bins starting at 0.1 (0.2 for μlike), 0.4, 1.33, and 3.0 GeV, where the last bin isunbounded from above. This binning was chosen following
the MC prediction, in order to show the zenith and energy
regions in which the predicted east-west effects havevarious different strengths, while keeping statistics highin each bin.
All data from the SK periods I–IV are summed together,
for a total of 13 079 e-like and 12 725 μ-like events. The
MC events, as described in Sec. II C, are binned in the sameway using their reconstructed variables, with weights forthe data live times, solar activity, and neutrino oscillation.
In this section we consider, unless otherwise stated, the MC
events generated according to the HKKM11 flux model,
which is the only model that publishes comprehensive
azimuthal distributions. The HKKM11 model uses the
IGRF-10 geomagnetic field model [76].

A¼

neast − nwest
neast þ nwest

ð4:1Þ

where nwest represents the total number of eventswith azimuth angle 0–180 degrees, i.e., westgoing events,
and neast represents eastgoing events with azimuth angle
180–360 degrees.
To calculate the significance of a nonzero east-west
effect, we test a reduced sample of events in the energy
range from 0.4 to 3.0 GeV, with j cos θrec j < 0.6. These
criteria are optimized by the MC method to select theevents giving the largest predicted value of A=ΔA, which isdefined as the significance in units of σ.
By examining the rigidity cutoff at SK as a function ofboth ϕ and θ, shown in Fig. 15, we can see that while thesimple dipole effect is clearly visible near the horizontal, in
1
10
10
20

cos θ

IV. AZIMUTHAL SPECTRUM ANALYSIS

Unlike the energy spectrum analysis, a forward-fitting
method is used instead of an unfolding procedure. This isbecause the azimuthal distributions are unique to the SKsite and will not be directly compared with the results ofother experiments, and as such we can compare our data
distributions with the MC predictions, without needing toreconstruct the true flux distributions. The azimuthal
distributions are plotted first using events in each energy
bin while summing over all zenith angle bins, then for eachzenith bin while summing over all energy bins. To quantify
the east-west dipole asymmetry in each plot, we define theparameter

0

50
10

20

30
5
5
20

10

5

5

10

-1
180

90

0
Azimuth φ [deg]

270

180

FIG. 15. Contour map of the rigidity cutoff, seen from theperspective of the SK detector. The red lines show the projectionof the Earth’s landmasses, and the blue contours and labels
represent the rigidity cutoff in GV. While the axes are flippedcompared with other plots in this paper, the labeling of theazimuth ϕ and zenith θ angles uses the same convention.
Reproduced with permission from T. Gaisser and M. Honda,
Annu. Rev. Nucl. Part. Sci. 52, 153 (2002). Copyright 2002,
Annual Reviews.

052001-15

E. RICHARD et al.

PHYSICAL REVIEW D 94, 052001 (2016)

general the structure of the rigidity contours is more
complicated. This is due to the difference between geomagnetic and geophysical north, and irregular contributions
from the nondipole components of the geomagnetic field[77]. In order to test if this predicted structure is visible at
SK within the limited statistics available, we define theparameter B by fitting the function
k1 sinðϕ þ BÞ þ k2

ð4:3Þwhich by Wilks’s theorem [78] should be distributed as a
χ 2k¼2 (where k denotes the number of degrees of freedom),since the two hypotheses have a difference in parameter
space dimension of 2. The significance of rejecting a
constant B parameter in favor of the MC expectation isthen defined in units of σ by comparing the obtained p
value
Z ∞
p¼
χ 2k¼2 ðxÞdx
ð4:4Þ
Δχ 2

with the normal distribution. By toy MC, we expect a2.0σ  0.9σ indication of a zenith-dependent B parameter,
as shown in Fig. 16.
C. Systematic uncertainties
Because the shape of SK detector is azimuthally
symmetric, any azimuthal asymmetry in the detector’sreconstruction abilities is expected to be small.
Nonetheless, the possible existence of reconstruction errors
directly dependent on the azimuthal angle is investigated.
The dependence of reconstructed momentum on the

Fraction of toy MCs

0.03
0.025
0.02
0.015
0.01

ð4:2Þ

to the azimuthal distribution of the data in each zenith bin,where k1 , k2 are free parameters. This parameter essentiallymeasures any phase shift in the azimuthal alignment of theasymmetry dipole, depending on zenith angle.
To test the significance of the existence of such zenithdependent effects, we perform a χ 2 test of the measured B
parameters across the zenith bins, first against a constant(zenith-independent) B parameter, then against the nominal
MC shape. For the first test, the normalization of theconstant B parameter is fitted to the data for both νe and νμevents separately; then, the χ 2 statistic is taken and denotedas χ 2flat . For the second test, the normalization of the MC
distributions for the νe and νμ events is again fitted, andfurthermore the amplitude of the effect (relative differenceof B in each zenith bin from the average in all bins) is freely
scaled for νe and νμ events, which allows for a stronger or
weaker zenith dependency of the dipole alignment than
predicted by MC simulations. The result of the MC shape
fit is denoted χ 2MC . We then defineΔχ 2 ≡ χ 2flat − χ 2MC

0.04
0.035

0.005
0
0

1

2

3

4
5
6
Significance [σ]

7

8

9

10

FIG. 16. The significance σ of rejecting a non-zenith-dependent
B parameter in favor of the MC expectation, estimated by the
result of 10 000 toy MC simulations.

azimuth angle is checked using the momentum of theMichel electrons coming from the decay of cosmic muons
which stopped in the detector, and is found to be uniform tobetter than 0.6%. While the energy range of Michelelectrons is lower than that of our event samples, thisvalue should give a reasonable indication of the maximum
energy scale asymmetry and is implemented as a systematic
error. Other possible azimuthal biases, such as particle
identification or ring counting, are tested by comparing theshape of the likelihood variable distributions between eastand westgoing single-ring events. The maximum shift inevents between classifications is found to be less than 0.1%
in each case, and is considered negligible [60]. Nonneutrino background events in the SK data, such ascosmic-ray muons and other noise events, are estimatedto represent 0.2% of events in the SK-I to SK-IV data set[60], and any azimuthal asymmetry present in this background is considered negligible.
Considering the neutrino flux at SK, the geomagnetic
effects are expected to cause changes in the energy
spectrum, neutrino path length, and neutrino/antineutrino
flux ratio, depending on the azimuth angle. Any systematic
error that depends on these factors should thus be considered. To give one clear example of this, we show in Fig. 17
the MC true energy spectra of μ-like sub-GeV events,
separated by east- and westgoing directions, where we
see that the eastgoing flux is relatively skewed towards
lower energies. The “NEUT axial mass parameter”
systematic is at the ∼10% level when averaging over all
events, but increases sharply for low-energy events in the
Oð100Þ MeV range; the eastgoing flux bins are thus
assigned a systematic error due to uncertainty in theneutrino cross section up to 3 percentage points higher
than the westgoing flux bins, due to their high proportionsof low-energy events.
The effects of all applicable errors from the SK
systematic error database (described in Sec. III C) arecalculated in the azimuthal binning, including cross section

052001-16

MEASUREMENTS OF THE ATMOSPHERIC NEUTRINO FLUX …

TABLE VI. Each systematic error source and the strength oftheir azimuthally dependent component, defined as the average
percentage shift from the nominal MC values in each bin whenbinned as in Fig. 20, for a 1σ shift in the error source. The
Monte Carlo statistical error does not cancel between east- andwestgoing events, but is small, less than 0.5%. “Others” represents the sum of 55 systematics with effects less than 0.05%,
combined in quadrature, and “Total” represents all errors combined in quadrature.

east-going

0.05

west-going
0.04

0.03

0.02

0.01

0
200

PHYSICAL REVIEW D 94, 052001 (2016)

Systematic error source
400

600

800

1000

1200

Eν [MeV]

FIG. 17. Energy spectra between eastgoing and westgoing
events, showing the MC true energy spectrum of events reconstructed as single-ring sub-GeV μ-like events. The area under bothcurves is normalized to unity.

errors (15 errors), reconstruction-related detector errors(53 errors), neutrino oscillation model errors (6 errors),
and flux model errors (21 errors). The detector-related
errors are numerous as most of them are modeled separately
for each SK period; however, because the selection of
single-ring events is simpler than the multiple samples usedin the energy spectrum study, fewer reconstruction-related
errors are applicable than in Sec. III C.
For the parameter A in Eq. (4.1), any systematic shift thatapplies equally to nwest and neast will cancel out; only
systematics that effect the two unequally, i.e., an azimuthally dependent component, will contribute to the final
systematic error. Similarly, when we plot the full azimuthal
distributions, we are interested in testing their shape as afunction of azimuth, rather than measuring the absolutenormalization of the neutrino flux. Thus, for each systematic error, we define its normalization component as theaverage effect on all azimuthal bins, and its azimuthally
dependent component as the effect on each bin relative tothe average effect. The normalization component is discarded, and the azimuthally dependent components of all
systematics are combined in quadrature to calculate thefinal systematic error in each azimuthal bin.
The systematic errors with the strongest effects on theazimuthal analyses are listed in Table VI. Cross-section
errors are still dominant, and there are some contributions
from relative normalizations, even though these errors are
only weakly azimuthally dependent. Some detector errorsof SK-IV and SK-I are also noticeable, due to their high live
times compared with other SK periods. However, the total
systematic effects are small, at the 1% level.
Finally, it is important to note that the dominant factor inmisreconstructing the neutrino direction does not comefrom detector inadequacies, but from the fact that incharged-current quasi-elastic (CCQE) interactions at lower

NEUT axial mass parameter
MC statistical error
Cross section ratio in CC quasielastic
Flux relative normalization (<1 GeV)
Cross section in single meson production
Flux relative normalization (>1 GeV)
Cross section ratio in neutral/charged-pion
production
Coherent pion production cross section
Deep inelastic scattering Q2 at low W
Relative normalization for sub-/multi-GeV FC
ν=ν̄ ratio in single-pion productionν=ν̄ ratio in CC quasielastic
Deep inelastic scattering model differences
Δm232 error (taken from T2K)
Azimuthal energy calibration (SK-IV)
Fiducial volume (SK-IV)
ν=ν̄ ratio (1 < Eν < 10 GeV)
Overall energy calibration (SK-IV)
Azimuthal energy calibration (SK-I)
Others
Total

Size of effect
(%)
0.59
0.49
0.46
0.34
0.34
0.25
0.22
0.17
0.10
0.10
0.08
0.08
0.08
0.07
0.06
0.06
0.05
0.05
0.05
0.15
1.01

energies the neutrino direction and the produced lepton
direction are poorly correlated. This is not modeled as asystematic error, but the effects are accounted for by thehigh statistics of the SK MC simulations model in theforward-fitting procedure. The correlation for reconstructed
single-ring events is plotted by energy in Fig. 18, whichshows that for interactions at less than 400 MeV there is
only a faint ability to even discriminate the neutrino’s
forward direction from its backwards one. Events in therange 400 < Eν < 1330 are generally correlated at least
within 90°, which allows a reasonable separation between
east- and westgoing events. Finally, multi-GeV events arevery well correlated, within ∼15°.
D. Results and discussion
To understand the origin of the reconstructed distributions from the true flux shape, we first show the azimuthal
distributions of MC events binned by true neutrino direction in Fig. 19, where we denote the true neutrino azimuthaland zenith angles as ϕν and θν respectively. The azimuthal
distributions binned by the lepton’s reconstructed direction

052001-17

E. RICHARD et al.

PHYSICAL REVIEW D 94, 052001 (2016)0.1

Eν < 400
400 < Eν < 700
700 < Eν < 1000
1000 < Eν < 1330
1330 < Eν

0.08

0.06

0.04

0.02

0

0

20

40

60

80

100

120

140

160

180

ν → l inner angle [deg]FIG. 18. Directional correlation between the true neutrino andlepton directions for the single-ring e-like events used in theazimuthal analysis, separated by true neutrino energy Eν given in
MeV. The correlation for μ-like events is not shown, but isgenerally similar. The area under each curve is normalized tounity.

ϕrec and θrec are shown in Fig. 20, for both data and MC
events. Both of the above figures bin events using thereconstructed particle energy Erec. The values of the dipole
asymmetry parameter A, corresponding to each plot ofFig. 20, are shown in Fig. 21.

600

0.1 < Erec < 0.4 [GeV]Considering first the MC events, we see that at
Erec < 0.4 GeV, the dipole asymmetry is high in the true
neutrino direction, but not in the reconstructed lepton
direction; this is due to the poor correlation between
the incoming neutrino and outgoing lepton directions atthese energies. The asymmetry in the reconstructed
lepton distributions is then expected to peak at around
Erec ∼ 1 GeV, and also for events at around cos θrec ∼ 0. The
prediction of the asymmetry strength is somewhat higher
than in the previous SK analysis [50], at which time the flux
simulation considered geomagnetic effects on the primarycosmic-ray particles only; the additional consideration of thebending of secondary cosmic rays in the atmosphere is seen
to further enhance the dipole asymmetry prediction [51]. In
the cos θ < −0.6 bins, the distribution of the true neutrino
direction also shows a different shape compared to thereconstructed lepton directional distribution. The shape ofthe true distribution is due to the rigidity cutoff having a
suppressive effect in both east and west directions for upgoing(cos θ < 0) neutrinos. However, as this effect is strongest for
sub-GeV events within ∼10° of upgoing neutrinos, it is also
washed out after reconstruction using the lepton direction.
Considering next the distributions of the data, there is in
general excellent agreement with the MC expectations in
the azimuthal plots (Fig. 20). A χ 2 value is calculated based
1.0 > cos θ ν > 0.6

300

µ
e

200

400

100

200
0
800

0
400

0.4 < Erec < 1.33 [GeV]

200

600

100

400

0

200

150

400
300
200
100
0
400

100

300

0
200

1.33 < Erec < 3.0 [GeV]

50

200

0
100

100

80

3.0 < Erec [GeV]

0.2 > cos θ ν > -0.2

-0.2 > cos θν > -0.6

0
300

60

-0.6 > cos θν > -1.0

200

40

100

20
0
0

0.6 > cos θ ν > 0.2

300

60

120
180
240
300
True Azimuth φν [deg]

360

0

0

60

120

180

240

300

360

True Azimuth φ ν [deg]FIG. 19. Azimuthal distributions of the forwardgoing true neutrino direction, for e-like (blue) and μ-like (red) MC events. The eventswere selected as described in Sec. IVA. The left plots show reconstructed energy ranges summing over all zenith angles, and the rightplots show zenith angle ranges summing over all energies. The MC truth information is used instead of the reconstructed lepton direction
for both zenith θν and azimuthal angle ϕν binning, although the energy binning still uses the reconstructed energy Erec.

052001-18

MEASUREMENTS OF THE ATMOSPHERIC NEUTRINO FLUX …
520

0.1 < Erec < 0.4 [GeV]

480

240

PHYSICAL REVIEW D 94, 052001 (2016)
1.0 > cos θ rec > 0.6

220
200
180

440

160
140

400

260

550
500

0.4 < Erec < 1.33 [GeV]

220

450

180

400

140

350

300

300

160

220
180
280

120

80

200

110

160

3.0 < Erec [GeV]

100
90

260
240

80

220

70

200

60

180
160

60

120
180
24 0
300
360 0
Reconstructed Azimuth φ rec [deg]0.2 < Erec < 0.4 [GeV]

320

340

240

300

200

260

300

0.4 < Erec < 1.33 [GeV]

60

120
180
240
300
360
Reconstructed Azimuth φ rec [deg]1.0 > cos θ rec > 0.6

0.6 > cos θ rec > 0.2

260
220
180

450
400
350
200

-0.6 > cos θ rec > -1.0

280

380

650
600
550
500

-0.2 > cos θ rec > -0.6

240

100

420

0.2 > cos θ rec > -0.2

260

1.33 < Erec < 3.0 [GeV]

140

00.6 > cos θ rec > 0.2

2600.2 > cos θ rec > -0.2

2201.33 < Erec < 3.0 [GeV]

180

180
140

160
140

220

-0.2 > cos θ rec > -0.6

120
180

100
80
60

3.0 < Erec [GeV]

140
220

50

200

40

180

30

160

20

140
120

0

-0.6 > cos θ rec > -1.0

60

120
180
24 0
300
Reconstructed Azimuth φrec [deg]

360 0

60

120
180
240
300
Reconstructed Azimuth φrec [deg]360

FIG. 20. Azimuthal distributions of the single-ring e-like (top) and μ-like (bottom) events for the SK-I–SK-IV data (points) and MC
simulations (boxes). The error bars on the data represent the statistical error, while those on the MC simulations represent the systematic
error. The binning is the same as in Fig. 19; however, events are binned by the reconstructed properties of the lepton. A dotted line isdrawn in each plot at the average value of the MC values as a visual guide, and the plots are zero suppressed to show the data/MC shapecomparison in detail.

052001-19

E. RICHARD et al.

PHYSICAL REVIEW D 94, 052001 (2016)eparameter A

parameter A

0.14

0.12

0.1

0.1

0.08

0.08

0.06

0.06

0.04

0.04

0.02

0.02

0

0

-0.02

-0.02

-0.04

-0.04

µ

0.14
0.12
0.1
0.08
0.06
0.04
0.02
0
-0.02
-0.04
0.5

1

1.5

e

0.14

0.12

µ

0.14
0.12
0.1
0.08
0.06
0.04
0.02
0
-0.02
-0.04

2 2.5
3
Erec [GeV]

1

0.6

0.2

-0.2
cos

-0.6

-1

rec

FIG. 21. The parameter A depending on the reconstructed energy (left) and zenith angle (right), for e-like (top) and μ-like (bottom)
events, from the SK-I–SK-IV data (points with statistical error) and MC simulations (boxes with systematic error). A grey dotted line isdrawn at A ¼ 0, the case of no east-west asymmetry, as a visual guide. The highest energy bin in the left plots is unbounded from above.

on the weighted-unweighted χ 2 test [79] for the energy(zenith) distributions; this value is 87.6 for 96 bins (106.6
for 120 bins), which gives a p value of 0.69 (0.79). The
values of the A parameters associated with these plots(Fig. 21) also agree within statistical error. The HKKM11
simulation thus models the geomagnetic effects wellenough to be consistent with the current data.
Azimuthal distributions were also published by the
Bartol group in 2003 [21] for a limited energy and zenith
range; by reweighting the MC histograms in Fig. 20 to the
Bartol predictions we find χ 2Bartol − χ 2HKKM11 ¼ 1.0, whichis too low to indicate any preference between the models
from the azimuthal distributions alone.
The azimuthal distributions of the subsample ofevents selected to optimize the significance of the eastwest dipole asymmetry are shown in Fig. 22. The result is
an asymmetry parameter of Aμ ¼ 0.108  0.014ðstatÞ 
0.004ðsystÞ for μ-like events, and Ae ¼ 0.153 
0.015ðstatÞ  0.004ðsystÞ for e-like events. The dipole
effect is thus seen at a significance level of 6.0σ (8.0σ)
for the μ-like (e-like) samples, which represents thediscovery of the effect in the μ-like sample.
Figure 23 shows the calculated B parameters for data andMC simulations. The data correspond to a 2.2σ indication
for the existence of a zenith-dependent change in theasymmetry dipole alignment, and the agreement between

data and MC simulations provides further confidence in thetreatment of geomagnetic effects in the HKKM11 model.
While various detailed measurements have been madeconsidering geomagnetic effects on other cosmic-ray particles, such as cosmic-ray protons [80], this is the firstneutrino measurement that indicates agreement with a more
complicated model than a simple east-west effect based on
the geophysical dipole approximation.
V. SEARCH FOR A CORRELATIONBETWEEN THE SOLAR CYCLE
AND THE NEUTRINO FLUX
A. Effect predictions
Data are provided by the HKKM flux group [74] thatpredict the relative normalization change of the νe , ν̄e , νμ ,and ν̄μ fluxes at the SK site depending on the solar cycle.
These data use the count rate of a specific neutron monitor,
the Climax NM [81], as the parameter corresponding to thedegree of solar activity. This parameter can be used as themodel assumes that, after corrections for the local airpressure at the monitor, all NM counts are linearly inversely
correlated with the solar flux modulation. In fact, theplasma wind should affect the low-energy primary flux
more strongly; then, areas with a low rigidity cutoff, whichhave relatively more flux at low energy, should experience

052001-20

MEASUREMENTS OF THE ATMOSPHERIC NEUTRINO FLUX …

PHYSICAL REVIEW D 94, 052001 (2016)
80

e

700
650

60

600

40

parameter B [deg]

550
500
450
400
350
0

60

120

180

240

300

e

20
0
-20
-40

360

-60

µ

850

-80

800

1

750

0.6

0.2

0.6

0.2

-0.2

-0.6

-1

-0.2

-0.6

-1

µ

80

700

60
650

40

parameter B [deg]

600
550
500
0

60

120

180

240

300

Reconstructed Azimuth φ rec [deg]360

FIG. 22. Azimuthal distributions of a subselection of e-like(left) and μ-like (right) events, from the SK-I–SK-IV data (points
with statistical error) and MC simulations (boxes with systematic
error). The subselection is optimized to obtain the highest
significance of the final A parameters, by using only events with
0.4 < Erec < 3.0 (GeV) and j cos θrec j < 0.6.

a larger decrease in flux. However, if two NMs are locatedat different rigidity cutoffs, it is assumed that while thegradient of the correlation will be different, a linear
correlation is still applicable. The HKKM model calculates
the effects at the SK site for a given Climax NM count,
extrapolating the solar effect on neutrinos coming from all
directions, and considering that rigidity cutoff is a functionof the direction as was shown in Fig. 15.
As an example, Fig. 24 shows the expected normalization change for the νe flux at SK site, as a relative decrease
from the solar minimum to the solar maximum. Here, solar
minimum is defined as a Climax NM count of4150 counts hr−1 × 0.01. While higher values are occasionally recorded by the Climax NM, the change inneutrino flux as a function of NM count is expected tobe negligible above this value. The normalization change
weights are then provided in intervals of 50 down to3500 counts hr−1 × 0.01, which is defined as the solar
maximum and is the highest solar activity considered by
the HKKM model. The factor of 0.01 is the common
notation for NM data. From here on, we generally refer tothe solar modulation effect as a “relative normalization
change,” taking the expected flux at solar minimum as

20
0
-20
-40
-60
-80
1

cos

rec

FIG. 23. The parameter B depending on the cosine of thezenith angle, for e-like (left) and μ-like (right) events, from the
SK-I–SK-IV data (points with statistical error) and MC simulations (boxes with systematic error).

the default normalization. In Fig. 24, we see that only theneutrino flux up to around 1 GeV is strongly affected by theeffect. It is also seen that there is a strong zenith dependence of the normalization change, which is due to the factthat neutrinos originating at the north and south polar
regions reach the SK detector from below, and the polar
regions are the areas of lowest rigidity cutoff on Earth.
B. Neutrino and neutron monitor data
The solar modulation analysis selects the single-ring
sub-GeV samples, which have reconstructed energies
Erec < 1330 MeV as described in Sec. II B. Other samples
are not used, since multi-GeV samples show a negligible
solar modulation effect, and the sub-GeV multiring and PC
samples have relative normalization systematic errors
between each SK period that are larger than the expectedsolar modulation effect. This analysis uses an additional six
months of data in addition to the data set described in

052001-21

E. RICHARD et al.

PHYSICAL REVIEW D 94, 052001 (2016)
1

0.8

0.98

0.6

0.96

0.4

0.94

0.2

0.92

0

0.9

-0.2

0.88

-0.4

0.86

cos θ

1

-0.6

0.84

-0.8

0.82

-1
10 -1

1

10

10 2
E ν [GeV]

10

3

10 4

10

5

FIG. 24. The predicted normalization change for the νe flux
due to solar activity, as a function of the true neutrino energy Eν
and zenith angle θ, as given by the HKKM model. The z axis
shows the expected relative change in flux from the solar
activity minimum, defined as the activity corresponding to a
Climax NM count of 4150 counts hr−1 × 0.01, to solar maximum, defined as corresponding to a Climax NM count of3500 counts hr−1 × 0.01. The normalization changes for theν̄e , νμ , and ν̄μ fluxes have the same form, with only small
variations within the 2% level.

Sec. II C, up to April 2015, in order to cover more of thecurrent period of high solar activity. This gives 10 892 (10264.3) events in the e-like single-ring sub-GeV sample, and
10 763 (10 895.3) events in the μ-like single-ring sub-GeV
sample, where numbers in brackets are the MC expectations
including oscillations. The purity of the samples, defined as
the fraction of correct-flavor CC or NC interactions, isestimated by MC simulations as 95.0% and 98.2%,
respectively. Studies were done to optimize the selectedenergy ranges of these samples, in case higher sensitivity
could be obtained by using a subselection of the data at
lower energy where the solar modulation effect is stronger,
but no significant benefits in sensitivity were found.
In addition to being categorized as either e like or μ like,
we divide the data into upwardgoing or downwardgoingsamples based on the fitted lepton direction, resulting in a
total of four data subsamples. For each sample, the SK-I–
SK-IV periods will be combined while considering theirrelative normalization errors, as explained later in Sec. V C.
The same categorizations are applied for MC events.
For our neutron data, we cannot directly use data fromthe Climax NM, which was shut down during the SK-II
period. We searched for monitors [82] that were active andwell maintained throughout the entire SK-I–SK-IV period,and found four such NMs—the Thule, McMurdo,
Kerguelen, and Newark monitors [69,83]. These stations
all monitor downgoing neutrons with relatively low rigidity
cutoffs (0.3, 0.3, 1.14, and 2.4 GV respectively), and, thus,
have good sensitivity to the changes in the solar cycle. In
particular, the Thule and McMurdo monitors in the northand south polar regions are very sensitive to solar activity.

To obtain an “equivalent Climax NM count” for each ofthese NMs, we compare the counts of the monitors during
the times when they were both operational, as shown in
Fig. 25. While the correlation between each monitor isalmost linear, we find that a small improvement in theχ 2 =DOF can be obtained using a second-order polynomial
fit; we therefore use second-order polynomials as theconversion functions between monitors. This nonlinearity
is mostly due to the outlier events during particularly strongsolar activity, where the assumption of uniformity in theneighborhood of the Earth may be less accurate. We thendefine a “NM parameter” C as the average of the four NM
monitor values after conversion to the equivalent Climax
NM counts. This parameter is thus directly comparable to
the Climax NM count, as used by the HKKM model. The
combination of four monitors also reduces the influence ofpossible systematic shifts in a monitor; the systematic error
on C is estimated by taking the average rms of the fourconverted counts, and found to be 15.8 counts hr−1 × 0.01
at 1σ. The variance of C over the SK-I–SK-IV period isshown in Fig. 26, showing that, with recent data included,
almost two solar maxima are covered by the SK data,
although there is some downtime between SK periods.
C. Analysis method
The effect of the solar modulation on the four reconstructed data samples is calculated using the SK MC
simulations, by reweighting events based on their true
neutrino properties corresponding to the relative normalization change predicted by the HKKM model (as describedin Sec. VA), and is shown in Fig. 27. It can be seen that thenormalization change in the atmospheric neutrino flux is
not linear in the NM parameter C, and only has a significant
effect at NM count values corresponding to particularlyhigh solar activities. Comparing with Fig. 26, we see that
only the data obtained during the solar maxima at the end of
SK-I and beginning of SK-II, and some of the recent SK-IV
data, will be strongly affected by the solar modulation. Of
course, the data in the low-solar-activity periods are alsoessential, as we are searching for a relative decrease.
The four functions in Fig. 27 give a single prediction of
the relative normalization change in more accurate way
than simple linear fitting would achieve; however, it ispossible that in reality the effect is stronger or weaker thanthe prediction. To test for this possibility, we allow thesefunctions to be scaled by a continuous parameter α, where
α ¼ 0 would represent no correlation between the neutrino
flux and the NM count, and α ¼ 1 represents the default
predicted correlation. These relative normalization functions are denoted as f s ðC; αÞ, for each data sample s. As an
example, Fig. 28 shows f s ðC; αÞ in the case of the e-likeupgoing sample for various fixed values of α. This set offunctions will be fitted simultaneously to the four samples,in order to measure the overall strength of the solar
modulation effect in terms of the single parameter α.

052001-22

4400

Climax NM [counts hr -1 × 0.01]

Climax NM [counts hr -1 × 0.01]

MEASUREMENTS OF THE ATMOSPHERIC NEUTRINO FLUX …
p = 390.90
p = 0.649
1
p = 1.24×10-4

4200
4000

2

3800
3600
3400
3200
3000
2800
2600
2400

2600

2800

3000

3200

3400

PHYSICAL REVIEW D 94, 052001 (2016)
14400

p = -29540
p = 2.26
1
p = -1.49×10-4

4200
4000

0.8

3600
3400
3200

0.7

3000
2800

0.6

2600

3600

3000

3500

4200
4000

4000

4500

5000

Thule NM [counts hr -1 × 0.01]
Climax NM [counts hr -1 × 0.01]

Climax NM [counts hr -1 × 0.01]

Newark NM [counts hr -1 × 0.01]
4400

0.9

2

3800

p = 6.0250
p = 0.375
1
p = 5.43×10-6
2

3800
3600
3400
3200
3000
2800
2600
6500 7000 7500 8000 8500 9000 9500 10000 10500 11000

4400

0.5
0.4

p = -313.40
p = 0.528
1
p = 4.52×10-6

4200
4000

0.3

2

3800
3600

0.2

3400
3200

0.1

3000
2800

0

2600
5500

McMurdo NM [counts hr-1 × 0.01]

6000

6500

7000

7500

8000

8500

Kerguelen NM [counts hr-1 × 0.01]

FIG. 25. Correlation of the neutron monitor counts between the Newark (upper left), Thule (upper right), McMurdo (lower left), andKerguelen (lower right) monitors to the Climax NM. The average pressure-corrected counts per hour were plotted and fitted with asecond-order polynomial, which are defined as the conversion functions between those NMs. The figure shows each point binned into a
2D histogram, where the z axis represents the arbitrarily normalized number of points, the conversion functions as black lines, and pi are
the polynomial parameters.

Higher values α > 1 or lower values α < 0 are allowed in
the fit, although α < 0 would represent the unexpected casefor which the atmospheric neutrino flux increases duringhigh solar activity.

In terms of systematic errors on the SK data, as the solar
modulation effect is a relative normalization shift depending on the solar activity, any systematic error that affects theoverall normalization of all SK data, such as neutrino cross

SK-I

NM parameter C [counts hr -1 × 0.01]

4400

SK-II
SK-III
SK-IV

4200

4000

3800

3600

3400
1996

1998

2000

2002

2004

2006

2008

2010

2012

2014

2016

Year

FIG. 26. The “NM parameter” C over the SK operational period, where low values of C correspond to high solar activity, showing that
almost two solar cycles are covered by the SK data.

052001-23

E. RICHARD et al.

PHYSICAL REVIEW D 94, 052001 (2016)

1
0.99

Relative normalization

0.98
0.97
0.96
0.95

µ -like down-going

0.94µ -like up-going

0.93

e-like down-going

0.92

e-like up-going

0.91
0.9
3500

3600

3700

3800

3900

4000

4100

NM parameter C [counts hr -1 × 0.01]

FIG. 27. The relative normalization change on the SK-I–SK-IV
reconstructed data samples due to the solar modulation
effect, calculated using the SK MC simulations and HKKM
model prediction. The normalization change is defined relative tothe expected event rate corresponding to a NM parameterC ¼ 4150.

section and oscillation parameters, may be ignored.
However, we should consider detector reconstruction errorsthat are dependent on the detector period, because the solar
activity is a function of detector period; thus, we must allow
for systematic shifts when combining the SK-I–SK-IV
data. Such errors can arise from the detector changes suchas the difference in configuration of the PMT coverage or
the replacement of the electronics. Table VII shows theeffect of all SK-period-dependent systematic errors for ourselected data samples. The dominant contributions, at the>1% level, are the fiducial volume cut and the ring
separation uncertainty; for the e-like sample the singlering π 0 rejection uncertainty also contributes, and for the

μ-like sample the Michel electron tagging uncertainty
contributes, although this error becomes smaller in
SK-IV due to the improved electronics. Each error sourceis combined in quadrature, such that there are a total ofeight systematic errors, representing the normalization of
the e-like and μ-like samples in each of the SK-I–SK-IV
periods. The same error is assigned for both upgoing anddowngoing samples.
We next discuss the fitting method between the neutrino
and neutron data. First, each hourly period from January
1996 onwards is associated with a NM count, which is theaverage NM parameter C during that hour. Neutrino
events are then binned by the parameter C at theirobservation time, in 25 equally sized bins from 3300 to
4550 counts hr−1 × 0.01. We define our hypothesis for the
expected number of events in each bin i in each of the four
samples s as
Hs;i ðαÞ ¼ tobs
i × rs × f s;i ðαÞ;

ð5:1Þrepresents the total observation time by thewhere tobsi
SK detector at each value of C, rs is the nominal event rateof that sample, and f s;i ðαÞ is the solar modulation
function f s ðC; αÞ for that sample (an example of which
was shown in Fig. 28) evaluated at the central value of C
for the bin i.
Figure 29 shows the binning as a function of C, andcompares the expected event rate H of the e-like upgoingsample in the case of no solar modulation effect to the caseof the standard solar modulation effect. The minimum
value on the Climax NM parameter in these plots is~ −1 × 0.01; while there is a small amount3300 counts hrof exposure time below this value, these data will be treated
Expected events (e-like up-going sample)

3

1.4

Relative normalization

1.2
1
0.8

= -1
0.6

=0
=1

0.4

=2

0.2
0
3200

10

102

α= 0

10

α= 1

3400

3400

3600

3800

4000

4200

3600

3800

4000

4200

4400

NM parameter C [counts hr -1 × 0.01]4400

NM parameter C [counts hr -1× 0.01]

FIG. 28. An example of the relative normalization functions
fs ðC; αÞ, showing the case where sample s is the e-like upgoingsample, for various values of α. These functions are based on the
e-like upgoing function shown in Fig. 27, but extended below
C ¼ 3500 down to C ¼ 3300 by a polynomial fitting, and above
C ¼ 4150 by assuming a relative normalization of unity.

FIG. 29. Monte Carlo predictions for the number of events in
the SK-I–SK-IV data, showing two example hypotheses Hs;i ðαÞfor the cases of no solar activity correlation (α ¼ 0, black), and acorrelation as predicted by the HKKM model (α ¼ 1, blue, withstatistical error band). The blue histogram is thus essentially theblack histogram scaled by the function f s ðC; α ¼ 1Þ, shown in
Fig. 28. The sample s is the e-like upgoing sample.

052001-24

MEASUREMENTS OF THE ATMOSPHERIC NEUTRINO FLUX …

PHYSICAL REVIEW D 94, 052001 (2016)

TABLE VII. Contributions from each systematic error source on the uncertainty in the number of observed single-ring e-like and μlike events in the solar modulation analysis, at the 1σ level, separated by detector period. As described in the text, overall normalization
errors which affect all SK periods in the same way are not considered. The “Total” row represents the combination of all sources inquadrature.
Single-ring sub-GeV μ-like (%)

Single-ring sub-GeV e-like (%)
Systematic error source

SK-I

SK-II

SK-III

SK-IV

SK-I

SK-II

SK-IIISK-IV

FC reduction
FC/PC separation
Non-νe backgroundNon-νμ background
Ring separation
Single ring particle ID
Absolute energy calibration
Up/down energy calibration
Azimuthal energy calibration
π 0 rejection in e-like sample
μ → e decay tagging
Fiducial volume cut
Total

0.20

0.50

1.49
0.23
0.12
<0.01
0.01
1.12
<0.01
2.00
2.80

0.20

0.20

1.47
0.66
0.18
<0.01
0.01
1.12
<0.01
2.00
2.82

0.80

0.10

1.88
0.26
0.30
0.01
0.01
0.83
0.01
2.00
3.00

0.30

0.10

0.51
0.28
0.26
<0.01
0.01
1.00
<0.01
2.00
2.35

0.2
0.01

0.21
0.70
0.18
0.40
0.03
0.03

1.29
2.00
2.53

0.20
<0.01

0.11
1.33
0.50
0.61
0.03
0.01

1.35
2.00
2.87

0.80
0.01

0.11
1.56
0.19
0.99
0.06
0.02

1.32
2.00
3.13

0.30
<0.01

0.12
0.84
0.22
0.85
0.01
0.02

0.78
2.00
2.48

separately as discussed later. We considered the possibility
that the small systematic error in the NM parameter could
cause events to migrate between bins and affect the final
result. Using the size of the systematic, estimated inSec. V B, to randomly shift the NM parameter for eachhour of observation, such migration was shown to have anegligible impact on the final measurement of α.
From Fig. 29 we see that some bins will have low
statistics, and the Gaussian approximation for statisticalerrors will be poor. We thus use a Poisson likelihood-ratio
method to compare the α ¼ 0 and α ≠ 0 hypotheses. To
take account of the eight systematic errors as shown in
Table VII, we first consider that we can vary the strength oftheir effects compared to their default strength by writing avector of error pulls ϵ, referring to the strength of eachsystematic in units of σ. We then modify our predictions forthe number of events as
H0s;i ðks ; ϵ; αÞ ¼ H s;i ðαÞ × Ss;i ðϵÞ × ks ;

ð5:2Þwhere Ss;i is the systematic modification based on thevector of systematic error pulls, and ks is a factor allowing afree renormalization of each of the four data samples.
While the modification Ss;i simply shifts the relative
normalization of each SK period, the shape of the modification on the hypotheses Hs;i is somewhat complicated, as
the relative contribution of each SK period varies depending on each bin i. If we denote our data results by N s;i,the Poisson log-likelihood including systematic error isdefined asln LðNjHÞ ¼s;i −H 0s;i 
X X H0Ns;i eþ ϵ2 ;ln
N s;i !
si

ð5:3Þwhere the second term is a penalty term introduced for thesystematic errors, to avoid unphysically large changes of
the relative normalization in each SK period. The final
likelihood ratio statistic is then defined asΛ ¼ 2 lnLM ðNjH0 Þ;
LM ðNjH 0 ðα ¼ 0ÞÞ

ð5:4Þ

where the subscript M denotes the maximum likelihood
estimator, i.e., the best fit over all parameters of H0 . In the
denominator, however, α is explicitly fixed to be zero. By
Wilks’s theorem, as the difference in parameter space is 1,
the test statistic Λ is expected to be distributed as χ 2 with
one degree of freedom. The significancepﬃﬃﬃﬃ of rejecting the
α ¼ 0 hypothesis is then given by Λ.
Toy data sets are generated for each of the four samples,for both α ¼ 0 and α ¼ 1 hypotheses, in order to test thesensitivity to observe a long-term solar activity correlation.
Each toy was created by generating a random set ofsystematics ϵ according to a Gaussian distribution, calculating the associated hypothesis H0s;i , then generating eventsin each bin based on a Poisson distribution. One thousand
toys are generated and passed through the fitting procedure,and several tests are performed on the results to ensure thatthe test statistic is distributed as expected and thepfitting
ﬃﬃﬃﬃ
procedure is unbiased. By taking the significance Λ foreach toy, we expect on average a 1.75σ sensitivity toobserve a nonzero solar activity correlation.
D. Seasonal correlation analysis
The HKKM group has predicted that at the SK site, thewinter decrease in the neutrino flux normalization is muchless than 1% [56]. The effect is stronger at higher energies,

052001-25

E. RICHARD et al.

PHYSICAL REVIEW D 94, 052001 (2016)

but even at the TeV scale the change is only at the 1% level.
As the SK detector has accumulated limited statistics atsuch energies, there is essentially no sensitivity to measure
the predicted changes.
Nonetheless, we perform a simple search for a seasonal
effect, as a sanity check of our data and the model
predictions. For this test, all of the sub-GeV data samples
used in the solar correlation analysis are combined into one
sample; however, as the effect is predicted to exist onlyat higher energies, we also combine and test all of themulti-GeV samples as defined in Table I in Sec. II C,
including PC and UPMU data. The search is done byplotting the average daily event rate over each month of theyear, then performing a χ 2 comparison with the best-fit
constant event rate. We also use the unbinned data fromthis plot and perform a Kuiper test, which is a type ofKolmogorov-Smirnov test that is invariant under cyclic
transformations, i.e., it is not biased by the fact that we start
counting in January and is able to detect a deviation from a
constant event rate with equal significance across the year.
The comparison function for the Kuiper test is taken as theaverage event rate over all the data, with exact consideration of the live time accumulated throughout theyear, accounting for the various downtime in each SK
period.
Applying the predicted 11-year solar modulation effectto the data when binned by month is calculated to cause lessthan a 0.3% difference between bins, so is considerednegligible in this test; similarly, the effects of SK-perioddependent systematics are disregarded, as each period
contributes roughly equally across the year. Conversely,
assuming that no unexpectedly strong seasonal correlationis seen, the effects from the seasonal correlation on the solar
correlation study were also calculated to be negligible anddisregarded.
E. Extreme solar events analysis
Some sharp downwards peaks can be seen in Fig. 26; thecause of these events are usually large coronal massejections, causing relatively short and high-intensity
increases in the solar plasma flux [84]. On Earth these
are associated with particularly large decreases in thecosmic-ray flux and are often termed Forbush decreases[85]. While there is some ambiguity in the term, a common
definition is any event that causes a decrease greaterthan 10% in the intensity of cosmic rays on the Earth’s
surface.
So far, we have discussed a search for a long-termcorrelation between the neutrino flux and the solar
activity, using events down to a NM parameter C of
3300 counts hr−1 × 0.01. In fact, C sometimes drops
significantly below this value, although such observations
are rare and spread over a range of values of C. These
periods are listed in Table VIII, and correspond to Forbush
decrease events occuring during solar maxima periods

TABLE VIII. The periods for which the NM parameter C drops
below 3300 counts hr−1 × 0.01, corresponding to strong Forbush
decrease events. Start and end times are given in UTC.
Start time

End time

Hour

15 July 2000, 18:00
11 April 2001, 23:00
29 October 2003, 11:00
01 November 2003, 00:00
19 January 2005, 00:00
Total

17 July 2000, 21:00
13 April 2001, 14:00
01 November 2003, 00:00
04 November 2003, 13:00
19 January 2005, 13:00

50
38
61
67
13
229

(Forbush decreases occuring during solar minima, on theother hand, did not generally cause C to fall below
3300 counts hr−1 × 0.01).
These time periods were not used in the long-termcorrelation analysis, as the HKKM simulation of the solar
modulation effect only provides a prediction down to a
Climax NM count of 3500 counts hr−1 × 0.01. While we
extended the correlation down to a NM parameter C ¼
3300 counts hr−1 × 0.01 by a polynomial fitting, we
suspect that extending the correlation further would be
an unrealistic model. In particular, a simple correlation
between the neutron and neutrino fluxes at various locations on Earth may not hold during such extreme
events.
We thus analyze separately these time periods where C
falls below 3300 counts hr−1 × 0.01. While no theoreticalprediction of the fractional neutrino flux decrease isavailable here, these may be the most sensitive times tomeasure a solar effect on the atmospheric neutrino flux. We
define a second search method by making a simple test
against the null hypothesis, by counting the number of
events observed during all such periods over all four data
samples, and comparing with the nominal event rate, takenfrom the fitted normalization constants ks × rs of the bestfit hypotheses in the long-term analysis.
F. Results and discussions
The data for the long-term solar cycle correlation search,summing over SK-I–SK-IV and dividing into the four
sample types (e-like or μ-like and upgoing or downgoing
events) are shown in Fig. 30. The test statistic Λ depending
on α is shown in Fig. 31. The best-fit value of alpha is
α ¼ 0.62þ0.57−0.58 , with errors given at 1σ. Although thesensitivity to measure α is fairly low, the data in Fig. 30
show no unexpected disagreements with the model; thefour plots together have a χ 2 statistic of 88.0 for 100 bins.
The rejection power of the null hypothesis (α ¼ 0) is 1.06σ.
This is lower than the mean value predicted by toy MC
simulations with α ¼ 1, but still a reasonably likely result
according to that MC simulation, with a p value ofp ¼ 0.26.

052001-26

MEASUREMENTS OF THE ATMOSPHERIC NEUTRINO FLUX …
e-like up-going

1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0

3400

3600

3800

4000

4200

4400

1.4
1.2
1
0.8
0.6
0.4
0.2
3400

3600

3800

4000

4200

4400

NM parameter C [counts hr -1 × 0.01]

4600

2

µ -like up-going

1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
3400

3600

3800

4000

4200

4400

4600

Neutrino events [counts day -1]

2

Neutrino events [counts day -1]

e-like down-going

1.8
1.6

0

4600

NM parameter C [counts hr -1 × 0.01]

0PHYSICAL REVIEW D 94, 052001 (2016)

2

Neutrino events [counts day -1]

Neutrino events [counts day -1]

2µ -like down-going

1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0

3400

NM parameter C [counts hr -1 × 0.01]

3600

3800

4000

4200

4400

4600

NM parameter C [counts hr -1 × 0.01]

FIG. 30. The test for a solar modulation correlation using the SK-I–SK-IV data (points). The solar correlation hypotheses H0s;i ðαÞ, as
explained in Sec. V C, are shown for no correlation (α ¼ 0, grey), best fit (α ¼ 0.62, red), and the default prediction (α ¼ 1, grey dotted)
for each of the four data samples. For clearer visualization, the H 0s;i are actually shown without the systematic pulls Ss;i ðϵÞ, but
equivalent and opposite shifts are instead applied to the data points. The statistical error bars are also drawn on the data instead of themodel, using Pearson’s-χ 2 -based Poisson errors as defined in [86].

We also performed the same analysis on subsamples ofthe data, as shown in Fig. 32; we test the periods SK-I to
SK-IV separately, then each data subsample individually(while including all SK periods). Data from the SK-III
period alone cannot give any result, as it has no observation

time above the minimum solar activity required to cause
any effect according to the HKKM model. Although SK-II
and SK-IV prefer a low value of α, the statistical power islower and not inconsistent with the overall result.
Somewhat interestingly, the e-like samples prefer no
3

20

2

18
16

1

14

0

12

3

10

-1

8
6
4
2
0
-2

-1

0

1

2

-2

1

-3

2

3

FIG. 31. The test statistic Λ as a function of the solar
modulation strength parameter α, both of which are defined in
the text. The significance levels are drawn assuming the validityof Wilks’s theorem. This statistic represents the combined
measurement of all four data samples.

SK-I

SK-II

SK-IIIe up
SK-IV
Sub-sample

e downµ upµ downFIG. 32. The best-fit α parameter obtained by performing theanalysis on subsamples of the data, showing SK periods on theleft and data subsample types on the right. The red band shows
the combined result using all data. The SK-III data point isabsent, due to a lack of any data in the sensitive solar-active
region and, thus, an inability to perform any fit by itself.

052001-27

E. RICHARD et al.

PHYSICAL REVIEW D 94, 052001 (2016)

2.8
Multi-GeV
Neutrino events [counts day -1]

2.6
2.4
2.2
Sub-GeV
4.6
4.4
4.2
4
Jan

Feb Mar Apr May Jun

Jul

Aug

Sep

Oct Nov

Dec

Month

FIG. 33. The average daily event rates for each month, using the
SK-I–SK-IV data (points with statistical error only). The dottedlines are the best-fit constant functions.

correlation, while the μ-like samples prefer more theexpected α ¼ 1 correlation; however, the significance is
not high enough to draw any strong conclusions.
The search for a seasonal correlation is shown in Fig. 33.
The χ 2 test statistic, comparing against the best-fit constantfunctions, is 12.8 for the sub-GeV data and 6.5 for themulti-GeV data. The Kuiper test on the unbinned data gives
p values of 0.76 and 0.62 respectively. As expected, nostrong indication is seen for any seasonal correlation at theSK site.
During the coronal mass ejection periods described inTable VIII, the SK detector was operational for a total of
7.21 days. Using the fitted normalization constants ks × rs
from the long-term analysis, we can calculate the number of
events expected in the case of no solar modulation effect;
we find that 31.80  0.17 sub-GeV single-ring events are
expected. The actual number of events recorded by thedetector was nF ¼ 20, which by consulting the Poisson
distribution gives PðnF ≤ 20Þ ¼ 0.017, corresponding torejection of the null hypothesis of no solar cycle effect at
the 98.3% (2.38σ) significance level.
This significance level is higher than the long-termcorrelation search and, while not high enough to claim
direct evidence of a solar activity effect, is nonetheless anindication that it may be possible to measure such effects
with higher accuracy in the next-generation neutrino
detectors.
VI. SUMMARY
The energy spectra of the directionally integrated atmospheric νe and νμ particle plus antiparticle fluxes at SuperKamiokande were measured using a Bayesian iterative
unfolding method. A wide energy range was covered from
100 MeV up to 10 TeV, with higher accuracy than previous

measurements. The measured νμ energy spectrum is consistent with the IceCube and AMANDA results, which
cover the energy range above 100 GeV.
The results were also compared with Monte Carlopredictions, including the HKKM11, FLUKA, andBartol models, by performing a χ 2 test incorporating boththe νe and νμ data. The models were found to be consistent
with the data. The flux normalizations and spectral inclinations of each model were also tested and found to beconsistent.
The azimuthal spectra of the fluxes were also comparedwith the Monte Carlo predictions, which is the first detailed
check of the azimuthal spectra including a full systematic
error analysis. The data and MC agreed well, which
provides confidence in the careful treatment of geomagnetic effects in the recent flux simulations. The existence ofan east-west dipole asymmetry in the flux was alsomeasured with greatly increased statistics compared toprevious measurements. The effect was seen at a significance of 6σ for νμ and 8σ for νe, which represents discoveryof the effect in the νμ flux. An indication that the angle ofthe dipole asymmetry shifts depending on the zenith angle
was found at the 2.2σ level, which is the first measurementthat explores the geomagnetic field effects beyond thesimple east-west asymmetry.
The expected correlation between the atmospheric neutrino flux and the solar magnetic activity was studied, using
the SK sub-GeV data sample and the ground-level neutron
flux measured at various neutron monitors. The predictedeffects on the neutrino flux based on the HKKM model
were found to be relatively small in the SK data set, but by
using data spanning 20 years a weak preference for somecorrelation was seen at the 1.1σ level. We also separately
examined several short periods of particularly intense solar
activity, for which no theoretical prediction is available,
corresponding to 7.1 days of detector uptime. During thisperiod, an indication for a decrease in the neutrino event
rate below the normal level was seen at the 2.4σ level.
The seasonal change in the neutrino flux was also
examined. While some seasonal correlation is expected
due to changes in the atmospheric density profile over theyear, the effect is predicted to be negligible at the SK siteand, as expected, no such correlation was seen.
ACKNOWLEDGMENTS
The authors would like to thank M. Honda for manysuggestions and discussions, and for providing the data on
the correlation between the atmospheric flux and neutron
count. The authors would like to thank S. Yoshida, K.
Mase, A. Ishihara, T. Kuwabara, and M. Schmitz for theuseful discussions on the atmospheric neutrino flux measurement. We acknowledge the NMDB database [82],founded under the European Union’s FP7 programme(Contract No. 213007) for providing data. The Newark,Thule, and McMurdo neutron monitors of the Bartol

052001-28

MEASUREMENTS OF THE ATMOSPHERIC NEUTRINO FLUX …

Research Institute are supported by the National Science
Foundation. Kerguelen neutron monitor data were kindly
provided by the French Polar Institute (IPEV, Brest) and by
the Paris Observatory. The authors gratefully acknowledge
the cooperation of the Kamioka Mining and Smelting
Company. Super-Kamiokande has been built and operated
from funds provided by the Japanese Ministry of
Education, Culture, Sports, Science, and Technology, theU.S. Department of Energy, and the U.S. National Science
Foundation. This work was supported by Grants-in-Aid for
Scientific Research on Innovative Areas, Grants
No. 25105004 and No. 25105005. This work was partially
supported by the joint research program of the Institute for
Cosmic Ray Research, at the University of Tokyo. This
work was partially supported by the Research Foundation
TABLE IX.

PHYSICAL REVIEW D 94, 052001 (2016)of Korea (BK21 and KNRC), the Korean Ministry of
Science and Technology, the National Science Foundationof China, the European Union H2020 Grant No. RISEGA641540-SKPLUS, the National Science andEngineering Research Council (NSERC) of Canada, theScinet and Westgrid consortia of Compute Canada, and the
National Science Centre, Poland (Grants No. 2015/17/N/
ST2/04064 and No. 2015/18/E/ST200758).

APPENDIX: CORRELATION MATRIXOF THE νe AND νμ ENERGY SPECTRUM
Please see Tables IX and X for the correlation matrix ofthe νe and νμ energy spectrum.

Values of correlation matrix element (Cˆij ) in the νe energy spectrum measurement.
νelog10 ðEν ½GeVÞ−0.8 to−0.6−0.6 to−0.4−0.4 to−0.2−0.2 to0.0

0.0 to0.2

0.2 to0.4

0.4 to0.6

0.6 to0.8

0.8 to1.0

1.0 to1.5

1.5 to2.0−0.8 to −0.6
−0.6 to −0.4
−0.4 to −0.2
−0.2 to 0.0
0.0 to 0.2

1.000
0.948
0.720
0.512
0.396

0.948
1.000
0.880
0.678
0.580

0.720
0.880
1.000
0.914
0.836

0.512
0.678
0.914
1.000
0.931

0.396
0.580
0.836
0.931
1.000

0.396
0.557
0.802
0.915
0.903

0.267
0.400
0.631
0.776
0.779

0.156
0.288
0.505
0.633
0.698

0.120
0.240
0.443
0.566
0.625

−0.079−0.034
0.064
0.131
0.187

−0.130−0.124−0.084−0.048−0.005νe

0.2 to 0.4
0.4 to 0.6
0.6 to 0.8
0.8 to 1.0
1.0 to 1.5
1.5 to 2.0−0.6 to −0.4
−0.4 to −0.2
−0.2 to 0.0
0.0 to 0.2
0.2 to 0.4
0.4 to 0.6

0.396
0.267
0.156
0.120
−0.079−0.130
0.843
0.692
0.520
0.423
0.400
0.291

0.557
0.400
0.288
0.240
−0.034−0.124
0.888
0.832
0.686
0.597
0.544
0.426

0.802
0.631
0.505
0.443
0.064
−0.084
0.809
0.940
0.896
0.839
0.771
0.659

0.915
0.776
0.633
0.566
0.131
−0.048
0.682
0.899
0.944
0.915
0.889
0.798

0.903
0.779
0.698
0.625
0.187
−0.005
0.537
0.842
0.900
0.941
0.874
0.824

1.000
0.909
0.753
0.700
0.195
−0.024
0.592
0.812
0.886
0.900
0.946
0.907

0.909
1.000
0.917
0.820
0.293
0.044
0.458
0.655
0.732
0.790
0.903
0.945

0.753
0.917
1.000
0.880
0.359
0.101
0.307
0.533
0.594
0.698
0.776
0.873

0.700
0.820
0.880
1.000
0.462
0.130
0.260
0.467
0.527
0.624
0.700
0.795

0.195
0.293
0.359
0.462
1.000
0.927
−0.019
0.085
0.100
0.169
0.172
0.227

−0.024
0.044
0.101
0.130
0.927
1.000
−0.110−0.069−0.072−0.030−0.052−0.024νμ0.6 to 0.8
0.8 to 1.0
1.0 to 1.5
1.5 to 2.0
2.0 to 3.0
3.0 to 4.0

0.252
0.240
0.181
0.069
0.034
0.027

0.380
0.360
0.282
0.106
0.051
0.026

0.603
0.565
0.445
0.167
0.077
0.031

0.732
0.681
0.531
0.200
0.092
0.040

0.774
0.715
0.555
0.207
0.090
0.037

0.849
0.793
0.589
0.215
0.099
0.048

0.929
0.888
0.615
0.230
0.121
0.084

0.906
0.891
0.615
0.260
0.163
0.132

0.837
0.849
0.660
0.368
0.285
0.254

0.281
0.409
0.687
0.763
0.747
0.730

0.021
0.153
0.509
0.692
0.700
0.705

052001-29

E. RICHARD et al.
TABLE X.

PHYSICAL REVIEW D 94, 052001 (2016)

Values of correlation matrix element (Cˆij ) in the νμ energy spectrum measurement.
νμlog10 ðEν ½GeVÞ−0.6 to−0.4−0.4 to−0.2−0.2 to0.0

0.0 to0.2

0.2 to0.4

0.4 to0.6

0.6 to0.8

0.8 to1.0

1.0 to1.5

1.5 to2.0

2.0 to3.0

3.0 to4.0−0.8 to −0.6
−0.6 to −0.4
−0.4 to −0.2
−0.2 to 0.0
0.0 to 0.2

0.843
0.888
0.809
0.682
0.537

0.692
0.832
0.940
0.899
0.842

0.520
0.686
0.896
0.944
0.900

0.423
0.597
0.839
0.915
0.941

0.400
0.544
0.771
0.889
0.874

0.291
0.426
0.659
0.798
0.824

0.252
0.380
0.603
0.732
0.774

0.240
0.360
0.565
0.681
0.715

0.181
0.282
0.445
0.531
0.555

0.069
0.106
0.167
0.200
0.207

0.034
0.051
0.077
0.092
0.090

0.027
0.026
0.031
0.040
0.037

νe

0.2 to 0.4
0.4 to 0.6
0.6 to 0.8
0.8 to 1.0
1.0 to 1.5
1.5 to 2.0−0.6 to −0.4
−0.4 to −0.2
−0.2 to 0.0
0.0 to 0.2
0.2 to 0.4
0.4 to 0.6

0.592
0.458
0.307
0.260
−0.019−0.110
1.000
0.840
0.699
0.574
0.596
0.462

0.812
0.655
0.533
0.467
0.085
−0.069
0.840
1.000
0.934
0.872
0.816
0.702

0.886
0.732
0.594
0.527
0.100
−0.072
0.699
0.934
1.000
0.942
0.889
0.788

0.900
0.790
0.698
0.624
0.169
−0.030
0.574
0.872
0.942
1.000
0.912
0.853

0.946
0.903
0.776
0.700
0.172
−0.052
0.596
0.816
0.889
0.912
1.000
0.952

0.907
0.945
0.873
0.795
0.227
−0.024
0.462
0.702
0.788
0.853
0.952
1.000

0.849
0.929
0.906
0.837
0.281
0.021
0.388
0.637
0.720
0.807
0.885
0.970

0.793
0.888
0.891
0.849
0.409
0.153
0.369
0.599
0.674
0.753
0.836
0.907

0.589
0.615
0.615
0.660
0.687
0.509
0.293
0.484
0.543
0.593
0.626
0.637

0.215
0.230
0.260
0.368
0.763
0.692
0.096
0.174
0.196
0.220
0.226
0.226

0.099
0.121
0.163
0.285
0.747
0.700
0.039
0.074
0.086
0.100
0.107
0.111

0.048
0.084
0.132
0.254
0.730
0.705
0.015
0.026
0.027
0.039
0.048
0.058

νμ

0.6 to 0.8
0.8 to 1.0
1.0 to 1.5
1.5 to 2.0
2.0 to 3.0
3.0 to 4.0

0.388
0.369
0.293
0.096
0.039
0.015

0.637
0.599
0.484
0.174
0.074
0.026

0.720
0.674
0.543
0.196
0.086
0.027

0.807
0.753
0.593
0.220
0.100
0.039

0.885
0.836
0.626
0.226
0.107
0.048

0.970
0.907
0.637
0.226
0.111
0.058

1.000
0.956
0.658
0.258
0.151
0.111

0.956
1.000
0.787
0.415
0.310
0.278

0.658
0.787
1.000
0.845
0.756
0.694

0.258
0.415
0.845
1.000
0.983
0.927

0.151
0.310
0.756
0.983
1.000
0.955

0.111
0.278
0.694
0.927
0.955
1.000

[1] T. K. Gaisser and M. Honda, Annu. Rev. Nucl. Part. Sci. 52,
153 (2002).
[2] C. V. Achar et al., Phys. Lett. 18, 196 (1965).
[3] F. Reines, M. F. Crouch, T. L. Jenkins, W. R. Kropp, H. S.
Gurr, G. R. Smith, J. P. F. Sellschop, and B. Meyer, Phys.
Rev. Lett. 15, 429 (1965).
[4] Y. Fukuda et al. (Super-Kamiokande Collaboration), Phys.
Rev. Lett. 81, 1562 (1998).
[5] J. Hosaka et al. (Super-Kamiokande Collaboration), Phys.
Rev. D 73, 112001 (2006).
[6] B. Aharmim et al. (SNO Collaboration), Phys. Rev. C 72,
055502 (2005).
[7] Y. Ashie et al. (Super-Kamiokande Collaboration), Phys.
Rev. D 71, 112005 (2005).
[8] R. Wendell et al. (Super-Kamiokande Collaboration), Phys.
Rev. D 81, 092004 (2010).
[9] S. Abe et al. (KamLAND Collaboration), Phys. Rev. Lett.
100, 221803 (2008).
[10] F. An et al. (DAYA-BAY Collaboration), Phys. Rev. Lett.
108, 171803 (2012).
[11] J. Ahn et al. (RENO Collaboration), Phys. Rev. Lett. 108,
191802 (2012).

[12] Y. Abe et al. (DOUBLE-CHOOZ Collaboration), Phys.
Rev. Lett. 108, 131801 (2012).
[13] M. Ahn et al. (K2K Collaboration), Phys. Rev. D 74,
072003 (2006).
[14] P. Adamson et al. (MINOS Collaboration), Phys. Rev. Lett.
106, 181801 (2011).
[15] K. Abe et al. (T2K Collaboration), Phys. Rev. Lett. 112,
181801 (2014).
[16] K. Abe et al. (T2K Collaboration), Phys. Rev. Lett. 112,061802 (2014).
[17] Z. Maki, M. Nakagawa, and S. Sakata, Prog. Theor. Phys.
28, 870 (1962).
[18] B. Pontecorvo, Sov. Phys. JETP 26, 984 (1968).
[19] M. Honda, T. Kajita, K. Kasahara, S. Midorikawa,and T. Sanuki, Phys. Rev. D 75, 043006 (2007).
[20] M. Honda, T. Kajita, K. Kasahara, and S. Midorikawa, Phys.
Rev. D 83, 123001 (2011).
[21] G. Barr, T. Gaisser, P. Lipari, S. Robbins, and T. Stanev,
Phys. Rev. D 70, 023006 (2004).
[22] G. Battistoni, A. Ferrari, T. Montaruli, and P. Sala,
Astropart. Phys. 19, 269 (2003).

052001-30

MEASUREMENTS OF THE ATMOSPHERIC NEUTRINO FLUX …
[23] T. Sanuki, M. Motoki, H. Matsumoto, E. Seo, J. Wang et al.,Astrophys. J. 545, 1135 (2000).
[24] S. Haino, T. Sanuki, K. Abe, K. Anraku, Y. Asaoka et al.,Phys. Lett. B 594, 35 (2004).
[25] J. Alcaraz et al. (AMS Collaboration), Phys. Lett. B 490, 27(2000).
[26] J. Alcaraz et al. (AMS Collaboration), Phys. Lett. B 494,
193 (2000).
[27] S. Roesler, R. Engel, and J. Ranft, in Advanced Monte Carlofor Radiation Physics, Particle Transport Simulation, andApplications. Proceedings, Conference, MC2000, Lisbon,
Portugal, 2000 (2000), pp. 1033–1038.
[28] K. Niita, T. Sato, H. Iwase, H. Nose, H. Nakashima, and L.
Sihver, Radiation Measurements 41, 1080 (2006).
[29] M. Catanesi et al. (HARP Collaboration), Astropart. Phys.
30, 124 (2008).
[30] V. Agrawal, T. K. Gaisser, P. Lipari, and T. Stanev, Phys.
Rev. D 53, 1314 (1996).
[31] R. Engel et al., in The flux of atmospheric muons(Copernicus Gesellschaft, 2001), Vol. 3, p. 1029.
[32] A. Ferrari and P. R. Sala, in Proceedings of the Workshop on
Nuclear Reaction Data and Nuclear Reactors Physics,
Design and Safety, ICTP, Miramare-Trieste, Italy, 1996,
edited by A. Gandini and G. Reffo (World Scientific,
Singapore, 1998), Vol. 2, p. 424 (ATLAS Report
No. ATL-PHYS-97-113, 1997).
[33] M. Honda, in Proceedings of the 14th Workshop on Elastic
and Diffractive Scattering (EDS Blois Workshop) on
Frontiers of QCD: From Puzzles to Discoveries (2011),
Vol. C111215.
[34] G. Barr, T. Gaisser, S. Robbins, and T. Stanev, Phys. Rev. D
74, 094009 (2006).
[35] K. A. Olive et al. (Particle Data Group), Chin. Phys. C 38,
090001 (2014).
[36] H. Nishino et al. (Super-Kamiokande Collaboration), Phys.
Rev. Lett. 102, 141801 (2009).
[37] K. Abe et al. (Super-Kamiokande Collaboration), Phys.
Rev. D 90, 072005 (2014).
[38] J. Billard, L. Strigari, and E. Figueroa-Feliciano, Phys. Rev.
D 89, 023524 (2014).
[39] K. Daum (Frejus Collaboration), Z. Phys. C 66, 417 (1995).
[40] R. Abbasi et al. (IceCube Collaboration), Phys. Rev. D 79,
102005 (2009).
[41] R. Abbasi et al. (IceCube Collaboration), Astropart. Phys.
34, 48 (2010).
[42] R. Abbasi et al. (IceCube Collaboration), Phys. Rev. D 83,
012001 (2011).
[43] R. Abbasi et al. (IceCube Collaboration), Phys. Rev. D 84,
082001 (2011).
[44] M. Aartsen et al. (IceCube Collaboration), Phys. Rev. Lett.
110, 151105 (2013).
[45] M. G. Aartsen et al. (IceCube Collaboration), Phys. Rev. D
91, 122004 (2015).
[46] M. Aartsen et al. (IceCube Collaboration), Phys. Rev. Lett.
111, 021103 (2013).
[47] M. Aartsen et al. (IceCube Collaboration), Science 342,
1242856 (2013).
[48] M. G. Aartsen et al. (IceCube Collaboration), Phys. Rev.
Lett. 113, 101101 (2014).

PHYSICAL REVIEW D 94, 052001 (2016)[49] B. Rossi, Phys. Rev. 36, 606 (1930).
[50] T. Futagami et al. (Super-Kamiokande Collaboration), Phys.
Rev. Lett. 82, 5194 (1999).
[51] M. Honda, T. Kajita, K. Kasahara, and S. Midorikawa, Phys.
Rev. D 70, 043008 (2004).
[52] G. Battistoni, A. Ferrari, P. Lipari, T. Montaruli, P. Sala, andT. Rancati, Astropart. Phys. 12, 315 (2000).
[53] P. Lipari, Astropart. Phys. 14, 153 (2000).
[54] M. S. Potgieter, Living Rev. Solar Phys. 10, 3 (2013).
[55] J. Simpson, Space Sci. Rev. 93, 11 (2000).
[56] M. Honda, M. S. Athar, T. Kajita, K. Kasahara, and S.
Midorikawa, Phys. Rev. D 92, 023004 (2015).
[57] Y. Fukuda et al. (Super-Kamiokande Collaboration),
Nucl. Instrum. Methods Phys. Res., Sect. A 501, 418(2003).
[58] H. Noumi et al., Nucl. Instrum. Methods Phys. Res., Sect. A
398, 399 (1997).
[59] 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).
[60] E. Richard, Ph.D. thesis, University of Tokyo, 2015.
[61] K. Abe et al. (Super-Kamiokande Collaboration), Phys.
Rev. Lett. 107, 241801 (2011).
[62] K. P. Lee, Ph.D. thesis, University of Tokyo, 2012.
[63] T. Irvine, Ph.D. thesis, University of Tokyo, 2014.
[64] T. Barszczak, Ph.D. thesis, University of California, Irvine,
2005.
[65] Y. Ashie et al. (Super-Kamiokande Collaboration), Phys.
Rev. Lett. 93, 101801 (2004).
[66] S. Desai et al. (Super-Kamiokande Collaboration),
Astropart. Phys. 29, 42 (2008).
[67] Y. Hayato, Acta Phys. Pol. B 40, 2477 (2009).
[68] K. Abe et al. (Super-Kamiokande Collaboration), Phys.
Rev. D 91, 052019 (2015).
[69] http://neutronm.bartol.udel.edu/.
[70] G. D’Agostini, Nucl. Instrum. Methods Phys. Res., Sect. A
362, 487 (1995).
[71] G. D’Agostini, arXiv:1010.0632.
[72] T. Adye, arXiv:1105.1160.
[73] G. Lafferty and T. Wyatt, Nucl. Instrum. Methods Phys.
Res., Sect. A 355, 541 (1995).
[74] M. Honda (private communication).
[75] S. Adrian-Martinez et al. (ANTARES Collaboration), Eur.
Phys. J. C 73, 2606 (2013); 73, 2606 (2013).
[76] http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html.
[77] R. T. Merrill, The Magnetic Field of the Earth (Academic
Press, 1998).
[78] S. S. Wilks, Ann. Math. Stat. 9, 60 (1938).
[79] N. Gagunashvili, Proc. Sci., ACAT2007 (2007) 060.
[80] J. Alcaraz et al., Phys. Lett. B 472, 215 (2000).
[81] http://cr0.izmiran.rssi.ru/clmx/main.htm.
[82] http://www.nmdb.eu/.
[83] http://cr0.izmiran.rssi.ru/kerg/baseinfo.htm.
[84] A. Babu et al., arXiv:1304.5343.
[85] S. E. Forbush, Proc. Natl. Acad. Sci. U.S.A. 43, 28(1957).
[86] J. Heinrich, Report No. CDF/MEMO/STATISTICS/
PUBLIC/6438, 2003.

052001-31