This is the accepted manuscript made available via CHORUS. The article has been
published as:

Diffuse supernova neutrino background search at SuperKamiokande
K. Abe et al. (Super-Kamiokande Collaboration)
Phys. Rev. D 104, 122002 — Published 10 December 2021
DOI: 10.1103/PhysRevD.104.122002

Diffuse Supernova Neutrino Background Search at Super-Kamiokande
K. Abe,1, 43 C. Bronner,1 Y. Hayato,1, 43 K. Hiraide,1 M. Ikeda,1 S. Imaizumi,1 J. Kameda,1, 43 Y. Kanemura,1
Y. Kataoka,1 S. Miki,1 M. Miura,1, 43 S. Moriyama,1, 43 Y. Nagao,1 M. Nakahata,1, 43 S. Nakayama,1, 43 T. Okada,1
K. Okamoto,1 A. Orii,1 G. Pronost,1 H. Sekiya,1, 43 M. Shiozawa,1, 43 Y. Sonoda,1 Y. Suzuki,1 A. Takeda,1, 43
Y. Takemoto,1 A. Takenaka,1 H. Tanaka,1 S. Watanabe,1 T. Yano,1 S. Han,2 T. Kajita,2, 43 K. Okumura,2, 43
T. Tashiro,2 J. Xia,2 G. D. Megias,3 D. Bravo-Berguño,4 L. Labarga,4 Ll. Marti,4 B. Zaldivar,4 B. W. Pointon,6, 47
F. d. M. Blaszczyk,5 E. Kearns,5, 43 J. L. Raaf,5 J. L. Stone,5, 43 L. Wan,5 T. Wester,5 J. Bian,7 N. J. Griskevich,7
W. R. Kropp,7 S. Locke,7 S. Mine,7 M. B. Smy,7, 43 H. W. Sobel,7, 43 V. Takhistov,7, 43 J. Hill,8 J. Y. Kim,9I. T. Lim,9 R. G. Park,9 B. Bodur,10 K. Scholberg,10, 43 C. W. Walter,10, 43 S. Cao,53 L. Bernard,11 A. Coffani,11O. Drapier,11 S. El Hedri,11, ∗ A. Giampaolo,11, † M. Gonin,11 Th. A. Mueller,11 P. Paganini,11 B. Quilain,11
T. Ishizuka,12 T. Nakamura,13 J. S. Jang,14 J. G. Learned,15 L. H. V. Anthony,16 D. Martin,16 M. Scott,16
A. A. Sztuc,16 Y. Uchida,16 V. Berardi,17 M. G. Catanesi,17 E. Radicioni,17 N. F. Calabria,18 L. N. Machado,18
G. De Rosa,18 G. Collazuol,19 F. Iacob,19 M. Lamoureux,19 M. Mattiazzi,19 N. Ospina,19 L. Ludovici,20
Y. Maekawa,21 Y. Nishimura,21 M. Friend,23 T. Hasegawa,23 T. Ishida,23 T. Kobayashi,23 M. Jakkapu,23
T. Matsubara,23 T. Nakadaira,23 K. Nakamura,23, 43 Y. Oyama,23 K. Sakashita,23 T. Sekiguchi,23 T. Tsukamoto,23
Y. Kotsar,24 Y. Nakano,24 H. Ozaki,24 T. Shiozawa,24 A. T. Suzuki,24 Y. Takeuchi,24, 43 S. Yamamoto,24 A. Ali,25
Y. Ashida,25, ‡ J. Feng,25 S. Hirota,25 T. Kikawa,25 M. Mori,25 T. Nakaya,25, 43 R. A. Wendell,25, 43 K. Yasutome,25
P. Fernandez,26 N. McCauley,26 P. Mehta,26 K. M. Tsui,26 Y. Fukuda,27 Y. Itow,28, 29 H. Menjo,28 T. Niwa,28
K. Sato,28 M. Tsukada,28 J. Lagoda,30 S. M. Lakshmi,30 P. Mijakowski,30 J. Zalipska,30 J. Jiang,31 C. K. Jung,31
C. Vilela,31 M. J. Wilking,31 C. Yanagisawa,31, § K. Hagiwara,32 M. Harada,32 T. Horai,32 H. Ishino,32 S. Ito,32
H. Kitagawa,32 Y. Koshio,32, 43 W. Ma,32 N. Piplani,32 S. Sakai,32 G. Barr,33 D. Barrow,33 L. Cook,33, 43
A. Goldsack,33, 43 S. Samani,33 D. Wark,33, 38 F. Nova,34 T. Boschi,22 F. Di Lodovico,22 J. Gao,22 J. Migenda,22
M. Taani,22 S. Zsoldos,22 J. Y. Yang,35 S. J. Jenkins,36 M. Malek,36 J. M. McElwee,36 O. Stone,36 M. D. Thiesse,36
L. F. Thompson,36 H. Okazawa,37 S. B. Kim,39 J. W. Seo,39 I. Yu,39 K. Nishijima,40 M. Koshiba,41, ¶ K. Iwamoto,42
K. Nakagiri,42 Y. Nakajima,42, 43 N. Ogawa,42 M. Yokoyama,42, 43 K. Martens,43 M. R. Vagins,43, 7 M. Kuze,44
S. Izumiyama,44 T. Yoshida,44 M. Inomoto,45 M. Ishitsuka,45 H. Ito,45 T. Kinoshita,45 R. Matsumoto,45
K. Ohta,45 M. Shinoki,45 T. Suganuma,45 A. K. Ichikawa,54 K. Nakamura,54 J. F. Martin,46 H. A. Tanaka,46
T. Towstego,46 R. Akutsu,47 V. Gousy-Leblanc,47, k M. Hartz,47 A. Konaka,47, k P. de Perio,47 N. W. Prouse,47
S. Chen,48 B. D. Xu,48 Y. Zhang,48 M. Posiadala-Zezula,49 D. Hadley,50 M. O’Flaherty,50 B. Richards,50
B. Jamieson,51 J. Walker,51 A. Minamino,52 K. Okamoto,52 G. Pintaudi,52 S. Sano,52 and R. Sasaki52(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, Japan
3
Institute for Cosmic Ray Research, University of Tokyo, Kashiwa, Chiba 277-8582, Japan
4
Department of Theoretical Physics, University Autonoma Madrid, 28049 Madrid, Spain
5
Department of Physics, Boston University, Boston, MA 02215, USA
6
Department of Physics, British Columbia Institute of Technology, Burnaby, BC, V5G 3H2, Canada
7
Department of Physics and Astronomy, University of California, Irvine, Irvine, CA 92697-4575, USA
8
Department of Physics, California State University, Dominguez Hills, Carson, CA 90747, USA
9
Institute for Universe and Elementary Particles, Chonnam National University, Gwangju 61186, Korea
10
Department of Physics, Duke University, Durham NC 27708, USA
11
Ecole Polytechnique, IN2P3-CNRS, Laboratoire Leprince-Ringuet, F-91120 Palaiseau, France
12
Junior College, Fukuoka Institute of Technology, Fukuoka, Fukuoka 811-0295, Japan
13
Department of Physics, Gifu University, Gifu, Gifu 501-1193, Japan
14
GIST College, Gwangju Institute of Science and Technology, Gwangju 500-712, Korea
15
Department of Physics and Astronomy, University of Hawaii, Honolulu, HI 96822, USA
16
Department of Physics, Imperial College London , London, SW7 2AZ, United Kingdom
17
Dipartimento Interuniversitario di Fisica, INFN Sezione di Bari and Università e Politecnico di Bari, I-70125, Bari, Italy
18
Dipartimento di Fisica, INFN Sezione di Napoli and Università di Napoli, I-80126, Napoli, Italy
19
Dipartimento di Fisica, INFN Sezione di Padova and Università di Padova, I-35131, Padova, Italy
20
INFN Sezione di Roma and Università di Roma “La Sapienza”, I-00185, Roma, Italy
21
Department of Physics, Keio University, Yokohama, Kanagawa, 223-8522, Japan
22
Department of Physics, King’s College London, London, WC2R 2LS, UK
23
High Energy Accelerator Research Organization (KEK), Tsukuba, Ibaraki 305-0801, Japan
24
Department of Physics, Kobe University, Kobe, Hyogo 657-8501, Japan2
25

Department of Physics, Kyoto University, Kyoto, Kyoto 606-8502, Japan
Department of Physics, University of Liverpool, Liverpool, L69 7ZE, United Kingdom
27
Department of Physics, Miyagi University of Education, Sendai, Miyagi 980-0845, Japan
28
Institute for Space-Earth Environmental Research, Nagoya University, Nagoya, Aichi 464-8602, Japan
29
Kobayashi-Maskawa Institute for the Origin of Particles and the Universe, Nagoya University, Nagoya, Aichi 464-8602, Japan30
National Centre For Nuclear Research, 02-093 Warsaw, Poland
31
Department of Physics and Astronomy, State University of New York at Stony Brook, NY 11794-3800, USA
32
Department of Physics, Okayama University, Okayama, Okayama 700-8530, Japan
33
Department of Physics, Oxford University, Oxford, OX1 3PU, United Kingdom
34
Rutherford Appleton Laboratory, Harwell, Oxford, OX11 0QX, UK
35
Department of Physics, Seoul National University, Seoul 151-742, Korea
36
Department of Physics and Astronomy, University of Sheffield, S3 7RH, Sheffield, United Kingdom
37
Department of Informatics in Social Welfare, Shizuoka University of Welfare, Yaizu, Shizuoka, 425-8611, Japan
38
STFC, Rutherford Appleton Laboratory, Harwell Oxford, andDaresbury Laboratory, Warrington, OX11 0QX, United Kingdom
39
Department of Physics, Sungkyunkwan University, Suwon 440-746, Korea
40
Department of Physics, Tokai University, Hiratsuka, Kanagawa 259-1292, Japan
41
The University of Tokyo, Bunkyo, Tokyo 113-0033, Japan
42
Department of Physics, University of Tokyo, Bunkyo, Tokyo 113-0033, Japan
43
Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University ofTokyo Institutes for Advanced Study, University of Tokyo, Kashiwa, Chiba 277-8583, Japan
44
Department of Physics,Tokyo Institute of Technology, Meguro, Tokyo 152-8551, Japan
45
Department of Physics, Faculty of Science and Technology, Tokyo University of Science, Noda, Chiba 278-8510, Japan46
Department of Physics, University of Toronto, ON, M5S 1A7, Canada
47
TRIUMF, 4004 Wesbrook Mall, Vancouver, BC, V6T2A3, Canada
48
Department of Engineering Physics, Tsinghua University, Beijing, 100084, China
49
Faculty of Physics, University of Warsaw, Warsaw, 02-093, Poland
50
Department of Physics, University of Warwick, Coventry, CV4 7AL, UK
51
Department of Physics, University of Winnipeg, MB R3J 3L8, Canada
52
Department of Physics, Yokohama National University, Yokohama, Kanagawa, 240-8501, Japan
53
Institute For Interdisciplinary Research in Science and Education, ICISE, Quy Nhon, 55121, Vietnam
54
Department of Physics, Faculty of Science, Tohoku University, Sendai, Miyagi, 980-8578, Japan(Dated: November 2, 2021)
26

A new search for the diffuse supernova neutrino background (DSNB) flux has been conducted at
Super-Kamiokande (SK), with a 22.5 × 2970-kton·day exposure from its fourth operational phase
IV. The new analysis improves on the existing background reduction techniques and systematic
uncertainties and takes advantage of an improved neutron tagging algorithm to lower the energy
threshold compared to the previous phases of SK. This allows for setting the world’s most stringent
upper limit on the extraterrestrial ν̄e flux, for neutrino energies below 31.3 MeV. The SK-IV results
are combined with the ones from the first three phases of SK to perform a joint analysis using
22.5 × 5823 kton·days of data. This analysis has the world’s best sensitivity to the DSNB ν̄e flux,
comparable to the predictions from various models. For neutrino energies larger than 17.3 MeV,
the new combined 90% C.L. upper limits on the DSNB ν̄e flux lie around 2.7 cm−2 ·sec−1 , strongly
disfavoring the most optimistic predictions. Finally, potentialities of the gadolinium phase of SK
and the future Hyper-Kamiokande experiment are discussed.

I.
A.

∗ Corresponding author: elhedri@llr.in2p3.fr† Corresponding author: giampaolo@llr.in2p3.fr‡ Corresponding author: assy@scphys.kyoto-u.ac.jp§ alsoat BMCC/CUNY, Science Department, New York, New
York, 1007, USA.
¶ Deceased.
k also at University of Victoria, Department of Physics and Astronomy, PO Box 1700 STN CSC, Victoria, BC V8W 2Y2, Canada.

INTRODUCTION

Diffuse supernova neutrino backgroundCore-collapse supernovae (CCSNe) are among themost cataclysmic phenomena in the Universe and are essential elements of dynamics of the cosmos. Their underlying mechanism is however still poorly understood,
as characterizing it would require intricate knowledge ofthe core of the collapsing star. Information about thiscore could be accessed by detecting neutrinos emitted bysupernova bursts, whose luminosity and energy spectra
closely track the different steps of the CCSN mechanism
in a neutrino-heating scenario (see Refs. [1–3], for exam-

3ple). Existing neutrino experiments, however, are mostly
sensitive to supernova bursts occurring in our galaxy andits immediate surroundings, that are extremely rare (afew times per century in our galaxy [4, 5]). Learning
about the aggregate properties of supernovae in the Universe hence requires observing the accumulation of neutrinos from all distant supernovae. The integrated fluxof these neutrinos forms the diffuse supernova neutrino
background, or DSNB (DSNB neutrinos are also referred
to as supernova relic neutrinos (SRNs) in various publications).

The DSNB is composed of neutrinos of all flavors
whose energies have been redshifted when propagating toEarth. Its spectrum therefore contains unique information not only on the supernova neutrino emission process
but also on the star formation and Universe expansion
history. Spectra predicted from various models areshown in Fig. 1. Their overall normalization is mostly
determined by the supernova rate, related to the cosmic
star formation rate. While this redshift-dependent rate
also impacts the DSNB spectral shape, the latter ismainly affected by the effective energies of supernova
neutrinos. Other factors shaping the DSNB spectrum
include the neutrino mass ordering, the initial massfunction of progenitors, the equation-of-state of neutronstars, the fraction of black-hole-forming supernovae, therevival time of the shock wave, etc. Possible combinations of these factors were systematically studied in the
Nakazato+15 model [6], the maximal and minimal fluxesof which are shown in the figure. The minimal prediction
from this model and the Malaney97 gas infall model [7]give the lower bounds on the DSNB flux. Conversely
the Totani+95 model [8] can be considered an upperbound on the expected DSNB flux, on par with the most
optimistic predictions of the Kaplinghat+00 model [9],also shown in the figure. Between these bounds, theDSNB flux can vary over one order of magnitude andits spectral shape bears the imprint of a wide range ofphysical effects. Black-hole-forming supernovae notably
make the DSNB spectrum harder, as can be seen with
the Lunardini09 model [10], that assumes that 17% of
core-collapse supernovae lead to black-hole formation.
Similarly, the Horiuchi+18 model [11] incorporates afraction of black-hole-forming supernovae determinedusing the stars’ compactness. Aside from accounting
for black-hole formation, the Kresse+21 model [12],on the other hand, uses state-of-the-art supernova
simulations to model contributions from helium stars
to the DSNB. Another recent study (Horiuchi+21 [13])discusses the impact of interactions in binary systems, such as mergers, and mass transfer, on the
DSNB flux. Detecting the DSNB would hence provide
valuable insights into a wide array of physical processes.

B.

Experimental searches

While a significant fraction of DSNB neutrinos are expected to have energies lower than 10 MeV —as shownin Fig. 1—, the O(1) MeV region is hard to probe bycurrent experiments due to overwhelming backgrounds
from reactor antineutrinos, solar neutrinos, and radioactivity. Experimentally the DSNB signal has therefore
been searched for at energies of O(10) MeV. At theseenergies, the dominant detection channel in most experiments is the inverse beta decay (IBD) of electron antineutrinos (ν̄e + p → e+ + n), and contributions fromsubleading channels such as electron neutrino elastic scattering or charged current interactions on oxygen can beneglected. This process produces an easily identifiable
positron and a neutron (Fig. 2). Here the neutron does
not emit Cherenkov or scintillation light directly but its
capture on a proton leads to emission of a 2.2 MeV photon. In pure water, the characteristic timescale for neutron moderation and capture is of 204.8 ± 0.4 µs [20].
Pairing the “prompt” positron signal with the “delayed”
photon emission from neutron capture is hence a key
component of many analyses.
Up to now, no evidence for the DSNB signal has been
confirmed and upper limits have been set by various
underground experiments, such as Super-Kamiokande(SK) [21–23], KamLAND [24, 25], SNO [26, 27] andBorexino [28]. Note that, unlike other searches, the SNO
analysis is sensitive to electron neutrinos and, due to theirreducible solar neutrino background, its effective energy
threshold lies around the hep solar neutrino flux endpoint, around 19 MeV. Among all past analyses, the SK
and KamLAND experiments placed the most stringent
upper limits on the DSNB ν̄e flux for neutrino energies
above about 9 MeV while Borexino set the tightest constraints at lower energies. At SK the first DSNB search
was carried out in 2003 using a 22.5 × 1496-kton·day
data set [21]. Using spectral shape fitting for signal
and atmospheric neutrino backgrounds, it placed an upper limit on the DSNB flux for a wide variety of models in the 19.3−83.3 MeV neutrino energy range. This
analysis already allowed to disfavor the most optimistic
DSNB predictions, in particular the Totani+95 model [8],and constrain the parameter space of the Kaplinghat+00
model [9]. In 2012, an improved analysis was performedat SK, using a 22.5×2853-kton·day exposure, a lower neutrino energy threshold of 17.3 MeV, and new event selection cuts allowing a 50% increase in signal efficiency [22].
Below 17.3 MeV, large backgrounds from cosmic muon
spallation and solar neutrino interactions make it extremely difficult to search for DSNB antineutrinos basedon positron identification alone. A search for extragalactic antineutrinos performed at KamLAND in 2021 [25]probed neutrino energies ranging from 8.3 to 30.8 MeV
by investigating coincident positron and neutron capture
signals. While neutron identification in a liquid scintillator detector such as KamLAND is significantly easierthan in pure water, KamLAND’s small fiducial volume

DSNB νe Flux [/cm2/sec/MeV]

4

10
Horiuchi+21 (Extrapolated, α λ =0.1, NO)
Tabrizi+21 (NS+BH, NO)
Kresse+21 (W20-BH2.7- α 2.0, NO)
Horiuchi+18 (ξ12.5,crit= 0.1)

Nakazato+15 (Max, IO)
Nakazato+15 (Min, NO)
Galais+10 (NO)
Horiuchi+09 (6 MeV, Max)
Lunardini09
Ando+03 (updated at NNN05)

10−1

Kaplinghat+00
Malaney97
Hartmann+97
Totani+95

10−2

10−3

10−4

10

20

30

40

50

60
70
80
νe Energy [MeV]

FIG. 1. DSNB ν̄e flux predictions from various theoretical models (Horiuchi+21 [13], Tabrizi+21 [14], Kresse+21 [12],Horiuchi+18 [11], Nakazato+15 [6], Galais+10 [15], Horiuchi+09 [16], Lunardini09 [10], Ando+09 [17], Kaplinghat+00 [9],Malaney97 [7], Hartmann+97 [18], and Totani+95 [8]). Refer to each publication for the detailed descriptions of model. In
the legend, “NO” and “IO” represent neutrino normal and inverted mass orderings assumed in the calculation, respectively.
For the Horiuchi+09 model with a 6 MeV temperature, only the maximal flux prediction is shown. The prediction for the
Galais+10 model here is extrapolated up to 50 MeV as the original publication was served up to 40 MeV. The prediction by
Nakazato+15 is only available up to 50 MeV. The values of the flux used in this analysis for the Ando+03 model are the ones
released at the NNN05 conference [19]. The corresponding flux is larger by a factor of 2.56 than in the original publication [17].

only allowed for an exposure of 6.72 kton·year. In 2015,a new SK analysis including neutron identification led
to stronger constraints down to 13.3 MeV neutrino energies [23]. Since the SK triggers did not allow to record
the neutron capture signal until 2008, this analysis could
only be performed on a small part of the SK data set,with a total livetime of 22.5 × 960 kton·days. Due tothis low exposure and the low neutron tagging efficiency
in water, this search however yielded weaker limits thanthe SK-I,II,III analysis [22] above 17.3 MeV.
In this study, we draw on the previous SK analyses

to present two DSNB searches for antineutrino energies
ranging from 9.3 to 81.3 MeV, with significantly improved background modeling and reduction techniques.
In the 9.3 to 31.3 MeV range, we derive differential upper limits on the ν̄e flux independently from the DSNB
model, following the strategy outlined in Ref. [23] using
a 22.5 × 2970-kton·day data set. In the 17.3 to 81.3 MeV
range we constrain a wide variety of DSNB models using
spectral fits analogous to the ones described in Ref. [22].
We then combine the results of this analysis with theones obtained in Ref. [22] for the former SK phases, thus

5
OD includes 1,885 8-inch PMTs, facing outwards, to detect the Cherenkov light from muons. Further detailed
descriptions of the SK detector and its calibration can befound in Refs. [29–33].

FIG. 2. Schematic illustration of an IBD process and the subsequent neutron capture on another proton. The characteristic neutron capture time in water is τ = 204.8 ± 0.4 µs [20].

analyzing 22.5 × 5823-kton·days of data. This unprecedented exposure will allow to probe the DSNB with anunmatched sensitivity.
The rest of this article proceeds as follows. First, theSK detector and the specific features of its data acquisition system are described in Section II. Then, details
about modeling of the DSNB signal and the different
backgrounds are given in Sections III and IV, respectively. We then describe the data reduction process in
Section V. The procedures associated with the DSNB
model-independent and spectral analyses are described in
Sections VI and VII, respectively. Finally we discuss thecurrent constraints on the DSNB flux and future opportunities at SK with gadolinium and Hyper-Kamiokande
in Section VIII before concluding in Section IX.

II.

SUPER-KAMIOKANDE

Super-Kamiokande is a 50-kton water Cherenkov detector located in the Kamioka mine, Japan. It is structured by a cylindrical stainless steel tank with a diameter of 39.3 m and height of 41.4 m and consists of twoparts: an outer detector (OD) that serves as a muon
veto and an inner detector (ID) where neutrino detection takes place. In order to reduce backgrounds due toradioactivity near the detector wall, most analyses consider only events reconstructed at least 2 m away fromthe ID wall, thus defining a 22.5-kton fiducial volume(FV). It is located 1000 m underground, that allows reduction of the cosmic ray muon flux by a factor of 105 .
In order to ensure a high-quality data taking, conditions
inside SK are tightly controlled; water is constantly recirculated and purified, and the ID wall is covered with
11,129 20-inch photomultipliers (PMTs) with a 3-ns time
resolution, corresponding to a 40% photocathode coverage. These features allow SK to detect particles withenergies ranging from a few MeV to a few TeVs. The

SK is currently undergoing its sixth data taking phasesince it started functioning in 1996. The first phase lasted
1497 days and ended for a scheduled maintenance. Dueto an accident following the maintenance, which resulted
in a loss of 60% of the ID PMTs, SK operated with areduced photocathode coverage for 794 days (phase II).
The coverage was brought back to its nominal value forphase III, that lasted 562 days. SK’s previous spectral
analysis of the DSNB [22] used data from all these threephases. Since 2008, the front-end electronics has been
replaced [34] and a new trigger system that allows forneutron tagging has been set up [35]. SK operated withthis new electronics during phase IV, for 2970 days, until
being stopped for refurbishment work in May 2018. Afterthe maintenance, SK operated for a short time with purewater for calibration and monitoring purposes (phase V)
since January 2019, before being loaded with gadolinium(phase VI) in July 2020. This study will primarily focus
on phase IV, and will present a combined analysis of the
SK-I to IV data. Additionally, we will briefly discuss theopportunities offered by the gadolinium lodaded SK and
by the future Hyper-Kamiokande experiments.
Data processing in SK makes use of multiple triggers,
corresponding to different thresholds on the number of
PMT hits found in a 200-ns window. In SK-I to III, all
PMT hits in a 1.3-µs window around the main activity
peak were stored, using hardware triggers. Since SKIV, the new data acquisition system acquires every PMT
signal and a software trigger system defines events. To
identify positrons from IBD processes, we use the superhigh-energy (“SHE”) trigger, which requires 58 PMT hits
(70 hits before September 2011) in 200 ns. The data ina [−5, +35]-µs window around the main activity peak iscollected with this trigger. This window is too short tocontain the delayed neutron capture signal after IBDs.
In order to identify this signal, a 500-µs (350 µs before November 2010) after trigger (“AFT”) window follows each SHE trigger that is not associated with a predetermined OD trigger. The trigger conditions for thedifferent periods in SK-IV are summarized in Table I. In
the rest of this paper, we will describe the analysis of hit
patterns in these large SHE+AFT windows to identify
the coincident production of a positron and a neutron.
The prompt positron event is reconstructed using thededicated solar neutrino [36–39] and muon decay electron vertex and direction fitter, which is then used toreconstruct the event energy [40]. For this SK phase, we
follow the convention introduced in Ref. [39] and subtract
the 0.511 MeV electron mass from this energy to obtain
the electron equivalent kinetic energy Erec . This quantity can also be interpreted as the total reconstructed
positron energy, and will be used to present the resultsof this study.

6
TABLE I. Summary on the SHE and AFT trigger conditions
and livetime for the different three periods in SK-IV.
SHE-threshold AFT-window
Livetime
70 hits
350 µs25.0 days
70 hits
500 µs
869.8 days
58 hits
500 µs2075.3 days

III.
A.

SIGNAL MODELINGDSNB spectra and kinematics

In this analysis we consider both the DSNB models
whose fluxes are shown in Fig. 1 and models where wevary specific physical parameters. For these parameterized models we compute the DSNB flux using the following formula:Φ(Eν ) =

c
H0

Z X
s

RSN (z, s)X

Fi (Eν (1 + z), s)

νi ,ν̄idz×p,ΩM (1 + z)3 + ΩΛ(1)where Eν is the neutrino energy, c is the speed of lightin vacuum, H0 is the Hubble constant, z is the redshift, RSN (z) is the redshift-dependent supernova rate,and Fi is the supernova neutrino emission spectrum fora given flavor i. The index s represents the different possible classes of supernovae —e.g. supernovae collapsing
into neutron stars or black-holes— associated with specific neutrino emission spectra. The last factor accounts
for the Universe expansion, with ΩM and ΩΛ being thematter and dark energy contributions to the energy density of the Universe, respectively. We evaluate the supernova rate using a phenomenological model describedin Ref. [16], which extracts the redshift dependence ofthe cosmic star formation history by fitting observations
and uses the Salpeter initial mass function [41] to obtainthe fraction of supernova progenitors. We then consider
blackbody neutrino emission spectra, where the effective
temperature Tν takes into account neutrino oscillation
effects. For a given neutrino flavor i the spectrum for amodel with effective temperature Tν is thus given by:

Fi (Eν ) = fν,i Eνtot

2 
−1
120 Eν,i
Eν,i /Tν,i
e
+
1
,
4
7π 4 Tν,i(2)where Eνtot is the total energy of the neutrinos emitted
by the supernovae and fν,i represents the fraction of neutrinos or antineutrinos with flavor i and can be roughly
approximated by 1/6 for a given species. The blackbody
model therefore only depends on the neutrino temperatures and luminosities. In what follows, we will refer toblackbody models with a neutrino effective temperature

T as Horiuchi+09 T MeV models. Note that, due to neutrino scattering in the densest regions of the collapsingstar, supernova neutrino emission is better modeled using a more sharply peaked, “pinched”, Fermi-Dirac spectrum [42]; however, the current exposure at SK does notallow to probe this pinching effect. In this analysis weneglect the effects of degeneracies between the pinchingparameter and e.g. the neutrino temperature on the final
constraints on neutrino emission spectra.
This analysis exclusively focuses on the detection ofelectron antineutrinos via IBD processes. In the rest of
this paper, we model these processes using the StrumiaVissani IBD calculation [43], that approximates the IBD
cross section more precisely than the Beacom-Vogel calculation used for previous SK analyses [22].

B.

Detector simulation

SK-IV has been the longest data taking phase of theexperiment, lasting about 10 years. During this term,
the PMT gain has steadily increased by about 15% over
the whole SK-IV period. Additionally, since the energy
scale change due to the gain change for this analysis was
about 3 %, that effect was corrected in order to maintainthe stability of the detector. A key ingredient of our
DSNB study is therefore a simulation of antineutrino IBD
processes that accounts for the evolution of the detector’s
properties throughout the entire SK-IV period.
In this study we simulate IBD processes uniformly distributed in the entire inner detector, in order to accountfor events close to the ID wall being misreconstructed
inside the fiducial volume. For each vertex we generate a set of three momentum vectors corresponding toan antineutrino, a positron, and a neutron. In a viewof the wide diversity of DSNB models, we generate uniformly distributed positron energies in 1−90 MeV andrenormalize events later to model specific spectra. This
procedure will also allow to use this simulation to model
backgrounds associated to other IBD processes or β decays, as will be mentioned later. We then simulate detection signals using a dedicated simulation, based on
GEANT3 [44], that reproduces the properties of the water in SK, as well as models the PMT properties andelectronic response. Detector properties such as watertransparency and PMT noise are being measured daily,
allowing to accurately model the time-dependent noisecontamination of the positron signal. Neutron capture,
however, produces a 2.2 MeV γ ray whose light yield
is below the SK trigger threshold and particularly difficult to distinguish from PMT dark noise and low energy
radioactive decays. In order to develop robust neutron
identification techniques we therefore collect noise samples from data using a random wide trigger, and inject
them into simulation results, after the positron activity
peak.

7
BACKGROUND SOURCES
A.

Atmospheric neutrinos

Atmospheric neutrinos are an important backgroundin the present DSNB searches. Below about 15 MeV,
neutral-current quasielastic (NCQE) interactions, which
induce nuclear γ ray emission [45–48], form the main
atmospheric neutrino background. On the other hand,charged-current quasielastic (CCQE) interactions andpion production dominate at higher energies. A typicalvisible signature for these interactions is the Cherenkov
light from electrons produced by muon and pion decays.
Hence, the associated reconstructed energies will follow
a Michel spectrum for reconstructed energies of about 15
to 50 MeV. Atmospheric neutrino events are simulatedusing the HKKM 2011 flux [49–51] as an input to theneutrino event generator NEUT 5.3.6 [52] to model interactions. Details about the model parameters used inthis analysis can be found in Ref. [53]. The nuclear deexcitation γ’s are simulated based on the spectroscopic
factors for the p1/2 , p3/2 , and s1/2 states calculated inRef. [45]. More detailed descriptions of the excited states
can be found in Refs. [46, 48] for the T2K NCQE measurement. One difference between our procedure andRef. [48] is treatment of the others state —a state affected by short-range correlations or with a very highexcitation energy—. This state is treated as the highest
excitation state (s1/2 ) in Ref. [48] while it is included in
the ground state (p1/2 ) in the present analysis as it isfrom the latest atmospheric neutrino analysis in SK [54].
The systematic uncertainty regarding this treatment is
considered in the scaling factor used in Section VI A.
The detector simulation for atmospheric neutrino
events is performed using the same GEANT3-based simulation tool as for the signal, but without corrections for
the PMT gain shift on the primary event. The associated systematic error is estimated in Section VI A andfound to be subdominant. Finally, random trigger data
are overlaid on top of neutron capture signals in order tosimulate background for neutron tagging, following theprocedure described in Section III B. Note that for thislast step the time evolution of the detector properties istaken into account.

B.

Cosmic ray muon spallationof this analysis. Since most spallation isotope decays
only produce β and γ rays, spallation backgrounds canbe significantly reduced using neutron tagging. Nonetheless, accidental pairing between prompt β or γ events and
PMT hits due to dark noise or intrinsic radioactivity willallow a sizable number of spallation events to pass as
IBDs. Furthermore, a few isotopes undergo a β+n decay that mimics the IBD signal, and thus cannot be removed using neutron tagging. Efficient spallation reduction therefore requires both dedicated spallation reduction techniques tailored to the properties of the different
isotopes, and neutron tagging cuts.
Relevant spallation isotopes and their visibility in water have been studied with the FLUKA simulation [55] inRefs. [56–58]. The isotope lifetimes and the end-point energies of their decays are summarized in Fig. 3. The highest end-point is 20.6 MeV, from 14 B and 11 Li. The spallation cuts used in this analysis will therefore be applied
up to 23.5 MeV reconstructed kinetic electron energy inorder to account for energy resolution effects. Note alsothat the isotope half-lives range from O(0.01) sec to a few
tens of seconds. Given the high rate of cosmic ray muons,pairing a given isotope decay with its parent process is aparticularly difficult endeavour.

11Be

101

16N
15C

Half-life [sec]

IV.

100

8B 8Li

16C

9Li

8He

10 1

12Be

10 2

9C

12B
13B

14B

12N
13O

8

10

12

14

16

18

11Li

20

22

End-point energy [MeV]In spite of its 2700-m water-equivalent overburden, SK
is still exposed to cosmic ray muons with a rate of ∼2 Hz.
Interactions of these muons in water produce electromagnetic and hadronic showers. Spallation of oxygen nuclei
induced by the muons or by secondary particles can thenlead to the production of radioactive isotopes, whose decays can be misidentified as IBD events in the DSNB
search window. Below about 20 MeV, the associatedbackground is 106 times higher than the DSNB flux predictions, making spallation reduction an essential aspect

FIG. 3. End-point energies and half-lives of the β decay isotopes produced by cosmic ray muon spallation in water, whose
decay products are found in the DSNB search window. The
isotopes represented by a red circle (a black square) decay
with (without) neutrons.

Isotopes that undergo β+n decays cannot be rejectedusing neutron tagging and therefore need to be modeled with particular care. This is notably the case for88He, 11 Li, 16 C, and 9 Li. Here 11 Li is particularly shortlived and therefore easy to eliminate using mild spallation cuts. Its production yield is also expected to beparticularly low (around 10−9 µ−1 ·g−1 ·cm2 [56]). In addition, 8 He and 16 C will remain subdominant due totheir low production yields (around 0.23 and 0.02 ×
10−7 µ−1 ·g−1 ·cm2 [56]) and endpoint β decay energies.
On the other hand, 9 Li has a non-negligible yield (around
1.9×10−7 µ−1 ·g−1 ·cm2 [56]), a medium half-life of about
0.18 sec, and decays into a β+n pair with a branchingratio of 50.8%. The associated β spectrum is shown in
Fig. 4, and extends up to 13.5 MeV reconstructed energy
when accounting for resolution effects. The β+n decayof 9 Li is hence similar to the IBD of a DSNB neutrino,
except that the β from 9 Li is an electron and the neutron energy is higher than that from IBDs. Since SK
does not distinguish electrons from positrons and does
not allow to measure the neutron energy, we model 9 Li
decay using the IBD Monte-Carlo simulation, renormalizing events to reproduce the 9 Li β+n spectrum. Finally,
we estimate the total expected number of 9 Li β+n decays in the DSNB search window by using an earlier SK
measurement [59] that estimated the 9 Li production rate
to be 0.86 ± 0.12(stat.) ± 0.15(syst.) kton−1 ·day−1 . Notethat the difference in neutron energies with respect to
IBD processes leads to a systematic uncertainty. However, this uncertainty is taken into account when calibrating neutron tagging efficiencies, as neutrons from thecalibration source have similar energies to the 9 Li decay
products. More details about this calibration procedureare given in Section V D.

the DSNB signal at low energies. We estimate the reactor antineutrino flux using the SKReact code [60] developed for the SK reactor neutrino analysis. This codeuses the reactor neutrino model from Ref. [61] based on
IAEA records [62] to evaluate the electron antineutrino
fluxes for any reactor in the world, and simulates oscillation effects. The resulting flux prediction also accounts
for the time dependence of the reactor antineutrino flux
due to changes in reactor activity. The expected total
reactor neutrino flux at SK can thus be predicted up toEνmax = 9.2 MeV. The expected reactor neutrino spectrum is obtained from this predicted flux using the IBD
simulation, and is shown in Fig. 5. Above 6 MeV, as
can be seen in this figure, the reconstructed spectrum
is primarily shaped by resolution effects. Since the reactor antineutrino spectrum has been constrained up to
12 MeV by Daya Bay [63], possible antineutrino contributions above Eνmax can be safely neglected in this study.
In the 7.5−9.5 MeV reconstructed kinetic electron energy range, the reactor antineutrino background is about
6 times higher than the DSNB signal predicted by the
Horiuchi+09 6 MeV model [16]. Reactor antineutrino
backgrounds thus effectively set a lower energy thresholdfor the DSNB search.

FIG. 5. True positron kinetic energy (blue) and reconstructedkinetic energy (orange) spectra from reactor neutrino IBD interactions. The y-axis shown the average number of events
per year at SK-IV. The event rate is not constant with time
due to reactors being powered on during the data taking period.

FIG. 4. True electron kinetic energy (blue) and reconstructedkinetic energy (orange) spectra from the 9 Li β+n decay.

C.

Reactor neutrinos

Antineutrinos emitted from nearby nuclear reactors
can undergo IBDs in SK and therefore be mistaken as

D.

Solar neutrinos

Electron neutrinos from the Sun cause elastic interactions with electrons in SK. In the DSNB search window,
the 8 B and hep solar neutrino fluxes represent an important background up to ∼20 MeV. This backgroundcan be drastically reduced by neutron tagging, althougha small fraction of electrons from solar neutrinos will be9accidentally paired with background fluctuations or radioactive decays, like for the spallation backgrounds. In
the absence of neutron tagging, a solar neutrino interaction can still be identified by comparing the reconstructed
directions of the scattered electron to the direction of the
Sun at the detection time. The angle between these two
directions, called the opening angle, is expected to beclose to zero for solar electron neutrinos scattering elastically. It is however smeared by reconstruction effects
that need to be modeled accurately. In this study, we use
the simulation designed for the latest SK-IV solar neutrino analysis [39] to model solar neutrino spectra andevaluate the impact of opening angle cuts.
Finally, note that in this analysis we do not study the
impact of solar antineutrino production, since the rate
predicted by the Standard Model (SM) is negligible. Antineutrino production via beyond the SM processes would
yield a signal that can only be distinguished from the
DSNB signature by spectral shape studies. So far, only
upper limits on the solar antineutrino flux have been derived at SK, notably in Ref. [53].

V.

DATA REDUCTION

We apply reduction cuts to the data taken by the
SHE trigger in SK-IV, corresponding to an exposure of
22.5 × 2970.1 kton·days. The SHE trigger efficiency isclose to 100% over the entire analysis window. The following sections describe the four reduction steps used inthis analysis: noise reduction in Section V A, spallation
reduction in Section V B, DSNB positron candidate selection in Section V C, and neutron tagging in Section V D.

A.

Search energy range and noise reduction

The lower energy threshold for this analysis is the SHE
trigger threshold, which corresponds to Erec > 7.5 or
9.5 MeV depending on the observation time (see Table I). The upper bound of the DSNB analysis window is
set to 29.5 MeV for the model-independent analysis and79.5 MeV for the spectral analysis.
We first select SHE-triggered events with Erec below
79.5 MeV and apply a set of cuts aimed at removing
PMT noise, calibration events, radioactivity from thesurrounding rock and the detector wall, and cosmic ray
activity. In particular, we require candidate events to belocated at least 2 m away from the ID wall, thus defining
a fiducial volume (FV) of 22.5 kton. Additionally, weremove events associated with an OD trigger, or within
50 µs of another event with reconstructed electron kinetic
energy larger than 6 MeV, in order to remove electrons
produced by the decay of low energy cosmic ray muons
stopping in the detector. In addition, we apply fit quality
cuts to eliminate noise-like events. Inside the FV thesenoise reduction cuts have a signal efficiency larger than
99% as confirmed by the IBD Monte-Carlo simulation.

Finally, we require all SHE events to be followed by an
AFT trigger to allow for neutron tagging. Due to occasional deadtime induced by trigger software failure, thisstep is associated with a ∼94% signal efficiency. Events
passing these criteria will be subsequently referred to as“DSNB candidates”.

B.

Spallation reduction

Spallation backgrounds largely dominate over putative
DSNB signals for reconstructed electron kinetic energies
lower than 23.5 MeV. Here, we present a set of dedicatedcuts that will allow for a substantial reduction of thesebackgrounds, even in the absence of neutron tagging.

1.

Spallation preselection

We apply two sets of preselection cuts to reduce contributions from the most energetic muons and some longlived spallation isotopes such as 16 N, with minimum
harm to the signal efficiency.
First, since energetic muons could induce the production of multiple radioactive isotopes, we remove all DSNB
candidates observed within 60 sec and 4.9 m from at least
one other low energy event with Erec in the 5.5-24.5 MeV
range [64]. The 60 sec time window has been chosen to
include decays of abundant and long-lived isotopes suchas 16 N while the 4.9-m distance cut provides optimal signal over background separation. We estimate the signal
efficiency of this so-called “multiple spallation cut” usinga sample of low energy events from radioactive decays
at SK (4 < Erec < 5 MeV), whose reconstructed vertices
have been replaced by random vertices inside the ID. The
cut parameters introduced above allow to remove about
45% of the background with a 98% signal efficiency, as
can be seen in Fig. 6. Since this cut removes low energy
events clustered in space and time, it also removes low
energy muons that are misclassified as electrons with tensof MeV energies. These muons are not targeted by theother spallation cuts applied in this study, which assume
that the prompt event is an electron or a positron. Consequently, in this study, we will apply this cut over thewhole energy range for each of our analyses, including
the spallation-free Erec > 23.5 MeV region.
In addition to removing multiple spallation events, welocate muon-induced showers by identifying neutron captures observed less than 500 µs after muons. These neutron clusters, called “neutron clouds”, are often produced
in muon-induced hadronic showers, that are the birthplaces of spallation isotopes. Neutron capture events following muons are collected separately from DSNB candidates, using a specific trigger called the Wideband Intelligent Trigger (WIT) [65–67]. This trigger, in addition
to setting a threshold on the number of PMT hits, applies cuts on the event quality and distance of the reconstructed vertex to the ID wall. Events triggered by WIT

10

FIG. 6. Minimal distance between a given DSNB candidate(blue) or SK low energy event with a random vertex (orange)
and a well-reconstructed event with Erec > 5.5 MeV observed
within 60 sec of it. The DSNB candidates are defined as
all the SK-IV events with Erec > 7.5 MeV passing the noise
cuts. Here we also show the 4.9 m separation required in thisanalysis with a red vertical line.

less than 500 µs after a given muon are considered as
neutron captures if their reconstructed vertex is within
5 m of the reconstructed muon track. If multiple neutron
captures are observed after a muon event, we define aneutron cloud using the criteria from Ref. [64]. We thenremove DSNB candidates that are found close in timeand space to these neutron clouds, following the procedure introduced in Ref. [64]. More details about the cut
criteria are given in Appendix A. To estimate the signal
and background efficiencies of the neutron cloud cut, wepair muons with DSNB candidates in data, defining two
samples: a “pre-sample”, including muons found up to
60 sec before a DSNB candidate, and a “post-sample”,
with muons found up to 30 sec after a candidate. While
the pre-sample contains a mixture of correlated and uncorrelated pairs, the post-sample contains only uncorrelated pairs. This sample can hence be used to both
readily estimate signal efficiencies and characterize theproperties of the pairs formed by each isotope decay andtheir parent muon in the pre-sample. An example of
how to extract the distribution of the distances between
spallation events and the neutron clouds of their parent
muons is shown in Fig. 7. We optimize the neutron cloud
cut by estimating the shape of the neutron clouds fromdata and maximizing the statistical significance of thesignal over the spallation background. The performanceof this cut is however currently limited by the weaknessof the neutron capture signal and the fact that the WIT
trigger has only been available only during the last 388
live days of SK-IV. While the signal efficiency is larger
than 99% this cut only removes about 10% of the spallation background in this analysis. Over the 388 days of
the WIT trigger livetime, however, it removes about 40%of spallation isotope decays.
Finally, it should be noted that the multiple spallation

FIG. 7. Distributions of distance between DSNB candidates
and neutron clouds for pre- and post-samples, where the selected neutron clouds have been observed up to 60 sec before
and after the DSNB candidates respectively. The distribution for the post-sample has been rescaled by assuming that
contributions from spallation pairs beyond 20 m are negligible. Contributions from pairs formed by spallation isotope
decays and neutron clouds associated with their parent muonsappear as a peak in the pre-sample below 10 m. The associated distribution can be extracted by subtracting the rescaled
post-sample from the pre-sample.

and the neutron cloud cuts described above are expected
to overlap since muons associated with large hadronic
showers are also likely to lead to the production of multiple isotopes. Accounting for this overlap, these preselection cuts allow to remove about 55% of the spallation
background when neutron cloud data are available.

2.

Searching for parent muons

Since preselection cuts remove only about half of thespallation background, a more in-depth study of the correlations between muons and DSNB candidates is necessary. In particular, in order to identify the decays of spallation isotopes, it is crucial to associate each isotope with
its parent muon. Due to the long half-lives of isotopes
such as 11 Be and 16 N, we need to investigate all possible
pairings between DSNB candidates and muons observed
up to 30 sec before them. This window size allows toaccommodate most decay products, without processing
a prohibitively high number of muons. As with neutron
cloud cut optimization, we define pre- and post-samples
in order to extract the observable distributions associatedwith spallation pairs and estimate signal efficiencies.
We define candidate muons by selecting events associated with more than 31 PMT hits (34 before May 2015)
in 200 ns and depositing more than 500 p.e. in the ID.
Additionally, since the SK trigger windows can containmultiple events, we also look for muons around the main
activity peak in these windows. For the same reason,

11
calibration trigger windows are also thoroughly investigated. The properties of the muon candidates, such as
the charge deposited in the detector, the entry point,
and the direction of the track, are then extracted using
a dedicated fitter [22, 68, 69]. This fitter classifies muons
into five categories: misfits (1.0%), single through-going
muons (82.2%), stopping muons (4.9%), multiple-track
muons (7.6%), and corner clippers (4.3%); the fractionsof the different categories in the SK-IV data are shownin parentheses. Misfit muons with a charge lower than
1000 p.e. are removed in order to reject non-muonic highenergy events. This cut does not affect removing spallation background because the track length of these low
charge muons through the ID is typically less than 50 cm.
After having paired DSNB candidates with the muonsselected above, we extract the following observables, related to the intrinsic properties of the muons and theircorrelations with DSNB candidates.

where qMI and L are the number of p.e. per centimeter expected from the minimum ionization andthe track length, respectively. Here qMI depends
on detector properties, such as water transparency,that vary with time and is therefore recomputed foreach run. For muons with multiple tracks, L is thelength of the first track. This observable allows todetermine how likely a given muon is to induce ashower.

muon track

position wheredE/dx is largestll• dt: time difference between a muon and the DSNB
candidate. For spallation events, this observablereflects the half-lives of the produced isotopes.
• `t : transverse distance of the DSNB candidate to
the muon track (see Fig. 8). This observable is related to the path lengths of the secondary particles
from muon showers. For well-fitted single throughgoing muons paired with their associated spallation
isotope decays, `t is typically no larger than a few
meters.
• `l : longitudinal distance to the point of the muontrack associated with the maximal energy deposition (see Fig. 8). The energy deposition along
the muon track is determined by projecting backthe light seen by each PMT to a point on themuon track, following the procedure described inRef. [22]. So `l provides an estimate of the distance between the low energy event and the originof the particle shower and, if the shower position iscorrectly identified, should also not be larger thana few meters for spallation pairs.
• Qµ : total charge deposited by muons in the detector. Higher values of this observable are expected when a shower is produced. While thisobservable also includes contributions from minimum ionization along the muon track —and hence
strongly depends on the muon track length, it provides the most robust shower predictor for poorlyfitted muons and muon bundles. Conversely, forwell-fitted single through-going muons, the shower
can be more accurately characterized by the Qres
observable described below.
• Qres : residual charge deposited by a muon compared to the value expected from the minimum ionization. This observable can be expressed as:
Qres = Qµ − qMI · L,DSNB candidate

lt

FIG. 8. Schematic of the spallation observables `t and `l . The
cylindrical shape represents the SK ID.

An example set of these spallation variables for single
through-going muons is shown in Appendix A. The distributions associated with corner clipping muons show
no evidence for spallation. Moreover, as discussed inRef. [64], spallation background events associated withcorner clipping muons are expected to represent less than
10−5 of the total number of spallation events. In the restof this study, we will therefore neglect contributions fromthese muons to spallation backgrounds.
Showers induced by muon spallation can be extremely
energetic and involve up to thousands of particles, notably neutrons and pions. Such neutrons cause nuclear
reactions and possibly produce deexcitation γ rays, andare later captured at a timescale up to ∼500 µs. The
MeV-scale γ rays or β particles —produced by the initial neutron, secondary nuclear reactions, or the decay
product of a short-lived isotope, followed by a 2.2 MeV
γ ray from neutron capture can mimic an IBD signature. To prevent these from leaking into the sample, we
require DSNB candidates to be more than 1 ms awayfrom any muon. The impact of this cut on the signal efficiency is negligible. Additionally, we apply a set of cutson dt and `t for different energy ranges, exploiting thedependence of the isotopes’ half-lives in their endpoint
energies shown in Fig. 3. Detailed descriptions of theserectangular cuts are given in Appendix A. These cuts
are particularly efficient above 15.5 MeV where shortlived isotopes dominate. In 15.5−19.5 MeV, notably,
they allow to remove about 85% of the spallation back-12
ground while keeping 88% of the signal events. Finally,in order to further eliminate spallation backgrounds, we
use distributions of the different observables consideredhere to define log-likelihood ratios. First, we prepare
two probability density functions (PDFs): the spallation
PDF (PDFispall ) and the random PDF (PDFirandom ), foreach variable i = dt, `t , `l , Qµ , Qres . We isolate contributions from spallation events by subtracting the postsample distributions from the corresponding pre-sample
distributions and obtain PDFispall after area normalization of these subtracted distributions. For PDFirandom wenormalize post-sample distributions by their areas. We
repeat this procedure for each category of muons. For single through-going and multiple track muons, that generate most of the spallation background, we tune the PDFs
in dt and `t bins. These bins, adjusted by considering theisotopes’ typical half-lives and the `t distributions, account for the correlations between these observables. For
each set of PDFs, we then define log-likelihood ratios as:

Lspall = log

Y PDFispalli

PDFirandom

!
.

(3)

Separate likelihoods are defined for each muon category,
except for corner clippers, whose contributions to spallation backgrounds are negligible as mentioned above.
Note that for misfit muons only dt is used to calculate
log-likelihood ratios, as the other observables are not reliable. An example distribution of the log-likelihood ratio
is given as Fig. 33 in Appendix A. We finally determine
cut conditions for each likelihood ratio, accounting forcomplex correlations between spallation observables, geometrical and muon reconstruction effects.
The signal efficiencies and background rejection ratesof the resulting cuts are estimated using the random
sample introduced in this section, as well as spallationdominated data samples. These samples as well as our
methodology to estimate the spallation cut performanceare described in detail in Appendix A. Note that the estimates of the spallation remaining rate presented there areused only for cut optimization and not for the final background predictions shown in Sections VI and VII. The
current cuts achieve significant improvement over theprevious cuts used in Refs. [22, 23], especially at low energies; for the same spallation remaining rate, the signal
efficiency is increased by up to 60% for Erec < 11.5 MeV,
20% for 11.5 < Erec < 13.5 MeV, and comparable forhigher energies.

we define the following discriminating observables, aimed
at characterizing the prompt event.

1.

Radioactivity near the detector wall as well as muon
spallation in the rock surrounding the detector, can leadto electrons or γ rays entering FV. Instead of tightening the FV cut, we consider the effective distance ofeach event to the ID wall deff [22, 36]. This observable
is computed by following the reconstructed direction ofeach event backwards from its reconstructed vertex to the
ID wall, and can be interpreted as the minimal distanceneeded for a radioactive particle produced near the wall
to travel to the event vertex, as shown in Fig. 9. The deff
distributions for the DSNB signal and for the events selected using the cuts outlined in Section V A are shown in
Fig. 10, for the reconstructed energy ranges corresponding to the model-independent analysis, 7.5 to 29.5 MeV,
and the spectral analysis, 15.5 to 79.5 MeV. Comparing
the deff distributions from the data and the Monte-Carlosimulation allows to estimate contributions from radioactivity near the wall, and determine a suitable deff cut. In
this study we impose lower thresholds on deff ranging
from 3 to 5 m depending on the reconstructed energy:Erec − 15.5 MeV
deff > max 300, 500 −
× 50
1 MeV

DSNB positron candidate selection

In the region above ∼20 MeV, spallation backgrounds
can be reduced to negligible levels using a series of spallation cuts while keeping most of the signal. However,significant backgrounds from atmospheric neutrino interactions and radioactive decays remain. To identify them,

cm.

FIG. 9. Views of the ID, showing the effective distance ofa reconstructed event to the detector wall. The deff is the
distance that a particle emitted near the ID wall —e. g. fromby a radioactive decay— needs to travel to give the observed
signal.

2.
C.

Incoming event cut

Pre- and post-activity cuts

Atmospheric neutrino interactions can produce bothprompt signals from photons, electrons, or energetic
heavy particles, and delayed signals from decay of muonsand pions to electrons. These signals will be typicallyseparated by a few microseconds and can therefore sharethe same trigger window. Depending on which particle deposits most light, unusually high activity can beEvents/100 cm

Events/100 cm

13

50000

40000

350
300
250
200

30000

150
20000
100
10000

0
effwall [cm]

0.8
0.6

500

1000

1500

2000

2500

3000

3500

4000

500

1000

1500

2000

2500

3000

3500
4000
deff [cm]

0.4
0.2
0

0
0.8
0.6
0.4
0.2
0
−0.2
−0.4

|data - MC|/MC|data - MC|/MC

50

effwall [cm]
500

1000

1500

2000

2500

3000

3500

4000

500

1000

1500

2000

2500

3000

3500
4000
deff [cm]FIG. 10. The deff distributions for the observed data and IBD Monte-Carlo simulation, for events with Erec ∈ [7.5, 29.5] MeV(left) and Erec ∈ [15.5, 79.5] MeV (right). The rate of IBD Monte-Carlo simulation events has been scaled to match the numberof events observed in the data for deff > 1000 cm. The lower panels show the relative difference between the numbers ofdata and IBD simulation events. We consider only events reconstructed in the FV and passing the noise reduction cuts. The
assumed DSNB spectrum is taken from the Horiuchi+09 6 MeV model [16].

noticed either before or after the main peak. For premax
activity we compute the maximal number of hits, Npre,
in a 15-ns time-of-flight subtracted window between 5 µs
and 12 ns before the main peak. We apply a similar procedure to post-activity, using an algorithm developed forother analyses in SK and T2K that computes the number of electrons from muon and pion decay Ndecay-e [70].
The Ndecay-e distributions for low and high energy events
are shown in Fig. 11. In the analyses discussed here we
max< 12 and Ndecay-e = 0.
require Npre
3.

Cherenkov angle

The opening angle of the Cherenkov light cone (θC )emitted by highly relativistic particles, electrons andpositrons in the current context, in pure water is around42◦ . Conversely, at the O(10) MeV energies considered
here, heavier particles like muons and pions will be typically observed near the Cherenkov threshold, leading tocones with smaller angles in the current analysis range.
In addition, NCQE interactions with multiple photon
emission can produce multiple overlapping Cherenkov
cones, that will be reconstructed as a single cone with aparticularly large opening angle [46, 48]. Consequently,
θC is one of the most powerful observables for reducing both NCQE and µ/π-producing atmospheric backgrounds. The Cherenkov angle distributions for the signal and atmospheric neutrino events are shown in Fig. 12
for the energy ranges of the model-independent analysis
and the spectral analysis. In what follows we will require the Cherenkov angle in the signal regions to be in[38, 50]◦ . In the spectral analysis, we will also use thisangle to define sidebands to evaluate atmospheric backgrounds.

4.

Ring clearness

Electrons and positrons do not follow a straight trajectory in SK due to scattering and bremsstrahlung, whichleads to fuzzy Cherenkov rings. Pions, on the other hand,lead to well-delineated ring patterns. In order to characterize this property we consider a 15-ns time-of-flight
subtracted window around the main activity peak andcompute the opening angles of the cones formed by the
directions of all possible 3-hit combinations. We thenidentify the peak of this opening angle distribution θ0
and estimate the “clearness” of the ring by computing:Lclear =Ntriplets (θ0 ± 3◦ ).
Ntriplets (θ0 ± 10◦ )(4)The Lclear distributions from the atmospheric and IBD
Monte-Carlo simulations are shown in Fig. 13. Eventswith lower Lclear have fuzzier rings. For the analyses
presented here we require Lclear < 0.36.
5.

Average charge deposit

Energetic muons scatter less than electrons and hence
often deposit more charge in individual PMTs than

600

CCQE-like
CC1π

EventsEvents14

CC other
NCQE
NC1π
NC other

500

400

CCQE-likeCC1π
CC other
NCQE
NC1π
NC other

1400
1200
1000
800

300

600
200
400
100

200

0
0

1

2

0
0

3
Ndecay-e123Ndecay-e

CCQE-likeCC1π
CC other
NCQE
NC1π
NC other
IBD signal

350
300
250

Events/2.7-degree

Events/2.7-degree

FIG. 11. The Ndecay-e distributions for the atmospheric neutrino Monte-Carlo simulation, for events with Erec ∈ [7.5, 29.5] MeV(left) and Erec ∈ [15.5, 79.5] MeV (right). We consider only events reconstructed in the FV and passing the noise reduction
cuts.

1200

CCQE-likeCC1π
CC other
NCQE
NC1π
NC other
IBD signal

1000
800

200
600
150
400
100
200

50
0

20

30

40

50

60

70

0

80
90
θC [degree]

20

30

40

50

60

70

80
90
θC [degree]

CCQE-likeCC1π
CC other
NCQE
NC1π
NC other
IBD signal

100

80

Events/0.005

Events/0.005

FIG. 12. The θC distributions for the atmospheric neutrino and IBD Monte-Carlo simulations, for events with Erec ∈[7.5, 29.5] MeV (left) and Erec ∈ [15.5, 79.5] MeV (right). We consider only events reconstructed in the FV and passing
the noise reduction cuts. The assumed DSNB spectrum is taken from the Horiuchi+09 6 MeV model [16]. The signal andbackground distributions are normalized to the same area.

CCQE-likeCC1π
CC other
NCQE
NC1π
NC other
IBD signal

300
250
200

60
150
40
100
20

0
0.25

50

0.3

0.35

0.4

0.45

0.5

0.55
0.6
ring clearness

0
0.25

0.3

0.35

0.4

0.45

0.5

0.55
0.6
ring clearness

FIG. 13. The Lclear distributions for the atmospheric neutrino and IBD Monte-Carlo simulations, for events with Erec ∈[7.5, 29.5] MeV (left) and Erec ∈ [15.5, 79.5] MeV (right). We consider only events reconstructed in the FV and passing
the noise reduction cuts. The assumed DSNB spectrum is taken from the Horiuchi+09 6 MeV model [16]. The signal andbackground distributions are normalized to the same area.

15
O(10) MeV electrons and positrons. We exploit thisfeature by considering a 50-ns time-of-flight subtracted
window around the main activity peak and calculate theaverage charge deposited per PMT hit (q50 /n50 ) in thiswindow. The distributions are given in Fig. 14. We require q50 /n50 to be less than 2 in the analysis.
6.

Cut efficiencies and systematic uncertainties

We estimate the signal efficiencies of most of the cutsdetailed above using the IBD Monte-Carlo simulation.
Hence, the main source of systematic uncertainties on
the cut efficiencies stems from the modeling of the SK
detector and particle propagation in water. We estimate
these uncertainties using calibration data taken by injecting a monoenergetic electron beam produced by a
LINAC at different positions inside the detector [32]. For
a given cut, we compare the measured and predicted efficiencies using a dedicated Monte-Carlo simulation, anddefine the systematic uncertainty as the maximum discrepancy over all calibration tests. LINAC calibration
thus allows to estimate uncertainties for ring clearness,
the average charge, and Cherenkov angle cuts. Conversely, for the incoming event cut —that requires considering uniformly distributed events— we use a strategy
developed for the SK solar neutrino analyses [36–39]: wemeasure the variation of the deff cut efficiency after shifting the event directions and vertices in the IBD MonteCarlo simulation. The vertex and direction shifts used
for these estimates are obtained from calibration studies
using a nickel source [31]. Finally, efficiencies for the preand post-activity cuts can be directly estimated using aspallation-dominated sample of low energy events, withnegligible systematic uncertainty. Overall, positron candidate selection cuts allow to remove a large fraction of
atmospheric backgrounds while keeping up to 85% of thesignal. The total systematic uncertainty computed for
positron candidate selection cuts using these procedures
are of a few percent.

D.

Neutron tagging

The introduction of a new trigger scheme at SK-IV has
made characterizing IBDs via neutron tagging possible,
by requiring the prompt and delayed signals be detected
within 500 µs of each other. In what follows, we will
devise a neutron tagging algorithm tailored to the DSNB
analysis, inspired by previous SK studies from Refs. [71,
72].

1.

DSNB candidate selection and neutron preselection

As explained in Section II, most SHE triggered eventsare followed by a 500 µs AFT trigger window in orderto record the neutron capture signal. Before attempting

to identify neutron candidates in this combined trigger
window, we apply a cut on the maximum hits inside a
200-ns window: N200 < 50 to eliminate SHE+AFT trigger windows containing muons. This cut is associatedwith a signal efficiency larger than 99%. After applying
the N200 cut we scan the SHE+AFT trigger windows of
the remaining events to select hit clusters that could be
associated with neutron captures. Here, we take advantage of the fact that a typical neutron capture vertex iswithin a few tens of centimeters of the well-reconstructed
IBD vertex. Since the SK resolution is about 50 cm forO(10) MeV events, this proximity allows to reliably estimate the required photon emission time for each PMT
hit to be the result of the neutron capture. A neutron
candidate is therefore defined as a group of hits that isclustered in emission time. We then look for clusters
maximizing N10 , the number of hits in a 10-ns window.
In previous SK-IV searches, neutron candidates were
defined as clusters of hits satisfying N10 ≥ 7 [71]. This
preselection cut has an efficiency of around 35% and was
motivated by the energy range of these searches, where
spallation backgrounds were largely dominating. In thisstudy, we modify this preselection step as follows. First,
we apply a cut to reduce the continuous PMT dark noise,likely due to scintillation in the PMT glass, that causesa single PMT to flash more than once in rapid succession with a timescale of ∼10 µs. To do so, for each PMT
we remove hit pairs separated by less than 6 (12) µs for
N10 > 6 (N10 ≤ 6). After applying this cut, we define neutron candidates as hit clusters verifying N10 ≥ 6.
This preselection cut will allow loosening of the neutron
tagging cut in energy regions where spallation no longer
dominates, in particular for the spectral analysis. Finally,
since the continuous dark noise cut relies on investigating PMT hit pairs, its efficiency sharply increases near
the edge of the SHE+AFT trigger window. We restrict
the neutron search window to [14, 523] µs ([14, 373] µs forevents with a 350-µs AFT trigger window). After thesedifferent steps, neutron candidate preselection is associated with an efficiency of 44.7%.

2.

Boosted Decision Tree

After preselection, we are still left with a large background from e. g. dark noise fluctuations and radon decays [73], totaling about 7 neutron candidates per event.
To further reduce this background, for each neutron candidate, 22 discriminating variables are calculated. Thesevariables are related to specific features of noise events(such as clustered PMT hits), the Cherenkov light pattern of the neutron candidate, and its position with respect to the primary event vertex. These 22 observables
are then used in a Boosted Decision Tree (BDT) implemented using the xgboost Python library [74]. The BDT
is trained using a data set of 2.8 × 108 neutron candidates. In this sample, about 2 × 106 candidates are neutron captures simulated using the procedure described

120

Events/0.5

Events/0.5

16

CCQE-likeCC1π
CC other
NCQE
NC1π
NC other
IBD signal

100
80

CCQE-likeCC1π
CC other
NCQE
NC1π
NC other
IBD signal

350
300
250
200

60

150
40
100
20
0

50
1

1.5

2

2.5

3

3.5

4
Q50/N50

0

1

1.5

2

2.5

3

3.5

4
Q50/N50

FIG. 14. The q50 /n50 distributions for the atmospheric neutrino and IBD Monte-Carlo simulations, for events with Erec ∈[7.5, 29.5] MeV (left) and Erec ∈ [15.5, 79.5] MeV (right). We consider only events reconstructed in the FV and passing the noise
reduction cuts. The assumed DSNB spectrum is taken from the Horiuchi+09 6 MeV model [16]. The signal and background
distributions are normalized to the same area.

in Section III B and the rest is composed of low energy
accidental coincidences from random trigger events. Table II summarizes the BDT training parameters. The
final BDT performance is shown in Fig. 15. Even accounting for preselection cuts, this neutron tagging algorithm allows reduction the background rate by a factorof 7 compared to the analysis presented in Ref. [72] forthe same signal efficiency.
TABLE II. Parameters used for the training of the neutron
tagging BDT in the xgboost Python library framework.
Parameter Valuelearning rate 0.02522
subsample
0.97
max depth
10tree method approx
max iterations
1500
best iteration
1170

3.

Systematic uncertainty on neutron tagging efficiency

To estimate the systematic uncertainty on the BDT
efficiency associated with mismodeling of neutron propagation in water and the 2.2 MeV photon emission, we
compare the efficiencies predicted by the IBD MonteCarlo simulation to measurements performed using anamericium-beryllium (241 Am/Be) radioactive source embedded in a BGO scintillator. A detailed description ofthis procedure is given in Ref. [35]. To calculate theneutron identification efficiency for data, we follow themethodology described in Ref. [76] and fit the timing distribution of tagged neutrons by a decaying exponential
plus a constant term as shown in Fig. 16. True neutron
capture event counts decay exponentially in time fromthe detection time of the prompt event, while uncorrelated backgrounds are constant in time. By comparing

FIG. 15. Neutron identification efficiency and mistag rate of
BDT used for the current analysis. The red band depicts theuncertainty associated with the time-dependence of the noise
in SK, as described in Section VI A 4. The fraction of remaining IBD neutrons and the number of accidental coincidences
per event after preselection, averaged over the entire SK-IV
period, are 44.7% and 7, respectively. These values depend
linearly on the PMT noise [75], and hence vary by about 15%
over the SK-IV period.

the efficiency in data and simulation for BDT cut points
relevant in our analysis, as described in Appendix B, weassign a relative uncertainty of 12.5% for all neutron tagging cuts.
E.

Solar neutrino cut

Since even mild neutron tagging cuts can reduce solar neutrino backgrounds to negligible levels, dedicatedsolar neutrino cuts are not implemented for the DSNB
model-independent analysis. For the spectral analysis,

Timing peaks / event / 20 s

17

Exp + constant fit,
= 195.512 ± 6.100= 0.229 ± 0.005
2/dof = 1.266
Constant fit
BG/event/ s = 9.575E­05 ± 5.36E­08
Am/Be data, Center location (2016)

0.020

0.015

0.010

0.005

0.000

0

100

200

300

400

T [ s]

500

600

700

800

FIG. 16. Exponential + constant fit of neutron candidate
timings in the 2016 241 Am/Be sample collected at the centerof the SK tank. The resulting time constant, 195.5 ± 6.1 µs,
is consistent with prior measurements.

however, events with no tagged neutrons are also considered in order to maximize the effective exposure. Whenconsidering these events, we hence apply additional cuts
specifically targeting solar neutrino backgrounds. These
cuts are based on the opening angle θsun between thereconstructed direction of candidate events and the direction of the Sun. The cosine of this angle is peaked at
1 for solar neutrino interaction products, and is smearedby kinematic and reconstruction effects. To estimate theimpact of the latter, an specific observable, the Multiple
Scattering Goodness (MSG) g has been designed for thesolar neutrino analysis [36–39]. Simulated opening angle
distributions for solar neutrino events in different MSG
bins are shown in Fig. 17.
In this analysis we use the method described inRef. [22] and apply cos θsun cuts in the four MSG binsshown in Fig. 17. We estimate the impact of these cutson DSNB signal events by assuming that these events
have a uniformly distributed θsun , and we obtain the associated MSG distributions using the IBD Monte-Carlosimulation. Conversely, to estimate the number of solar
neutrino events above 15.5 MeV, we evaluate the numberof solar events in 13.5-15.5 MeV from data and extrapolate the result to higher energies using the solar neutrino MC simulation from Ref. [77], a procedure alreadyused in the DSNB search at SK-I,II,III [22]. Here, we
can safely assume the cos θsun distribution for non-solar
event to be uniform [39], and therefore approximate thenumber of solar events in 13.5-15.5 MeV by:
solar
Nlow
= Nlow (cos θsun > 0) − Nlow (cos θsun < 0),(5)

using the data passing the noise, spallation, and positron
selection cuts described in Secs V A, 3, and V C. Using

FIG. 17. Opening angle distributions for simulated solar neutrino events for the different MSG bins used in this analysis.
Events with low MSG have wider opening angle distributions
due to direction reconstruction issues.

the solar neutrino MC simulation from Ref. [77], the predicted number of solar neutrino events above 15.5 MeV
will then be:solar
N solar ≈ 0.12Nlow.

(6)

The impact of solar and neutron tagging cuts can thenbe assessed by rescaling N solar by the corresponding efficiencies.
Finally, we evaluate the systematic uncertainties on thesignal efficiency with the solar event cut using the LINAC
calibration results. Specifically, we compute the scaling
factor that allows to minimize the difference between thepredicted and observed MSG distributions for LINAC
events. We then evaluate the impact of this rescaling on
the final efficiency to be of about 1%.

F.

Cut optimization

In this study, we perform two types of analyses. Here
we adopt a common cut optimization procedure for bothanalyses.

1.

DSNB positron candidate selection cuts

We determined the positron candidate selection cutsby comparing the IBD Monte-Carlo simulation with theatmospheric neutrino simulation for the ring clearness,average charge deposit, and Cherenkov angle, and theobserved data for the effective distance deff . Since theimpact of the positron candidate selection cuts depends
weakly on the DSNB model considered, we assumed the
Horiuchi+09 spectrum [16]. As discussed in Section V C,
we estimate systematic uncertainties using results from

18



Erec − 15.5 MeV
deff > max 300, 500 −× 50 cm,
1 MeV
max
Npre < 12,
Ndecay-e = 0,
Lclear < 0.36,
q50 /n50 < 2,38◦ < θC < 50◦ .

uncertainties, computed as described in Sections V B,
V C, and V D, in Table III.

Signal Efficiency [%]the solar neutrino analysis [36–39] and the LINAC calibration [32]. The final cuts are the following:

Unlike positron candidate selection cuts, spallationand neutron tagging cuts need to be optimized in multiple energy regions to account for the important variationsof the spallation and atmospheric background rates in theanalysis windows. While only events with one tagged
neutron are used in the DSNB model-independent analysis, the spectral analysis also considers events with zeroor more than one tagged neutron. Here, we first optimize
spallation and neutron tagging cuts for the signal region
where exactly one tagged neutron is required. Using thethus optimized neutron tagging cuts, we then derive optimal spallation cuts for events with 6= 1 tagged neutron.
a. Events with one tagged neutron:
We devise acut optimization scheme common to both the modelindependent and spectral analyses by simultaneously
adjusting spallation and neutron tagging cuts in Erecbins. For each energy bin, we select events with exactly one tagged neutron and find the spallation and neutron tagging cuts that yield optimal sensitivity under the
background-only hypothesis. We define the sensitivity as
the predicted number of signal events after cuts over the
90% C.L. upper limit on the number of signal events,
computed using the Rolke method [78]. To simplify theinterpretation of the results presented in this paper, weapply the same cuts for both the model-independent andspectral analyses in the 15.5−29.5 MeV reconstructed energy range, where the two analyses overlap. In particular,
we require spallation and solar neutrino backgrounds tobe negligible above 15.5 MeV. The remaining rates ofthese backgrounds can be reliably evaluated by rescaling the number of events observed after noise reduction,
spallation, and positron candidate selection cuts by theneutron mistag rate.
The signal efficiencies for our choice of cuts are shownin the top panel of Fig. 18 for each energy region. For
each set of cuts we also present the associated systematic

SHE+AFT, precut
spallation cut
effwall, pre-/post-activity cuts

60

ring clearness cut
Cherenkov angle cut

40

neutron tagging (tagged-n =1)

20

0

10

15

20

25

30
Erec [MeV]

Signal Efficiency [%]Spallation and neutron tagging cuts

80

charge/hit cut

The efficiencies of these cuts are shown in Fig. 18 as afunction of Erec and vary between 70 and 85%. Overall,positron candidate selection cuts allow to remove 80%of the atmospheric neutrino backgrounds, and bring contamination from natural radioactivity to negligible levels.

2.

100

100

80

60

SHE+AFT, precut

40

spallation cut
effwall, pre-/post-activity cuts
ring clearness cut
charge/hit cut

20

Cherenkov angle cut
neutron tagging (tagged-n =0 or >1)
solar neutrino cut

0

16

18

20

22

24

26

28

30
Erec [MeV]

FIG. 18. DSNB signal efficiencies for the signal regions where
one tagged neutron (top) and 6= 1 tagged neutrons (bottom)
are required, as started after the SHE+AFT trigger requirement and the noise reduction cuts. While the one tagged
neutron region is common to both the model-independent andspectral analyses above 15.5 MeV, the 6= 1 tagged neutron region is used only in the spectral analysis. Here, we applied
noise reduction, spallation, third reduction, neutron tagging,and solar cuts sequentially and show cumulative efficiencies
at each stage. The result in the final bin above 29.5 MeV
corresponds to the efficiency for the range through 29.5 to79.5 MeV, as calculated assuming the Horiuchi+09 6 MeV
spectrum [16]. Note also that efficiencies for the solar neutrino cut in the bottom are shown in 1-MeV bins.

b. Events with zero or more than one tagged neutrons:
As previously discussed, the events rejected by
the neutron tagging cuts obtained above can still be considered in the spectral analysis. In the 23.5−79.5 MeV
region, that covers most of the spectral analysis window, spallation and solar neutrino backgrounds are negligible and only noise reduction and positron candidate
selection cuts need to be applied to these events. In
the 19.5−23.5 MeV region, solar neutrino backgrounds

19
TABLE III. Systematic uncertainties for the reduction cuts,
computed as described in Sections V B, V C, and V D. The
systematic uncertainties for the AFT trigger, noise reduction,pre- and post-activity cuts are negligible. The solar neutrino
cut is used only for the spectral analysis, when events have
6=1 tagged neutrons.
Cut
Relative systematic error
Spallation cut
0.1%
deff cut
0.1%
Lclear cut
0.2%
q50 /n50 cut
1.2%
θC cut0.7%
Neutron tagging
12.5%
Solar neutrino cut
1.0%

are still negligible and spallation backgrounds are mostly
linked to short-lived isotopes, that are easy to reject using the cuts from Section V B. For these energies, wetherefore apply the spallation cut derived above, whose
efficiency is close to 90%. In the 15.5−19.5 MeV range,
however, both solar neutrino and spallation backgrounds
largely dominate and need to be reduced using the targeted reduction strategies described in Sections V B and
V E. In this analysis, we remove the solar neutrino backgrounds using the cuts in Section V E with the cut points
from Ref. [22]. The cut points and the corresponding
signal efficiencies are shown in Table IV for the different
multiple scattering goodness and reconstructed kinetic
electron energy regions. The remaining number of solar
neutrino events after all cuts is less than 2 for the entire
energy region.
TABLE IV. Cuts on cos θsun (upper limits), and signal efficiencies for SK-IV. These cuts are only applied to events with
6= 1 tagged neutron, in the spectral analysis. They depend on
the value of the multiple scattering goodness, g, and on thereconstructed kinetic electron energy. They are the same as
the ones applied for SK-I and III in Ref. [22]. The signal efficiency for SK-IV is within 1% of the one obtained in Ref. [22]for SK-I and III, which is similar to the associated systematic
uncertainty.
Erec [MeV] 15.5−16.5 16.5−17.5 17.5−18.5 18.5−19.5
g < 0.4
0.05
0.35
0.45
0.93
0.4 < g < 0.5
0.39
0.61
0.77
0.93
0.5 < g < 0.6
0.59
0.73
0.81
0.93
g > 0.6
0.73
0.79
0.91
0.93
Efficiency
73.1%
82.2%
88.3%
96.6%

In contrast, the impact of spallation cuts is difficult toevaluate in the 15.5−19.5 MeV range . To optimize spallation cuts for the events considered here, we therefore
proceed as follows. First, we bring contributions fromshort-lived isotopes —with O(0.01) sec half-lives— to anegligible level, using the cuts shown. in Appendix A. As
can be seen in Fig. 3, the remaining event sample will bedominated by decays from 8 B and 8 Li, that have a halflife close to 1 sec. We then optimize the spallation like-lihood cuts using the method described in Appendix A 4
to compute the remaining
amount of 8 B and 8 Li decays,√
and maximize S/ B where B and S are the numbersof spallation and non-spallation events respectively. The
efficiencies of the optimal cuts are of 65% and 63% in
15.5−17.5 and 17.5−19.5 MeV, for a 8 B+8 Li remaining
fraction of about 2.5%. The corresponding remaining
number of spallation events, as well as the performanceof this cut will be discussed and compared to results fromthe previous SK DSNB analyses in Section VII E.

G.

Selected events in the final data samples

Here we discuss the behavior of the observed data after the cuts. In particular, the impact of spallation, deff ,and the other positron candidate selection cuts on datais shown in Fig. 19. As shown in this figure, for reconstructed energies below 20 MeV contributions fromspallation and radioactivity near the wall are significant.
To validate our background modeling and neutron tagging procedure we also compare the observed data to thebackground expectations after the noise reduction, multiple spallation cuts, DSNB positron candidate selection,and neutron tagging cuts. To simplify the interpretationof the data, we impose the same neutron tagging cut,with a 20% signal efficiency, for all energies. While mostspallation cuts are removed in order to focus on neutron
tagging, multiple spallation cuts are conserved in order toavoid contaminating the neutron detection window withisotope decays —an effect not accounted in the mistag
rate predictions. The predicted and observed event spectra are shown in Fig. 20 and agree within uncertainties.
After the cuts derived in Section V F are applied, we find 102 events with exactly one neutron in
7.5−79.5 MeV. For these events, the distribution of thetime difference between the neutron and the prompt
event is shown in Fig. 21. Fitting this distribution by an
exponential plus a constant, as is done in the 241 Am/Be
calibration, yields a time constant of τ = 154 ± 104 µs,
which is compatible with the expected neutron capture
time in water. In the 15.5−79.5 MeV region, that will
be used for the spectral analysis, we also find 248 eventswith zero neutrons and 9 events with more than one neutron. The observation times and reconstructed vertices
for these events are shown in Figs. 22 and 23. No event
time cluster has been observed and the time-dependenceof the event rate is consistent with the livetime variations over the SK-IV period. Events are also uniformly
distributed all over the FV, irrespective of the numberof neutrons. The space and time distributions of theselected events are therefore consistent with what is expected for DSNB candidates.

Events/2-MeV

20

6

10

Noise reduction
Spallation
deff
Positron ID

5

10

4

10

3

10

2

10

10
10

15

20

25

30
Erec [MeV]Events/2-MeV

FIG. 19. Data spectrum after noise reduction, spallation, deff ,and all other positron candidate selection cuts. Here the cuts
are applied in the order shown above.

700

A.

SK-IV Data
Atmospheric-ν (non-NCQE)
Atmospheric-ν (NCQE)
9
Spallation Li
Reactor-ν
Accidental coincidence

600
500
400
300
200
100
0

8

10

12

14

16

18

20

22

24

26

28
Erec [MeV]

FIG. 20. Predicted and observed spectra for events passing
noise reduction, multiple spallation, DSNB positron candidate selection, and neutron tagging cuts. Here we apply auniform neutron tagging cut with a 20% signal efficiency for
all energies. In the absence of most spallation cuts, the dominant backgrounds are accidental coincidences and 9 Li decays.
The systematic uncertainties are shown by the hatched lines.

VI.

FIG. 21. Time differences between reconstructed neutron captures and their corresponding prompt events (positron candidates). Here we show the data for one tagged neutron. The
solid orange curve shows the fit result by an exponential plus
a constant, giving a time constant of τ = 154 ± 104 µs.

DSNB MODEL-INDEPENDENT ANALYSIS

In this section we perform a DSNB model-independent
search. Specifically, we divide the 7.5−29.5 MeV reconstructed energy range into 2-MeV bins and, in each ofthem, search for DSNB IBDs. In what follows, we describe how to estimate the various backgrounds encountered in this analysis, model their energy dependence,and estimate the corresponding uncertainties.

1.

Background estimation
Atmospheric neutrinos

In the high end of the analysis window (Erec >19.5 MeV), the atmospheric neutrino spectrum is dominated by decay of invisible muons or pions from CC or
NC interactions, forming the well-known Michel spectrum. Systematic uncertainties in this region will therefore only affect the total number of predicted atmospheric
neutrino events. Here, we estimate this number using the
sideband region of 29.5 < Erec < 49.5 MeV, that has asimilar background composition. The number of eventsin this sideband is 88.8% of the atmospheric neutrinoMonte-Carlo prediction, and we therefore use this factor
to rescale the non-NCQE backgrounds. The associateduncertainty is due to statistical fluctuations in the number of events in the sideband and is around 19%. The
resulting spectrum is labeled as “non-NCQE” in Fig. 24.
Below 19.5 MeV, the dominant atmospheric neutrino
background arises from NCQE interactions. The crosssection measurements from T2K [48] are used to estimate this background. In particular, the MonteCarlo simulation predictions of neutrino and antineutrino NCQE interactions are renormalized using thescaling factors from this measurement; fν-NCQE =0.80 ± 0.08(stat.)+0.24−0.18 (syst.) and fν̄-NCQE = 1.11 ±0.18(stat.)+0.29(syst.),
respectively.
Note that this
−0.22
rescaling provides an inclusive estimation of the NCQE
and NC 2p2h interactions. The errors of these factors
are taken as cross section uncertainties.
Two additional flux uncertainties are considered for
the NCQE background: the atmospheric neutrino flux
uncertainty and the uncertainty due to the flux difference

21

FIG. 22. Observation times for events with one neutron (top),
and with zero or more than one neutron (bottom) after thecuts. Overlaid is the expected time distribution, where variations of the data-taking rate over the SK-IV period are taken
into account.

between the T2K beam and atmospheric neutrinos. For
the former, a 15% uncertainty is taken from Refs. [49, 50].
The latter uncertainty arises from the fact that the scaling factors above are measured for the T2K fluxes while
the current focus is the atmospheric neutrino flux andthen the effect of cross section model uncertainties is different. To estimate this uncertainty, the ratios of fluxaveraged events in the T2K beam and atmospheric neutrino fluxes for different cross section models are calculated. Here six models are considered: Spectral Function(SF) [79, 80], Relativistic Mean Field (RMF) [81–84], Superscaling (SuSA) [85, 86], Relativistic Green’s Function(RGF) [84, 87–89] with two functional forms for the potential (EDAI, Democratic), and Relativistic Plane Wave
Impulse Approximation (RPWIA) [84]. The maximum
differences are taken as systematic errors, 5% for neutrinos and 7% for antineutrinos.
There are two systematic error sources related to neutron tagging for the NCQE background: tagging effi-

ciency and neutron multiplicity. For the uncertainty on
the tagging efficiency, 12.5% is employed in this analysis
as the maximum difference between the measured andpredicted efficiencies in the 241 Am/Be calibration. The
neutron tagging efficiency will also depend on the distance between the neutron emission and capture vertices.
Indeed, neutrons produced in atmospheric interactions
can travel up to a few meters and are hence more difficult to identify than neutrons produced in DSNB IBDs.
The effect of the uncertainties in modeling the neutron
travel distance on the efficiency is of around 7% [47].
In this analysis, the number of tagged neutrons is required to be one, therefore the model uncertainty affecting the neutron multiplicity has to be estimated. Herethe recent measurement of the neutron multiplicity after (mostly CC) interactions of neutrinos from the T2K
beam inside the SK tank is used [75]. This measurement allows to compare observations with Monte-Carlosimulation predictions for neutron multiplicity as a function of the reconstructed squared momentum transferredto the nucleus, Q2 , and could therefore be adjusted to
study the NCQE background. We incorporate the maximum difference between observed and predicted neutron
multiplicity in T2K for each Q2 range, and renormalize
the present Monte-Carlo simulation prediction to accommodate this gap. This procedure gives variances of approximately 40% in the number of events with the tagged
neutron being one, both for neutrinos and antineutrinos.
The NCQE spectral shape is sensitive to the γ emission
model, affecting the current differential limit extraction.
It is complicated to estimate the associated systematic
uncertainty, as γ rays are produced at every stage of the
NCQE process, from the primary neutrino interaction to
the secondary nuclear reaction. The modeling of these
γ emission is in particular a probable cause of the currently observed discrepancy between the observed andpredicted Cherenkov angle distributions using neutrinos
from the T2K beam [48]. One notable effect of thismismodeling is a smearing of the reconstructed energy
analogous to a resolution effect. To evaluate the associated systematic uncertainty we model this smearing byconvolving the original spectrum with a Gaussian distribution and using the energy-dependent ratio between
the smeared and unsmeared spectra to renormalize theMonte-Carlo events. We then compute the number ofrenormalized events in each energy bin after the reduction steps described in Section V. The difference between
this result and the nominal prediction is then taken as
the systematic uncertainty associated to the γ emission
modeling. For this study, we set the smearing parameterσ to 3 MeV for all energies; this value allows to cover
the discrepancy observed between the predicted and observed θC distributions in the T2K measurements mentioned above. The resulting uncertainty ranges from 30%
to 60% in the 7.5−29.5 MeV range.
Finally, we combine all the systematic uncertainties
described above by adding them in quadrature. The total
systematic uncertainty on the NCQE background then

22

FIG. 23. Top (top) and transverse (bottom) views of reconstructed vertices for events with one neutron (left), and with zeroor more than one neutron (right) after the cuts. Events with more than one neutron are shown as red crosses. The solid anddashed gray lines show the ID and the FV, respectively.

varies from 60% to 80% for low to high energies. The
predicted NCQE spectrum with the T2K scaling is shownin Fig. 24.
2.

Spallation 9 Li

The background rate of spallation 9 Li can be estimatedby the following formula:
B9 Li = R9 Li × Tlive × 22.5 kton
9×Br[ Li → β + n] × fwindow × reduc ,(7)

where R9 Li is the production 9 Li rate measured at SK(0.86 ± 0.12(stat.) ± 0.15(syst.) kton−1 ·day−1 ) [59], Tlive
is the operational livetime (2075.3 and 2970.1 days below and above Erec = 9.5 MeV, respectively), Br[9 Li →β + n] = 0.508 is the branching ratio for β+n decay,
fwindow is the fraction of the 9 Li decay energy spectrum
above the search energy threshold, and reduc is the reduction efficiency. We estimate fwindow using the MonteCarlo simulation elaborated for IBD events, renormalized
to fit the 9 Li spectrum. This same simulation can be used

to estimate the efficiencies for the noise and positron reductions as well as the neutron tagging, that are the same
as for the DSNB signal. The spallation cut efficiency isestimated using the procedure outlined in Section V. The
dominant systematic uncertainty on the number of 9 Li
events in each energy bin after cuts arises from estimating the impact of spallation cuts. Indeed, as discussed
in Appendix A 4, our method assumes that the time difference between an isotope decay and its parent muonis the only isotope-dependent spallation observable. In
this study, we account for the impact of other spallation
observables by defining a 50% uncertainty on the number of 9 Li events. The uncertainty on the event rate istaken from the previous SK study [59] as 22%. Other uncertainties come from the reduction, especially from thespallation cut and neutron tagging, that is 20% in total.
The total uncertainty assigned for the 9 Li backgroundestimation is therefore 60%.

3.

Reactor neutrinos

Reactor neutrino backgrounds are estimated by renormalizing the IBD Monte-Carlo simulation. They popu-

23
lates only in the lowest energy bin, as seen in Fig. 24. In
this analysis, a conservative 100% uncertainty is assigned
to their estimated rate.

4.

Accidental coincidences

A prompt SHE signal accidentally paired with a subsequent dark noise fluctuation or a low energy radioactive
decay can look like an IBD signature. Since accidental
pairing can occur for both signal and background events,
the resulting background can be readily estimated bycomputing the number of events left after noise reducdata
tion, spallation, and positron reduction cuts Npre-ntag,
and rescaling it by the neutron mistag rate:

data
Bacc = mis × Npre-ntag.

(8)

The associated systematic errors are highly cutdependent and arise from the time dependence of the
PMT dark noise and natural radioactivity. In order toevaluate its effect on the background rejection, we separate the random trigger data mentioned in section III B
over all the SK-IV period into 10 time bins of about eightmonths each. Then we evaluate the BDT mistag rate foreach bin separately and, for a given signal efficiency, takethe standard deviation from the average mistag rate assystematic error. The typical size of uncertainty is a few
to 30%, in the region of interest for the DSNB analyses,
as shown in Fig. 15.

B.

φlimit
=
90

limit
N90
,
t · Np · σ̄IBD · sig

(9)

where t is the operational livetime [sec], Np is the number of free protons in the SK’s FV, σ̄IBD is the IBD
cross section [10−41 cm2 ] at a mean neutrino energy inthe corresponding region (Ēν ), and sig is the signal efficiency. Note that the neutrino energy is obtained asEν = Erec + 1.8 MeV in IBD. For the expected sensitivity, the same procedure is applied while the numberof observed event is replaced with the number of nominal background prediction. Here statistical uncertainties
on the backgrounds are considered instead. The expected
sensitivities and observed upper limits from this work aresummarized in Table V. These results are also comparedwith the previously published results and theoretical predictions in Fig. 25. Note that the SK-I,II,III limits presented in this figure have been obtained using a different
methodology, more similar to the one used for the spectral analysis in Section VII. In addition, the SK-I,II,III
analysis not only used a higher energy threshold but also
did not involve neutron tagging [22, 71]. Hence, we do
not combine the SK-IV result with the SK-I,II,III one forthe model-independent differential limit. The sensitivities obtained with this analysis are the world’s tightest
above neutrino energies of 11.3 MeV. The current limit
disfavors the Totani+95 model [8] and the most optimistic predictions of the Kaplinghat+00 model [9], and isreaching close to several other model predictions. Due tonot only a higher exposure but also higher cut efficiencies
and more precise background estimation, this new analysis considerably improves on the previous SK-IV DSNB
model-independent search [71].

Differential upper limit

The observed data and expected background spectra
are shown in Fig. 24. The corresponding p-values arecalculated for each energy bin by performing pseudo experiments as follows. We vary the number of expectedbackground events randomly based on Gaussian distributions with their widths being the 1σ systematic uncertainties from each background source. We then calculate
p-values from the resulting distributions and the observed
number of events in each bin, as shown in Table V. Herethe most significant bins are found to be at 2σ level. We
therefore conclude that no significant excess is observedin the data over the background prediction in any energybin, and place upper limits on the extraterrestrial ν̄e flux
using the pseudo experiments above. Here we obtain the
limit)
90% C.L. upper bound on the number of signal (N90
as an excess of the observation over the background expectation, by varying the number of observed events ineach energy bin by their statistical uncertainties as wellas varying the number of expected background events by
their systematic uncertainties from each source. The 90%
C.L. upper limit on the ν̄e flux is then calculated as:VII.

SPECTRAL FITTING

In this section, we derive model-dependent limits on
the DSNB flux by fitting signal and background spectral
shapes to the observed data in the 15.5−79.5 MeV range,
and combine the results with the ones from previous SKphases [22] to achieve a 22.5 × 5823 kton·day exposure.
While atmospheric neutrino backgrounds, notably fromdecays of invisible muons and pions, will dominate overmost of the analysis window, we also account for possible
residual spallation in 15.5−19.5 MeV.

A.

Signal and sideband regions

In order to evaluate the number of atmospheric neutrino background events, we define three regions of parameter space based on the reconstructed Cherenkov angle: one signal region with θC ∈ [38, 50]◦ , that will contain most of the signal and the irreducible backgrounds,
and two sidebands with θC ∈ [20, 38]◦ and θC ∈ [78, 90]◦ .

24
TABLE V. Summary on the 90% C.L. expected sensitivities and observed upper limits as well as the corresponding p-values
in each electron antineutrino energy bin (Eν = Erec + 1.8 MeV).

Events/2-MeV

Eν [MeV]
9.3−11.3
11.3−13.3
13.3−15.3
15.3−17.3
17.3−19.3
19.3−21.3
21.3−23.3
23.3−25.3
25.3−27.3
27.3−29.3
29.3−31.3

25

Expected [cm−2 sec−1 MeV−1 ]4.44 × 101
1.14 × 101
4.17 × 100
1.87 × 100
8.48 × 10−1
4.64 × 10−1
3.28 × 10−1
2.11 × 10−1
2.13 × 10−1
1.98 × 10−1
1.50 × 10−1

SK-IV Data
Atmospheric-ν (non-NCQE)
Atmospheric-ν (NCQE)
9

Spallation Li

20

Reactor-ν
Accidental coincidence
DSNB (Horiuchi+09 6-MeV, Maximum)

15

10

5

0

10

20

30

40

50

60

70

Observed [cm−2 sec−1 MeV−1 ]3.71 × 101
2.04 × 101
9.34 × 100
3.29 × 100
5.08 × 10−1
6.84 × 10−1
1.27 × 10−1
3.75 × 10−1
7.77 × 10−2
2.42 × 10−1
7.09 × 10−2p-value
0.346
0.886
0.938
0.830
0.243
0.686
0.073
0.597
0.051
0.605
0.126

The low Cherenkov angle region will be populated withmostly atmospheric backgrounds involving visible muonsand pions while the high angle region will be mostly populated by NCQE atmospheric neutrino events with multiple γ rays. Finally, we separate events with exactly
one identified neutron from the others, thus defining an“IBD-like” and a “non IBD-like” region. Note that due
to the low efficiencies of the neutron tagging cuts the non
IBD-like region is expected to contain a sizable amountof signal. Our analysis will hence involve six regions ofparameter space: two signal regions with intermediate
values of the Cherenkov angle, and four sidebands withlow and high Cherenkov angle values, as summarized inTable VI.

Events/2-MeV

Erec [MeV]

SK-IV Data
Atmospheric-ν (non-NCQE)
Atmospheric-ν (NCQE)
9

Spallation Li
Reactor-ν

10

Accidental coincidence
DSNB (Horiuchi+09 6-MeV, Maximum)

TABLE VI. Overview of the regions used in the spectral analysis. We split the parameter space according to the reconstructed Cherenkov angle and the number of tagged neutrons.
Regions with small and large Cherenkov angles are dominatedby the interactions producing visible muons and pions and by
NCQE interactions, respectively. We assign numbers for eachregion.
θC
Ntagged-n20−38◦ 38−50◦ 78−90◦

1
6=1

I
IV

II
V

B.

Spectral shape fitting

III
VI

1

10

20

30

40

50

60

70
Erec [MeV]

FIG. 24. Reconstructed energy spectra after data reductions
for the background expectation and observation including thesignal window and the sideband region with a linear (top)and a log (bottom) scale for the vertical axis. Color-filled histograms correspond to each background source with hatchedlines in the signal window representing the total systematic
uncertainty. The red dashed line represents a DSNB signal
expectation from the Horiuchi+09 model [16] shown only forthe signal window.

We perform a simultaneous fitting of the signal andbackground spectra to the observed data in all six regionsof parameter space defined in Table VI using an extended
maximum likelihood method. Performing this type ofanalysis requires knowing the shapes of the signal andbackground spectra in each of the regions. While the signal spectrum can be reliably predicted by the IBD MonteCarlo simulation for a given DSNB model, the treatmentof the atmospheric neutrino and spallation backgrounds
is more complex. For this study, we follow the methodνe Flux Upper Limit [/cm2/sec/MeV]

25

SK-IV 2970 days, Observed 90% C.L. (This work)
SK-IV 2970 days, Expected 90% C.L. (This work)
SK-IV 960 days (2015)

2

SK-I/II/III 2853 days (2012)10

KamLAND 4529 days (2021)
DSNB Theoretical Predictions

10

1

10−1
10

15

20

25

30

νe Energy [MeV]
FIG. 25. The 90% C.L. expected and observed upper limits on the extraterrestrial electron antineutrino flux from the present
work, in comparison with previously published results from SK [22, 23] and KamLAND [25] and DSNB theoretical predictions
from Fig. 1 (in gray). The upper limit from Ref. [22] (blue) has been derived in Ref. [23].

described in Ref. [22] and divide these backgrounds into
the following five categories.
a. Invisible muons and pions: this category regroups events with electrons produced by the decays ofinvisible muons and pions. The energy distribution ofthese electrons follows a Michel spectrum, whose shapeis independent of the Cherenkov angle and the neutron
multiplicity and will therefore be the same across all six
regions. For this study, we estimate the shape of theMichel spectrum at SK directly from data, using a sample of electrons produced by cosmic ray muon decays(a similar sample for atmospheric and accelerator neutrino oscillation and proton decay analyses is describedin Refs. [90, 91]). We then use the atmospheric neutrino Monte-Carlo simulation to compute the fractionsof background events in the different signal and background regions. Since electrons cannot be distinguishedfrom positrons at SK this background dominates in thesignal regions II and V and are negligible everywhere else.

b. νe CC interactions: in this category we find backgrounds arising from CC interactions of electron neutrinos and antineutrinos, with no visible muons and pions
in the final state. Their contributions will dominate in
the signal regions II and V above 50 MeV. We estimate
the associated spectral shapes in all regions using theatmospheric neutrino Monte-Carlo simulation. Similarly
to Michel electrons, this background is negligible outside
the regions II and V.
c. µ/π-producing interactions: visible muons andpions will be associated with low Cherenkov angles, asthese particles are significantly heavier than electrons.
The associated background will therefore dominate in
the low Cherenkov angle regions I and IV, and, after
positron candidate selection cuts, will be negligible in thesignal regions. We extract the associated spectral shapesby considering an atmospheric Monte-Carlo sample with
only CC interactions, visible muons and pions, and noelectrons.

26d. NCQE interactions: unlike the other types of atmospheric neutrino events the NCQE spectrum peaks at
low energies, similarly to the DSNB spectrum. Here we
define spectral shapes using simulated NCQE events withno electrons. This background, often involving multiple
γ rays, will dominate in the high Cherenkov angle regions
III and VI, but will also contribute to all other regions.
e. Spallation backgrounds: these backgrounds are
negligible in the one-tagged-neutron regions I, II, andIII after the cuts. However, in the other regions, especially in signal region V, residual spallation backgrounds
from 8 Li, 8 B, and 9 C decays could remain. Since 9 C has
a shorter half-life and is produced closer to the muontrack than the other two isotopes, it is expected to be
subdominant. To model the spallation spectrum and its
variability we parameterize this spectrum and vary its
isotope composition as follows:
Sspall (E) =f0 S9 C (E)(10)+ (1 − f0 ) (f1 S8 B (E) + (1 − f1 )S8 Li (E)) ,with f0 ∈ [0, 0.5] and f1 ∈ [0, 1]
where Si (E) refers to the decay spectrum for a given
isotope i, normalized to one in the spectral analysis
energy window, and f0 and f1 are the isotope fractions.
Here, we constrain 9 C to always be less abundant thanthe other two isotopes after spallation cuts are applied.
In what follows, we will take the nominal spallation
spectrum to be the average between the steepest and theleast steep spectra and parameterize possible deviations
from this hypothesis as systematic uncertainties. We
evaluate the β spectra using external measurements [92]for 8 B decays and GEANT4 for 8 Li and 9 C. As with 9 Li,
we then use these spectra to renormalize events in the
IBD Monte-Carlo simulation and thus model detector
resolution effects.
Although the first four background categories do notinclude all possible types of atmospheric neutrino interactions, their associated spectra can be linearly combined
to model other background categories. Pion-producing
NC backgrounds can notably be treated as a superposition of NCQE and decay electron backgrounds. The
five categories described here hence allow to map all thebackground configurations relevant to this analysis.
For each background category we define probability
distribution functions (PDFs) which model the corresponding energy spectra in each of the six Cherenkov angle and neutron regions. These PDFs are normalized to 1
across all regions. We apply a similar procedure to obtainsignal spectral PDFs associated to different DSNB models, using the IBD Monte-Carlo simulation. For Nevents
observed events with energies {E1 , ..., ENevents } we thenextract the most probable numbers of signal and background events by performing an extended unbinned maximum likelihood fit, using the following likelihood function:

L({Nj }) = e−P5

j=0 NjNY
5events X
i=1(r)

Nj PDFj (Ei ). (11)j=0Here, the index j = 0, ..., 5 refers to the signal and the fivebackground categories described above. Then {Nj } designates the numbers of events for the signal and the different background types, across all six Cherenkov angle(r)and tagged neutron regions. For each event, PDFj (Ei )is the PDF for the signal or background category j in theregion r (for the regions defined in Table VI).
C.

Systematic uncertainties

One considerable advantage of using the likelihood defined in Equation 11 is that the sizes of the different background contributions can be treated as nuisance parameters. This likelihood, however, does not account for theconsiderable uncertainties on the spallation and atmospheric background spectral shapes in the energy windowof the DSNB analysis. In this section, we outline how
to incorporate systematic uncertainties on these shapes
as well as on detector modeling and predicted cut efficiencies. We account for systematic uncertainties on thespallation and atmospheric neutrino background spectral shapes, energy reconstruction, and reduction efficiencies by convolving the likelihood introduced in Equation 11 with suitable probability distributions. Here weneglect the systematic uncertainties associated with theµ/π background. The amounts of these background in
the signal regions II and V are indeed constrained to be
negligible by observations from the low angle regions I
and IV, a result that is robust under large variations of
the assumed spectral shapes in these regions. We also neglect the systematic uncertainties on the spectral shapeof the background associated with invisible muon andpion decays since this shape has been directly extracted
from a large sample of SK-IV data. We then consider
five sources of systematic uncertainties: energy scale andresolution, spectral shape predictions for the spallation,
NCQE, and νe CC backgrounds, and reduction cut efficiencies.
1.

Energy scale and energy resolution

To estimate the impact of energy scale and energy
resolution uncertainties on the signal and atmosphericbackground spectra, we follow the procedure described inRef. [22]. The energy scale and energy resolution uncertainties are obtained from the solar neutrino analysis andincreased to 2% and 3% respectively to account for thehigher energies studied here. We model the effect of these
uncertainties by distorting the signal and atmosphericbackground PDFs in Equation 11 using the energy resolution functions from Refs. [36–39]. We parameterize the27amount of distorsion applied to these PDFs using thevariable  such that  = 1 represents a 1σ deviation in
the energy scale and resolution. The final likelihood canthen be expressed as:

L0 =

Z +∞

2

e− /2 L() d,(12)vary in the two signal regions. We then integrate the distorted PDF using pseudo-Gaussian weights. A detailed
description of this procedure is provided in Ref. [22].
At this step of the analysis, we assume the spectral
shape distortions in the =1 and 6=1 tagged neutron regions to be fully correlated and model the effect of neutron multiplicity uncertainties separately.

−∞

where L() represents the likelihood from Equation 11
with the PDFs distorted by an amount . These uncertainties have a negligible impact on the final result andare therefore not taken into account in the final fit.
While energy scale and resolution uncertainties on signal and atmospheric background spectra can be neglectedin this analysis, the spallation background spectrum isparticularly steep and hence more sensitive to the modeling of energy reconstruction. In this analysis, we will
therefore account for the associated modeling uncertainties when estimating spallation spectral shape systematic
uncertainty.

2.

Spallation background spectral shapes

As discussed above, we consider two sources of systematic uncertainties here: energy scale and resolution
uncertainties, and isotope composition. We estimate thelatter by varying the 8 B, 8 Li, and 9 C fractions as shown inEquation 3. Then, we distort the steepest and least steep
spectra found using this method in order to account forthe 2% and 3% energy scale and resolution uncertainties
discussed above. We take these distorted spectra as our±1σ variations from the nominal spallation spectrum.
Finally, we parameterize the total spectral distorsion associated with the systematic uncertainty by using a thirdorder polynomial P3 such that:
PDFnew = PDFold × N [1 + P3 (E)] ,(13)

where  is expressed in units of σ. More details about
these estimates are given in Appendix D. We then convolve the likelihood L from Equation 11 with a Gaussianprior on :

L0 =

Z 2

2

L()e− 2 d.

(14)

−2

3.

Atmospheric νe CC spectral shapes

The shape of the νe CC background for electron energies below about 100 MeV has not been studied in-depth
to date, although Monte-Carlo simulations predict a linear increase in the number of events with the electronenergy. In our analysis, we model this spectral shape uncertainty by allowing the slope of the νe CC spectrum to4.

Relative normalization for NCQE events

The Cherenkov angle distribution of atmospheric
NCQE events is associated with large uncertainties stemming from the modeling of γ ray emission. We incorporate this uncertainty in the current analysis by allowing
the relative normalization of the NCQE PDFs between
the different Cherenkov angle regions to vary. Following the procedure outlined in Ref. [22], we define a 1σvariation of this normalization to correspond to a 100%change of the number of NCQE events in the intermediate Cherenkov angle region. Additionally, in order for theoverall normalization of the NCQE PDFs across regions
to remain constant, we reweight these PDFs as follows:
NCQE
PDF’NCQE
med, n = PDFmed, n × (1 + ),

(15)

NCQE
PDF’NCQE
high, n = PDFhigh, n × (1 − fhigh ),

(16)

NCQE
PDF’NCQE
low, n = PDFlow, n × (1 − flow ),(17)

where the “low”, “med”, and “high” indices refer tothe low (20−38◦ ), medium (38−50◦ ), and high (78−90◦ )Cherenkov angle ranges from Table VI and n = 0, 1 forthe tagged neutron regions respectively. flow , fmed , fhigh
are scaling factors introduced to preserve the overall normalization of the PDFs. The  parameterizes the PDF
distortion, in units of σ. We finally inject these distorted
PDFs into Equation 11 and convolve the new likelihood
with a weighted Gaussian prior, as with νe CC background uncertainties.

5.

Uncertainty on the neutron multiplicity

The fractions of signal and background events observedin the =1 and 6=1 tagged neutron regions depend on boththe neutron tagging efficiency and the true neutron multiplicity of the observed events. As shown in the 241 Am/Be
calibration study, the efficiency of the neutron tagging
cut for the DSNB spectral analysis is associated with asystematic uncertainty of around 6%. Conversely, the actual neutron multiplicity for atmospheric neutrino interactions is poorly known. As discussed in Section VI A,
neutron multiplicity uncertainties of up to ∼40% havebeen measured for reconstructed energies above 100 MeV
using the T2K neutrino beam [75]. In what follows, wetherefore consider a 40% systematic uncertainty on theneutron tagging.

28
We model the effect of neutron multiplicity uncertainties by varying the relative normalization of the PDFs
in the =1 and 6=1 tagged neutron regions for the νe CC,
NCQE, and invisible muon and pion backgrounds. As
mentioned at the beginning of this section we neglect thesystematic uncertainty associated with µ/π backgrounds,
which are associated with negligible contributions in thesignal region. As discussed above, we set the 1σ variation of the fraction of events with one tagged neutron toa conservative value of 40%. We then distort the PDFs
in the =1 and 6=1 tagged neutron regions for each background category as follows:
IBD
PDFIBDj,new = PDFj,old × (1 + αj ),(18)

PDFnon-IBD= PDFnon-IBD
× (1 − αj f ),
j,new
j,old

(19)

where αj is the relative 1σ systematic uncertainty forbackground category j,  measures the magnitude of the
distortion in units of σ and f preserves the overall normalization of the PDFs. Finally, we incorporate this uncertainty into the final likelihood by integrating the initial
likelihood similarly to Equation 14.

6.

Combining the background uncertainties

The three categories of systematic uncertainties onbackgrounds described above —spallation and νe CC
spectral shapes, NCQE Cherenkov angle distribution,and neutron multiplicity— are incorporated into Equation 11 using the modified likelihood defined in Equation 14 and the modified PDFs defined in Ref. [22] aswell as in Eqs 17 and 19 with distorsion parameters νe CC ,
NCQE , and n respectively. The final likelihood Lsyst canthen be written as:

through the detector, R, requires rescaling Ns by the average value of the signal efficiency over all analysis regions
and their associated energy windows, sig . We therefore
need to incorporate the systematic uncertainty on sig
into Lopt , following a procedure similar to the treatmentof the background systematic uncertainty. Since the neutron tagging efficiency only affects the relative normalization between regions, the main sources of uncertainties
on sig will arise from noise reduction and positron candidate selection cuts. These uncertainties are summarizedin Table III. We also incorporate the uncertainties associated to the fiducial volume cut, inverse beta decay crosssection, and livetime calculation, that have been takenfrom the solar neutrino analyses and are summarized inTable VII. Adding all these uncertainties in quadrature,
we then rewrite the likelihood as:

0

Z 1

L (R) =

Lopt (R)e−(−sig )2
2σ 2
sig

d,

(21)

0

where σsig is the 1σ systematic uncertainty on sig .
TABLE VII. Uncertainties associated with the fiducial volume cut, the Strumia-Vissani IBD cross section calculation,and the livetime calculation for SK-I to IV. The uncertainties
for SK-I,II,III have been taken from the SK solar neutrino
analyses [36–39]. For SK-IV the uncertainties have been obtained from the IBD Monte-Carlo simulation.
SK phaseI
II
III IV
Fiducial volume 1.3% 1.1% 1.0% 1.5%
Cross section 1.0% 1.0% 1.0% 1.0%
Livetime
0.1% 0.1% 0.1% 0.1%

Systematic uncertainties associated with reduction efficiencies weaken the final limits by around 1%. Similarly,
while the uncertainties associated with spallation and atmospheric neutrino backgrounds are of order 100%, their
Z 3
Z 3
Z 3
Z 2
overall effect on the final limits is of about 5%. The limLsyst ({Nj }) =dνe CCdNCQEdndspallited impact of spectral shape uncertainties is due to the
−1
−1
−2
−2
large dominance of backgrounds associated with invisible
× L({Nj }, νe CC , NCQE , n , spall )(20)
muon and pion decays, whose spectrum is well-known.
2
2
The final analysis will hence be largely dominated by− spall− 2ne 2e, statistical uncertainties, which will allow to easily com× G(νe CC ) × G(NCQE ) × √ × √ππbine SK-IV observations with results from previous SKphases.
where L({Nj }, νe CC , NCQE , n , spall ) is the likelihood
defined in Equation 11 with distorted PDFs and G is
the weighted Gaussian introduced in Ref [22].
D. Combination with SK-I,II,III

7.

Energy-independent efficiency systematic error

Maximizing the likelihood Lsyst defined in Equation 20
over the numbers of events in all five background categories results in a marginalized likelihood Lopt (Ns ) where
Ns is the number of signal events. Using this likelihood toconstrain the total number of DSNB events that passed

The methodology described above, based on an extended maximum likelihood fit over signal and sideband
regions, is similar to the one described in Ref. [22], that
performed a DSNB spectral analysis over the 22.5 ×
2853 kton·days of data corresponding to the first three
SK phases. While this previous analysis did not include
neutron tagging and used slightly different backgroundcategories and reduction cuts, it is also dominated by29
statistical uncertainties. Hence, the SK-I,II,III results
can be readily combined with the SK-IV observations.
In order to study a wide variety of the discrete andparameterized models introduced in Section I, we build
signal PDFs for all four phases of SK using Monte-Carlosimulation to model detector resolution effects over mostof the energy range. For SK-IV, we renormalize the IBD
Monte-Carlo simulation by the predicted DSNB spectrum for all energies. For SK-I,II,III, we renormalize thesimulations used in Ref. [22] —with events generated in
Ee+ in 10−90 MeV— and use the resolution functions
from Refs. [36–38] to extrapolate them to lower energies.
For the Galais+10 [15] and the Nakazato+15 [6] models,
for which DSNB fluxes are available only up to antineutrino energies of 50 MeV, we assume the predicted flux
to be zero in the 50−81.3 MeV range. This assumption
is not expected to affect the spectral fit, as backgrounds
from νe CC interactions largely dominate at these energies. Finally, for each DSNB model, we fit the data samples obtained in Ref. [22] for the first three SK phases.
Here, although we use the same cuts and search regions
as in Ref. [22], we account for possible residual spallation
by modeling the spectra and associated systematic uncertainty using the same procedure as at SK-IV. We then
combine the four fit results by summing the corresponding likelihoods maximized over the nuisance parameters.

E.

Results

We apply the procedure described in the previous sections to the observed data collected during phases I to IV
of SK for the Horiuchi+09 DSNB model [16]. The associated spectral fits for SK-IV and for SK-I,II,III are shownin Figs. 26 and 27 for the DSNB signal and for the spallation and atmospheric neutrino background categories introduced in Section VII B. Since neutron tagging was notpossible in the first three phases of SK, the corresponding
fits span only three regions, defined by the low, medium,and high Cherenkov angle ranges defined in Table VI. For
SK-IV, the low energy bump in the predicted decay electron spectrum for IBD-like events is due to the energy
dependence of the neutron tagging cut, that strongly depends on the background composition. As expected, thedominant backgrounds in the low and high Cherenkov
angle sidebands are backgrounds with visible muons andpions and NCQE backgrounds, respectively. While mostof the signal region is dominated by backgrounds frominvisible muons and pions, νe CC contributions become
sizable above about 50 MeV. For SK-IV, the backgroundrates are extremely low in the IBD-like region, that will
hence be particularly sensitive to IBD signals in spite of
the low efficiency of the neutron tagging cut.
In spite of introducing spallation backgrounds in the
spectral fits, the results shown in Fig. 27 for SK-I,II,III
are highly similar to the ones shown in Ref. [22]. The
number of spallation events for each phase never exceed
2, within the upper limit obtained in Ref. [22]. Since

the current analysis uses a similar procedure with more
inclusive muon selection criteria and 20−25% lower signal efficiencies, the spallation cut effectiveness at SK-IV
should remain on par with the ones observed at previous SK phases. Yet, an important amount of spallation
backgrounds (∼8 events) is predicted below 19.5 MeV
at SK-IV. The origin of this important spallation background will be further studied in future SK analyses withan enhanced performance by gadolinium. Note that thiseffect was not observed in other SK searches probing asimilar energy range. These searches however used specific cuts on e.g. the event timings or opening angles that
we do not apply in this analysis.
Even accounting for possible residual spallation, thenumber of events observed in 17.5−19.5 MeV is higher
than the best-fit predictions from Fig. 26. While this canbe explained by statistical fluctuation, it led us to investigating possible non-spallation sources for the excess
observed over atmospheric neutrino background predictions in 15.5−19.5 MeV. We first ruled out the possibilityof a strong DSNB signal, as it would be associated witha steeply falling spectrum and hence to an unrealistically high supernova rate. The hypothesis of an isolated
astrophysical event is also strongly disfavored since theobserved events are uniformly distributed over the SK-IV
period, as discussed in Section V G. A radioactive origin
for the observed events is also highly unlikely, given theirhigh energies and their uniform spatial distribution. The
event positions are notably not correlated with the locations of calibration sources. Residual spallation therefore
remains the most plausible explanation for the high number of events observed in 15.5−19.5 MeV at SK-IV.
The likelihoods associated to the SK-I to IV fits for
the Horiuchi+09 model are shown in Fig. 28 as a function of the DSNB rate for the reconstructed equivalent
electron kinetic energies Erec > 15.5 MeV. This figure
also shows the likelihood for the combined analysis, witha total exposure of 22.5×5823 kton·days. The associatedlimits on the DSNB flux are shown in Table VIII. The
expected sensitivities of the SK-I,II,III analysis change
by only ∼3% with the inclusion of the spallation backgrounds while the associated observed upper limits on the
DSNB flux become tighter. The expected sensitivities of
the SK-IV and the combined analyses to the DSNB flux
at 90% C.L. are 2.2 ν̄e cm−2 ·sec−1 and 1.5 ν̄e cm−2 ·sec−1 ,respectively. Hence, the predicted flux for this model,
1.9 ν̄e cm−2 ·sec−1 , is well within the reach of the SKI,II,III,IV analysis. The observed upper limits on the
DSNB flux for SK-IV and for the combined analysis areof 4.6 ν̄e cm−2 ·sec−1 and 2.6 ν̄e cm−2 ·sec−1 , respectively.
The best-fit flux for the combined analysis is now of−2
1.30+0.90·sec−1 . The 1.5σ excess observed over
−0.85 ν̄e cm
the background prediction is higher than the 0.9σ excessfrom the previous combined analysis [22]. While not statistically significant, this result is compatible with a wide
range of DSNB predictions with a flux comparable to theone of the Horiuchi+09 model [16].

30

Sideband (20 < C < 38 )17.5

Number of events after cuts
0 / >1 neutron tags

Signal region (38 < C < 50 )Sideband (78 < C < 90 )

e CC

NCQE
Decay e/15.0

Spallation
DSNB
All backgrounds

12.5
10.0
7.5
5.0
2.5
0.0
20

40

50

60

70

80

20

30

40

50

60

70

80

20

30

40

50

60

70

80

50

60

70

80

20

30

40

50

60

70

80

20

30

40

50

60

70

80

Erec [MeV]

Erec [MeV]

Erec [MeV]

e CC
NCQE
Decay e17.5

Number of events after cuts
1 neutron tag

30

/

15.0

Spallation
DSNB
All backgrounds

12.5
10.0
7.5
5.0
2.5
0.0
20

30

40

Erec [MeV]

Erec [MeV]

Erec [MeV]

FIG. 26. Best-fit signal and background spectra for SK-IV, assuming the DSNB spectrum predicted by the Horiuchi+09
model [16]. The six regions presented here include two signal regions and four sidebands in Table VI.

VIII.
A.

DISCUSSION

Constraints on the DSNB parameter space

Over the two search strategies presented in this paper, the analysis described in Section VI gives a modelindependent differential limit on the ν̄e flux while themodel-dependent spectral analysis uses data from allphases of SK and yields the tightest constraints on DSNB
models. In this section, we present the results of thisanalysis using both discrete DSNB models —most ofthem based on supernova simulations— and simplified
parameterizations.
The 90% C.L. expected and observed upper limits on
the DSNB rate and flux for theoretical models are shownin Fig. 29 for the combined SK-I,II,III,IV analysis, withthe corresponding results being tabulated in Table VIII.
With an exposure of 22.5×5823 kton·days, the sensitivityof the combined analysis at 90% C.L. is the tightest to
date and is on par with predictions from the Ando+03,
Horiuchi+09 (Tν = 6 MeV), Galais+10, and Kresse+21
models [12, 15–17]. One common feature of these optimistic models is that they use the highest cosmic star formation history allowed by observations. Additionally, theKresse+21 predictions [12] are based on neutrino emission models derived from state-of-the-art supernova simulations, and take a wide range of progenitors into account. The ability of this analysis to probe such realistic
models makes the prospects for future experiments particularly promising. Currently, an excess of about 1.5σ
has been observed in the data and the resulting 90% C.L.
limits on the DSNB flux improve on the constraints obtained in Ref. [22] for SK-I,II,III. These results confirm
the exclusion of the Totani+95 model [8] and of the most
optimistic predictions of the Kaplinghat model [9].
Another striking feature observed in Fig. 29 is the stability of the observed and expected limits on the DSNB
flux over all models. This stability results from the limited sensitivity of the current analysis to the DSNB spectral shape. In order to translate flux limits into constraints on astrophysical parameters, we therefore ignore
subtle effects associated with e.g. black-hole formation or
multiple classes of progenitors, and use the simple FermiDirac neutrino emission model described in Section III,as in Ref. [22]. The associated two-dimensional 90% C.L.
limit on the neutrino emission temperature and on the
DSNB rate for Eν > 17.3 MeV is shown in Fig. 30, forthe combined SK-I,II,III,IV analysis. This figure also

31

Sideband (20 < C < 38 )17.5

Number of events after cuts

Signal region (38 < C < 50 )Sideband (78 < C < 90 )

e CC
NCQE
Decay e20.0/

Spallation
DSNB
All backgrounds

15.0
12.5
10.0
7.5
5.0
2.5
0.0
25

50

Erec [MeV]75

Sideband (20 < C < 38 )25

50Erec [MeV]

75

25

50

Erec [MeV]75

Signal region (38 < C < 50 )Sideband (78 < C < 90 )

20

20

e CC

NCQE
Decay e12

Number of events after cuts/

Spallation
DSNB
All backgrounds

10
8
6
4
2
0
20

40

60

Erec [MeV]80

Sideband (20 < C < 38 )40

60

Erec [MeV]80

Signal region (38 < C < 50 )40

60Erec [MeV]80

Sideband (78 < C < 90 )

e CC

8

NCQE
Decay eNumber of events after cuts/

Spallation
DSNB
All backgrounds

6

4

2

0
25

50

Erec [MeV]

75

25

50

Erec [MeV]

75

25

50

Erec [MeV]75

FIG. 27. Best-fit signal and background spectra for SK-I (top), II (middle), and III (bottom), assuming the DSNB spectrum
predicted by the Horiuchi+09 model [16]. The three regions presented here correspond to different range of the Cherenkov
angle, as shown in Table VI. Since neutron tagging is not possible in SK-I,II,III, only these three Cherenkov angle regions are
considered.

32

SK­I
SK­II
SK­III
SK­IV
Combined

0.00
0.25

Likelihood

0.50
0.75
1.00
1.25
1.50
1.75
2.00

0.0

2.5

5.0

7.5

10.0

12.5

DSNB [events/year]

15.0

17.5

20.0

FIG. 28. Likelihoods associated with phases I to IV of SK, as well as the combined likelihood. As mentioned in the main text,
we can neglect correlations between background systematic uncertainties when combining all phases of SK.

FIG. 29. The 90% C.L. upper limits, best-fit values, and expected sensitivities for the DSNB fluxes (left) and rates (right)
associated with the models described in Section I. The best fit rates and fluxes are shown with their associated 1σ error bars.
This figure also shows the predictions for each models, either as a range or as one value. Note the stability of the expected andobserved flux limits across all models.

33

9event /22.5kton yr > 16 MeV

3.5

90% C.L. upper limits for individual T s
2D Blackbody DSNB exclusion

SK 1497+794+562+2970 Days
Excluded (E > 17.3 MeV)
90% C.L.

3.0
2.5
2.0
1.5

IMB

1.0

Kamioka

0.5
0.0
2

3

4

5

T [MeV]

6

7

8

FIG. 31. The 90% C.L. excluded region for the total ν̄e energy
emitted during a supernova and the effective neutrino temperature, for the blackbody models described in Section III (redsurface). To guide the eye, we also show the one-dimensional
exclusion limit for individual neutrino temperatures. We also
show the 90% C.L. contours corresponding to the range ofparameters allowed for SN1987A.

90% C.L. upper limits for individual T s
2D Blackbody DSNB exclusion
Predicted rates

8

4.0 1e53

SN e energy in 1053erg

shows the expected rates for different effective neutrino
emission temperatures. Here, the band accounts for thecurrent uncertainty on the cosmic star formation rate,using the low, fiducial, and high estimates from Ref. [16].
Since we used two-dimensional exclusion contours insteadof computing limits for individual models, the limits appear weaker than the ones shown in Fig. 29. These limits
are on par with the SK-I,II,III analysis [22]. In addition
to the rate limits, whose interpretation is complicatedby the high SK analysis threshold, we also show 90%
C.L. exclusion contours for the total ν̄e energy and theeffective neutrino temperature in Fig. 31. This figure
also shows the regions of parameter space allowed by the
Kamiokande II and IMB experiments for the supernova
SN1987A [93]. While the neutrino temperature for thissupernova is particularly low, future experiments shouldbe able to probe at least the associated IMB region.

7
6
5
4
3
2
1
0
2

3

4

5

T [MeV]

6

7

8

FIG. 30. The 90% C.L. excluded region for the DSNB rate forEν > 17.3 MeV and the effective neutrino temperature, forthe blackbody models described in Section III (red surface).
To guide the eye, we also show the one-dimensional exclusion
limit for individual neutrino temperatures. The blue bandshows the predictions for blackbody models, for the low, fiducial, and high star formation rates used in Ref. [16].

B.

Prospects at future experiments

Recent doping of SK with gadolinium (SK-Gd) as wellas the advent of Hyper-Kamiokande (HK) in the moredistant future open new promising perspectives for the
DSNB search. Using gadolinium will notably allow to
considerably improve neutron tagging and hence eliminate most spallation backgrounds and lower the threshold of the analysis. As shown in Fig. 24, however, multiple backgrounds involving neutrons can still dominate
over the DSNB up to about 15.5 MeV. Reactor antineu-

trinos, in particular, provide an irreducible background
up to 9.5 MeV. Up to about 11.5 MeV, 9 Li decays mightthen dominate, and spallation reduction will remain a key
component of future analyses. Locating muon-induced
showers using neutrons will become especially powerful at
SK-Gd. After eliminating spallation, atmospheric NCQE
interactions will dominate over a wide energy range andwill call for new dedicated reduction techniques. At SKGd, notably, the high visibility of the neutron capture
signal will allow to locate neutron capture without using
the positron vertex, which could allow to use the neutron travel distance as a discriminating observable. Finally, while HK may not be doped with Gd, this considerably larger detector coupled with the accelerator neutrino beams will allow to improve the characterization of
the NCQE rate, spectral shape, and neutron multiplicity, that are currently associated with large systematic
uncertainties [48]. In the future, further studies using
the upcoming Intermediate Water Cherenkov Detector(IWCD) will allow to considerably reduce neutron multiplicity uncertainties by studying monochromatic beams,
and will hence be a key piece of the DSNB search program at HK. The DSNB study will also be extended bysearches at lower energies, achieved at future large-scaleliquid scintillation detectors, such as JUNO [94], giving
deeper insights into spectrum shape over the wide energy
range in combination with HK [95].

34

TABLE VIII. Best-fit values and the 90% C.L. upper limits on the DSNB fluxes (in cm−2 ·sec−1 ) for the theoretical models
for phases SK-I to IV as well as for the combined analysis. Here the upper limits are given for Eν > 17.3 MeV. For the
Kresse+21 models, the “High”, “Fid”, and “Low” predictions correspond to the “W20-BH2.7-α2.0”, “W18-BH2.7-α2.0”, and“W18-BH2.7-α2.0-He33” models from Ref. [12] respectively.
Best fit
90% CL limit
Pred.
SK4
All SK1 SK2 SK3 SK4 All+0.9
Totani+95 Constant
2.5+1.4
2.3 6.3 7.0 4.5 2.6
4.67−1.3 1.3−0.9
+0.9
+1.5
Kaplinghat+00 HMA (max) 2.6−1.3 1.3−0.9
2.3 6.7 7.1 4.7 2.6
3.00
+0.9
Horiuchi+09 6 MeV, max
2.6+1.4
2.4 6.0 7.0 4.6 2.6
1.94−1.3 1.3−0.9+0.9
1.4
Ando+03 (updated 05)
2.7+1.5
2.3
6.6
7.2
4.7
2.7
1.74
−0.9
−1.4
+0.9
Kresse+21 (High, NO)
2.7+1.5
2.3 6.7 7.2 4.7 2.7
1.57−1.3 1.4−0.9+0.9
Galais+09 (NO)
2.5+1.4
2.3 6.3 7.0 4.5 2.6
1.56−1.3 1.3−0.9+0.9
1.3
Galais+09 (IO)
2.6+1.4
2.3
6.4
7.0
4.5
2.6
1.50
−0.9
−1.3
+0.9
2.4 6.1 7.1 4.6 2.7
1.23
Horiuchi+18 ξ2.5 = 0.1
2.6+1.4−1.3 1.4−0.9+0.9
Kresse+21 (High, IO)
2.7+1.5
2.3 6.7 7.1 4.7 2.7
1.21−1.3 1.4−0.9+0.9
1.4
Kresse+21 (Fid, NO)
2.7+1.5
2.3
6.8
7.2
4.7
2.7
1.20
−0.9
−1.3
+0.9
Kresse+21 (Fid, IO)
2.7+1.5
2.3 6.8 7.2 4.7 2.7
1.02−1.3 1.4−0.9+0.9
Kresse+21 (Low, NO)
2.7+1.5
2.3 6.8 7.2 4.8 2.7
0.96−1.4 1.4−0.9+0.9
Tabrizi+21 (NO)
2.7+1.5
2.4 6.6 7.1 4.7 2.7
0.92−1.3 1.4−0.9+0.9
Kresse+21 (Low, IO)
2.7+1.5
1.4
2.3
6.8
7.2
4.8
2.7
0.84
−1.4
−0.9
+0.9
Lunardini09 Failed SN
2.8+1.5
2.4 6.8 7.3 4.8 2.8
0.73−1.4 1.4−0.9+0.9
Hartmann+97 CE
2.6+1.4
2.3 6.5 7.1 4.6 2.6
0.63−1.3 1.3−0.9+1.0
Nakazato+15 (max, IO)
2.7+1.5
1.4
2.4
6.5
7.2
4.8
2.7
0.53
−1.4
−0.9
+0.9
Horiuchi+18 ξ2.5 = 0.5
2.7+1.5
2.2 7.1 7.1 4.8 2.6
0.55−1.4 1.3−0.9+0.9
Horiuchi+21
2.1+1.3
3.4 4.3 5.9 3.9 2.5
0.28−1.2 1.2−0.9+0.9
Malaney97 CGI
2.7+1.5
1.3
2.3
6.8
7.1
4.7
2.6
0.26
−1.3
−0.9
+1.0
2.3 6.8 7.2 4.8 2.7
0.19
Nakazato+15 (min, NO)
2.8+1.5−1.4 1.4−0.9
Model

TABLE IX. Best-fit values and the 90% C.L. upper limits on the DSNB rates (in events·year−1 ) for the theoretical models
for phases SK-I to IV as well as for the combined analysis. Here the upper limits are given for Eν > 17.3 MeV. For the
Kresse+21 models, the “High”, “Fid”, and “Low” predictions correspond to the “W20-BH2.7-α2.0”, “W18-BH2.7-α2.0”, and“W18-BH2.7-α2.0-He33” models from Ref. [12] respectively.
Best fit
90% CL limit
Pred.
SK4
All SK1 SK2 SK3 SK4 All+1.6
Totani+95 Constant
4.5+2.5
4.2 11.2 12.2 7.9 4.6
8.35−2.2 2.3−1.5
+2.4
+1.5
Kaplinghat+00 HMA (max) 4.4−2.2 2.2−1.4
3.9 11.2 11.7 7.7 4.4
5.09
+1.7
Horiuchi+09 6 MeV, max
4.6+2.5
4.4 11.7 12.5 8.2 4.8
3.54−2.3 2.4−1.6+1.6
Ando+03 (updated 05)
4.7+2.5
4.2 11.8 12.4 8.2 4.7
3.09−2.4 2.4−1.6+1.6
Kresse+21 (High, NO)
4.5+2.5
4.1
11.7
12.2
8.1
4.6
2.76
2.4
−2.3
−1.5
+1.6
4.0 11.0 11.9 7.8 4.5
2.74
Galais+09 (NO)
4.4+2.4−2.2 2.2−1.4+1.6
4.0 11.0 11.9 7.7 4.5
2.62
Galais+09 (IO)
4.4+2.4−2.2 2.2−1.4+1.7
Horiuchi+18 ξ2.5 = 0.1
4.8+2.6
4.5 12.1 12.8 8.4 4.9
2.29−2.4 2.5−1.6+1.6
Kresse+21 (High, IO)
4.6+2.5
2.4
4.1
11.7
12.2
8.1
4.6
2.14
−2.3
−1.5
+1.5
Kresse+21 (Fid, NO)
4.5+2.5
4.0 11.7 11.9 7.9 4.5
2.06−2.2 2.3−1.5+1.5
Kresse+21 (Fid, IO)
4.5+2.5
4.0 11.6 11.9 7.9 4.5
1.75−2.2 2.3−1.5+1.6
Kresse+21 (Low, NO)
4.5+2.5
2.2
3.9
11.6
11.8
7.9
4.5
1.65
−2.2
−1.4
+1.6
Tabrizi+21 (NO)
4.6+2.5
4.2 11.8 12.3 8.2 4.7
1.64−2.3 2.4−1.5+1.6
Kresse+21 (Low, IO)
4.5+2.5
4.0 11.6 11.9 7.9 4.5
1.43−2.2 2.2−1.4+1.7
Lunardini09 Failed SN
5.0+2.7
2.6
4.5
12.7
13.1
8.8
5.1
1.36
−2.5
−1.7
+1.6
Hartmann+97 CE
4.4+2.4
4.0 11.2 11.9 7.7 4.5
1.09−2.2 2.2−1.4+1.8
Nakazato+15 (max, IO)
5.0+2.7
4.5 12.1 13.0 8.7 5.0
0.99−2.5 2.5−1.6+1.5
Horiuchi+18 ξ2.5 = 0.5
4.2+2.3
3.5 11.4 10.9 7.4 4.1
0.87−2.1 2.1−1.4+1.6
Horiuchi+21
3.7+2.2
2.0
6.0
7.5
10.2
6.8
4.2
0.49
−2.1
−1.6
+1.5
Malaney97 CGI
4.3+2.4
3.8 11.2 11.3 7.6 4.3
0.44−2.1 2.2−1.4+1.6
Nakazato+15 (min, NO)
4.7+2.5
4.0 11.8 12.1 8.2 4.6
0.34−2.4 2.4−1.6
Model

35
IX.

Appendix A: Spallation cut

CONCLUSION

In this paper, we presented two analyses taking advantage of the longest data taking phase of SK, the improved data reduction, and better background estimates.
We developed improved spallation reduction and neutron
tagging algorithms with respect to the previous analyses
of SK-IV data, and constrained the DSNB flux in the
7.5−29.5 MeV and 15.5−79.5 MeV energy ranges. With
the model-independent analysis, we achieved the world’s
tightest upper limit on the extraterrestrial electron antineutrino flux. For the latter region, we notably extended the spectral analysis developed for SK-I,II,III [22],and combined the results of the first four SK phases toreach an exposure of 22.5 × 5823 kton·days. This exposure yields SK a world-leading sensitivity to the DSNB
flux at 90% C.L., comparable to the fluxes of several realistic models. A 1.5σ excess has been observed across all
the models probed and the final 90% C.L. limits on the
DSNB flux are on par with the ones obtained in previous
SK analyses [22, 71]. The last 20 years of data taking at
SK have hence brought us to the doorstep of the DSNB
parameter space. The enhanced neutron tagging capabilities at SK-Gd and the larger detector volume achievedat HK will allow to explore this parameter space, setting
meaningful constraints on astrophysical observables, withthe realistic prospect of a groundbreaking discovery.

1.

For each cloud, we define a specific coordinate system
whose origin is the projection of the cloud barycenteron the muon track. For a given DSNB candidate reconstructed vertex, `l is defined as the distance to this origin
along the muon track while `t is the transverse distanceto the track. For a given cloud multiplicity, a DSNB candidate observed less than dt after the associated muon,
and within the ellipse defined by `l and `t , is classified as
spallation background. Conditions for the neutron cloud
cut are summarized in Table X for different multiplicitiesof the cloud.
TABLE X. Conditions for the neutron cloud cut as a function of the cloud multiplicity. A number followed by the ‘+’
symbol represents a lower bound on this multiplicity.
Multiplicity 2+ 2+ 2
3 4−5 6−9 10+
dt [sec]
0.1 1 30 60 60 60 60`l [cm]1200 800 383 548 603 712 766`t [cm]1200 800 219 268 379 490 548

2.

ACKNOWLEDGMENT

We gratefully acknowledge the cooperation of the
Kamioka Mining and Smelting Company. The SuperKamiokande experiment has been built and operated
from funding by the Japanese Ministry of Education,
Culture, Sports, Science and Technology, the U.S. Department of Energy, and the U.S. National Science Foundation. Some of us have been supported by funds fromthe National Research Foundation of Korea NRF-20090083526 (KNRC) funded by the Ministry of Science,
ICT, and Future Planning and the Ministry of Education (2018R1D1A3B07050696, 2018R1D1A1B07049158),
the Japan Society for the Promotion of Science, the National Natural Science Foundation of China under Grants
No. 11620101004, the Spanish Ministry of Science, Universities and Innovation (grant PGC2018-099388-B-I00),
the Natural Sciences and Engineering Research Council(NSERC) of Canada, the Scinet and Westgrid consortia of Compute Canada, the National Science Centre,
Poland (2015/18/E/ST2/00758), the Science and Technology Facilities Council (STFC) and GridPPP, UK, theEuropean Unions Horizon 2020 Research and Innovation Programmeunder the Marie Sklodowska-Curie grant
agreement no. 754496, H2020-MSCA-RISE-2018 JENNIFER2 grant agreement no.822070, and H2020-MSCARISE-2019 SK2HK grant agreement no. 872549.

Neutron cloud cut

Rectangular cuts

As mentioned in Section V B, we applied a series ofrectangular cuts based on spallation variables. The underlying idea is to remove isotopes with their half-lives
dependent on the energy region and muon type. In particular, for events with a number of tagged neutrons different from one in the spectral analysis, we apply specific rectangular cuts to entirely O(0.01) sec half-lives.
The criteria are summarized in Tables XI and XII for
events with and without exactly one tagged neutron respectively.

3.

Spallation PDFs and log-likelihood ratio

When producing the spallation and random PDFs, we
separate the dt and `t bins to account for the possible correlations, 0−0.05, 0.05−0.5, and 0.5−30 sec fordt (determined by the isotope half-lives in Fig. 3), and
0−300, 300−1000, 1000−5000 cm for `t . Fig. 32 shows
PDFs for each spallation variable from single throughgoing muons. Here the results from the most spallationrich region (0 < dt < 0.05 sec and 0 < `t < 300 cm) are
shown. We calculate the spallation log-likelihood ratiowith the procedure described in Section V B. An example distribution is shown for single through-going muonsin Fig. 33. The cut criteria are tuned depending on thereconstructed energy as well as the dt and `t bins. We
use the common binning for dt, 0−0.05, 0.05−0.5, and36
TABLE XI. Summary on the rectangular spallation cuts for events with exactly one tagged neutron. When a given cut has
requirements on both dt and `t , at least one should be fulfilled.
Muon type

Energy range
dt requirement `t requirement
Erec < 23.5 MeV
dt > 0.001 sec
All
Erec < 23.5 MeV
dt > 0.1 sec`t > 400 cm
Misfit15.5 < Erec < 19.5 MeV dt > 1.5 sec
Well-fitted single through-going 15.5 < Erec < 17.5 MeV
dt > 7 sec`t > 150 cm
Stopping
15.5 < Erec < 19.5 MeV dt > 0.05 sec
Poorly-fitted stopping
15.5 < Erec < 17.5 MeV
dt > 6 sec
All multiple
15.5 < Erec < 19.5 MeV dt > 0.05 sec-

TABLE XII. Summary on the rectangular spallation cuts for events with tagged neutrons different from one. When a given
cut has requirements on both dt and `t , at least one should be fulfilled. Here tighter rectangular cuts are applied to themuon categories generating most of the spallation —single through-going and multiple muons— to remove contributions fromshort-lived isotopes.
Muon type

Energy range
dt requirement `t requirement
Erec < 23.5 MeV
dt > 0.001 sec
All
Erec < 23.5 MeV
dt > 0.1 sec`t > 400 cm
Misfit15.5 < Erec < 19.5 MeV dt > 1.5 sec
15.5 < Erec < 17.5 MeV
dt > 7 sec`t > 150 cm
Well-fitted single through-going
17.5 < Erec < 19.5 MeV
dt > 7 sec`t > 100 cm
Stopping
15.5 < Erec < 19.5 MeV dt > 0.05 sec
All multiple
15.5 < Erec < 19.5 MeV dt > 0.05 sec
15.5 < Erec < 19.5 MeV dt > 0.1 sec
Well-fitted multiple
15.5 < Erec < 19.5 MeV
dt > 1 sec`t > 400 cm

0.5−30 sec, independent of the muon type. In contrast,for `t , we separate bins only for most impactful muontypes, that is, single through-going and multiple muons:0−200, 200−300, 300−500, 500−1000, and >1000 cm
for single through-going, and 0−100, 100−200, 200−300,
300−500, 500−700, 700−1000, 1000−2000, 2000−3000,>3000 cm for multiple muons.
4.

Spallation cut performance estimation

In this section we present our methodology to estimate
the impact of spallation cuts on non-spallation “signal”
events (DSNB, solar, reactor, and atmospheric neutrino
interactions), as well as on spallation backgrounds. In
what follows, we treat 9 Li decays separately as they provide the dominant β+n signature associated with spallation.
As explained in Section V B, random event efficiencies
for multiple spallation and neutron cloud cuts are estimated to be 98% and 99.8%, respectively. We estimate
the signal efficiencies for the rectangular and likelihood
cuts described above by evaluating their impact on the
post-sample, as described in the previous section. The
associated efficiency (random ) can thus be expressed as:
random = Npost,after /Nbefore ,(A1)where Nbefore and Npost,after are the numbers of eventsin the post-sample before and after spallation cuts. We
validate this efficiency by evaluating the impact of spallation cuts on events from solar neutrino elastic interactions, using the opening angle θsun defined in Section IV
as a discriminator.
To estimate the background rejection power of thespallation cuts, we use a sample of DSNB candidate
events verifying cos θsun < 0. In order to suppress remaining radioactivity and atmospheric neutrino events,
we also apply cuts on the effective distance to the walland the Cherenkov angle, that will be defined in Section V C, and require Erec < 15.5 MeV. The spallation
remaining rate (spall ) is then estimated as the fraction of
events remaining after spallation cuts, after the expected
number of remaining atmospheric neutrino events fromMonte-Carlo simulation has been subtracted.
The procedure outlined above allows to accurately estimate the spallation remaining rate for Erec < 15.5 MeV—where spallation largely dominates even after cuts, but
not to assess the impact of spallation cuts on individual isotopes. Using the spallation-dominated sample described above, we therefore refine our approach as follows. To compute the remaining fraction of a specificisotope, we attribute a random dt, drawn from the isotope decay time distribution, to each pair in the sample,
leaving all other observables untouched. This is equivalent to renormalizing the spallation likelihood from Equation 3 by the isotope’s exponential decay time distribution. Then, we apply the likelihood cuts described aboveprob.

prob.

37

0.1

0.4

0.35

Spallation PDF

0.3

Random PDF

0.08

0.25
0.06
0.2
0.15

0.04

0.1
0.02
0.05
0
10−2

10−1

1

0
0

10
|dt| [sec]

prob.

prob.

0.16

500

1000 1500 2000 2500 3000 3500 4000 4500 5000
lt [cm]

−1

10

0.14
0.12

10−2

0.1
0.08

−3

0.06

10

0.04
10−4

0.02

3

0

1000 2000 3000 4000 5000
ll [cm]

0

500

1000

1500

2000

×10
2500
Qµ [p.e.]prob.

0−5000 −4000 −3000 −2000 −1000

10−1

−2

10

−3

10

10−4
−500

3

0

500

×10
1000 1500 2000 2500 3000 3500 4000 4500
Qres [p.e.]

FIG. 32. Spallation and random PDFs for dt (top left), `t (top right), `l (middle left), Qµ (middle right), and Qres (bottom),for the pre (blue) and post (red) samples. Results from the single through-going muons for the dt region of 0−0.05 sec and the`t region of 0−300 cm are shown here.

to the pairs with the modified dt, and compute the fraction of remaining events, pre,iso . For each event, thisfraction can be expressed as:(pair)

pre,iso = iso,likeli × random,iso(A2)

be computed using a post-sample analogous to the one
described in Section V B 2 and, since we reassigned dt for
all pairs in each event, will be artificially lower than thesignal efficiency random . The probability for a given isotope decay to pass the likelihood cuts (iso,likeli ) is thus
given by:(pair)

where iso,likeli is the probability for a pair formed by aspallation event and its parent muon to pass the likelihood cuts. Conversely, random,iso is the probability for
all other (uncorrelated) pairs between the event and preceding muons to pass these cuts. This probability caniso,likeli = pre,iso /random,iso × random .

(A3)

The fraction of remaining events after both preselection
and likelihood cuts (iso,spall ) can then be estimated as:

Events

38

underflow bin

104

overflow bin

Pre sample
Post sample

2-MeV bins from 7.5 to 15.5 MeV, that will have different
isotope compositions. For a given spallation cut, we thencompute the spallation remaining rates for the four subsamples and take the largest deviation from the averagevalue as our systematic uncertainty. This uncertainty depends on the isotope considered but typically lies around
50%.

3

10

Appendix B: 241 Am/Be calibration study forneutron tagging

102

-30

-20

-10

0

10

20

30
LspallFIG. 33. Spallation log-likelihood ratio of the single throughgoing muon for 7.5 < Erec < 9.5 MeV. The values smaller orlarger than −30 or 30 are filled in the underflow or overflow
bins, respectively.

iso,spall = pre,iso /random,iso × random × presel ,spallwhere presel =.
likeli(A4)(A5)

Here, spall is the spallation remaining rate for all isotopes
after spallation cuts and likeli is the rate after likelihood
cuts only. This ratio yields presel ∼ 70% while the preselection cut alone removes about half of the background,highlighting the significant overlap between preselection
and likelihood cuts.
This procedure is used throughout this study to estimate the remaining fractions of various isotopes. For
Erec <15.5 MeV, we will notably use it to estimate theimpact of spallation cuts on the 9 Li background. Above
15.5 MeV, we will estimate the remaining fractions of 8 Li
and 8 B, two isotopes that are particularly difficult to remove due to their long half-lives. This estimate will allow
us to model the spallation background spectrum for the
spectral analysis.
In order to estimate the remaining amounts of individual isotopes without relying on a simulation we made two
major assumptions. First, that the time difference between spallation events and their parent muon is not correlated with other spallation observables. Second, thatthe impact of spallation cuts above 15.5 MeV can be estimated using a sample of lower energy events. Bothassumptions would be verified if the `t , `l , Qµ , and Qres
spallation observables described in Section V B 2 weresimilarly distributed for all isotopes. This is however
not the case as their values will depend on the isotope
production mode. To estimate the effect of this isotope
dependence, we split the sample described above into four

As described in Section V D 3, in order to validate theefficiency of the neutron tagging procedure on real data,
we applied neutron tagging to data collected in the presence an 241 Am/Be source, which, surrounded by a BGO
scintillator, produces a prompt signal and a delayed neutron capture signal as part of its radioactive decay chain.
However, the data collected with the 241 Am/Be source
is not fully modeled by the simulation software; notably,
no modeling of the BGO scintillator is available. The
Monte-Carlo simulation that we use to model 241 Am/Be
events is hence similar to the IBD simulation. Becauseof this limitation, a direct comparison between the neutron tagging efficiency in simulation and calibration data
would be overly pessimistic, as scintillation in the BGO
can produce prompt events that are not accompanied by
the production of a neutron, thus artificially reducing
the observed preselection efficiency of the neutron candidates (see Section V D 1). To account for this, we treateduncertainty associated with the preselection separatelyfrom that of the BDT selection, following the proceduredescribed in Ref. [76].
The preselection efficiency is based on the number ofhits in a given time window and is hence particularly
sensitive to the PMT photon detection efficiency (quantum efficiency, QE). Hence, we determine the associatedsystematic uncertainty by determining the value of QE
for which the predicted N10 distribution fits the observations best. The N10 predicted and observed distributions
are shown in Fig. 34 for the best-fit QE, that is found tobe 2.1% larger than the QE used in the IBD simulationfrom Section III B. The preselection efficiency associatedto this best-fit QE is 3.0% larger than the efficiency predicted by our simulation. Moreover, the 1σ uncertainty
on this best-fit value, found using a χ2 fit, leads to a
2.1% uncertainty on the value of the preselection efficiency. The combined preselection uncertainty is henceof 3.7%.
After the uncertainty on the preselection efficiency IS ,we evaluate the uncertainty associated with the BDT selection cuts by computing the relative efficiencies rel =BDT+IS /IS for the data and the Monte-Carlo simulation for different values of the BDT discriminant. We
determine these efficiencies by fitting the timing distributions of neutron candidates by an exponential plus a
constant. Results for different BDT cuts are shown in
Table XIII for the different 241 Am/Be datasets consid-

39ered here. The maximum discrepancy between data andMonte-Carlo lies around 12% irrespective of the BDT
cut.

4000
3000
2000
1000
0
2.0
1.5
1.0
0.5
0.0

6

8

10

12

14

N10

16

18

2

The neutron multiplicities predicted by the MonteCarlo simulation for atmospheric neutrinos are shownin Fig. 35. Note that a significant number of atmospheric neutrino interactions produce one or more neutrons. However, these neutrons are more energetic thanfor the IBDs of DSNB neutrinos and sometimes travel
further away from their production vertex. Since thecurrent neutron tagging algorithm requires neutrons to
be close to their associated positron vertex, it can also
be used to reduce neutron-producing atmospheric backgrounds. A detailed study of neutrons associated withatmospheric neutrino backgrounds has been performed
in Ref. [76].

22CCQE-likeCC1π
CC other
NCQE
NC1π
NC other

3

10

10

Appendix C: Neutrons produced from atmospheric
neutrino interactions

20

FIG. 34. N10 observed and expected distributions for neutrons. The expected distribution corresponds to the PMT
photon detection efficiency for which the N10 fits the observations best. The observed distribution has been obtainedby subtracting the N10 background distribution from the total observed distribution. This background distribution has
been obtained from data, using random triggers. The lower
predicted number of events for N10 = 6 is due to harsher
continuous dark noise reduction cuts, as mentioned in Section V D 1.

Events

BDT cut Position Year Data MC (Data − MC)/Data
2009 69.4 76.2
−9.8
Center2016 80.0 72.3
9.5
2009 74.6 77.4
−3.7
0.40
Y12
2016 83.8 77.6
7.2
2009 87.6 77.1
12.0
Z152016 80.2 74.1
−7.7
2009 50.0 55.3
−10.6
Center2016 56.1 52.8
5.8
2009 56.4 59.6
−5.7
0.90
Y12
2016 60.8 59.2
2.6
2009 65.8 58.3
11.4
Z15
2016 57.0 55.4
2.8
2009 33.6 35.5
−5.8
Center2016 34.9 33.2
4.8
2009 36.9 37.8
−2.4
0.99
Y12
2016 38.2 37.1
2.8
2009 38.8 34.5
10.5
Z15
2016 31.9 32.2
−0.1

Neutron candidates after preselection

TABLE XIII. Observed (‘Data’) and predicted (‘MC’) relative efficiencies (in %) for different BDT cuts and for the241
Am/Be calibration runs considered in this study. Data
was taken in 2009 and 2016, at three different locations in thetank: in the center of the detector, near the barrel wall (Y12)
and the top wall (Z15).

5000

Data / MC

Combining the 12% BDT reduction uncertainty to the
3.7% preselection uncertainty leads to an overall uncertainty of 12.5% on the neutron tagging efficiency. This
uncertainty also include possible discrepancies in neutron
kinetic energy distributions, like the ones between 9 Li decays and IBD events.

MC neutrons
Data (Am/Be 2016), background subtracted

6000

10

1
−1

10

0

1

2

3

4

5

6

7

8

9
10
Ntagged-n

FIG. 35. Predicted neutron multiplicity for atmospheric neutrino backgrounds. Here we use the atmospheric neutrinoMonte-Carlo simulation in the 7.5−79.5 MeV range, after thenoise reduction cuts.

40
Appendix D: Spectral fitting
1.

Modeling of atmospheric neutrino backgrounds

The probability distribution functions of the atmospheric neutrino background spectra used in the spectral
analysis are displayed in Fig. 36. As mentioned in Section VII, we divide atmospheric neutrino backgrounds
into four categories: µ/π, νe CC, decay electrons, and
NCQE interactions. The decay electron spectrum is obtained from data while the other spectra are obtained
from the atmospheric neutrino Monte-Carlo simulation.
In the extended maximum likelihood fit described in
Section VII, atmospheric background fluxes are treatedas nuisance parameters. In Table XIV, we compare
the predicted numbers of events in the four atmosphericbackground categories to the best-fit values obtained
when considering the Horiuchi+09 model [16]. Since,
as discussed in Section VII, atmospheric events frome.g. pion-producing NC interactions could be categorized
as either NCQE or µ/π, we quote ranges rather thannumbers for these two categories. For decay electrons andνe CC interactions, the predicted and observed numbersof events are within at most 23% of each other.
TABLE XIV. Predicted and best-fit numbers of atmospheric
events using the Horiuchi+09 model [16]. The predicted numbers of events in the NCQE and µ/π categories are quoted as
ranges, since certain types of NC events could belong to either
category after reconstruction.
Predicted
Observed2.

Decay-e νe CC NCQEµ/π
218.5 121.4 81.9−138.4 10.6−67.2
169.5 109.9
88.2
27.9

Modeling of spallation backgrounds

To model the spallation background PDF we vary theisotope composition of the background as described in
Equation 3 and take the nominal spectrum as the average
between the two extreme slopes. We parameterize thisaverage by the following function:(Erec + 0.511 MeV)α
Sspall (Erec ) = N exp −, (D1)
βwhere N is a normalization factor and α, β are determined by a fit and depend on the SK phase considered.
Then, as described in Section VII C 2, we both vary theisotope composition of the spallation spectrum and distort it to account for energy scale and resolution uncertainties. We then fit the resulting extreme spectra by thefunction from Equation D1 times a third order polynomial:(Erec + 0.511 MeV)α
0
Sspall(Erec ) =N 0 exp −β
× [1 + P3 (Erec )] ,(D2)
and thus parameterize the 1σ systematic error distorsions. The nominal spectrum and its associated analytical fit are shown in Fig. 37 for SK-IV with its 1σ uncertainty range. Note that, in order to account for the
steepness of the spallation spectrum while avoiding un0
physical results, we impose Sspall(Erec ) = 0 for energies
larger than the minimal Erec for which the distorted PDF
is zero.

41

38◦ < θC < 50◦

78◦ < θC < 90◦

20◦ < θC < 38◦

38◦ < θC < 50◦

78◦ < θC < 90◦

20◦ < θC < 38◦

38◦ < θC < 50◦

78◦ < θC < 90◦

20◦ < θC < 38◦

38◦ < θC < 50◦

78◦ < θC < 90◦

0 or >1 ntags

0 or >1 ntags

1 ntag1 ntag1 ntag1 ntag0 or >1 ntags

0 or >1 ntags

20◦ < θC < 38◦

FIG. 36. Probability distribution functions used for the fitting of atmospheric backgrounds in the unbinned spectral fit of thedata. As detailed in Section VII, the backgrounds are separated into four categories, and for each category a PDF is definedover three Cherenkov angle regions (two sidebands and one signal region) and two neutron tagging regions (corresponding to=1 or 6=1 tagged neutrons). Each PDF is normalized to 1 over all regions. The blue histograms show the background spectra
obtained from Monte-Carlo simulation or directly from data for the decay electron, while the red lines are the fitted PDFs,
obtained by fitting the binned spectra with a piecewise polynomial. Top left: electrons from µ or π decays. Top right: νe CCinteractions. Bottom left: NCQE interactions. Bottom right: µ/π-producing backgrounds.

42

FIG. 37. Nominal spallation spectrum (solid red) and its
associated analytical fit (dashed black) for SK-IV. The 1σ
systematic uncertainty is shown as a red band.

43[1] M. Fukugita et al., Springer ISBN 978-4-431-67029-2(1994).
[2] K. Kotake et al., Reports on Progress in Physics 69, 971(2006).
[3] K. Nakazato et al., The Astrophysical Journal Supplement Series 205, 2 (2013).
[4] G. Tammann et al., The Astrophysical Journal Supplement 92, 487 (1994).
[5] R. Diehl et al., Nature 439, 45 (2006).
[6] K. Nakazato et al., The Astrophysical Journal 804, 75(2015).
[7] R. A. Malaney, Astroparticle Physics 7, 125 (1997).
[8] T. Totani and K. Sato, Astroparticle Physics 3, 367(1995).
[9] M. Kaplinghat et al., Physical Review D 62, 043001(2000).
[10] C. Lunardini, Physical Review Letters 102, 231101(2009).
[11] S. Horiuchi et al., Monthly Notices of the Royal Astronomical Society 475 (2018).
[12] D. Kresse et al., The Astrophysical Journal 909, 2 (2021).
[13] S. Horiuchi et al., Physical Review D 103, 043003 (2021).
[14] Z. Tabrizi and S. Horiuchi, Journal of Cosmology andAstroparticle Physics 2021, 05 (2021).
[15] S. Galais et al., Physical Review D 81, 053002 (2010).
[16] S. Horiuchi et al., Physical Review D 79, 083013 (2009).
[17] S. Ando et al., Astroparticle Physics 18, 307 (2003).
[18] D. H. Hartmann and S. E. Woosley, Astroparticle Physics
7, 137 (1997).
[19] Next Generation of Nucleon Decay and Neutrino Detectors: Proceedings, Aussois, Savoie, France, 7-9 April
2005 (2005).
[20] D. Cokinos and E. Melkonian, Physical Review C 15,1636 (1977).
[21] M. Malek et al. (Super-Kamiokande Collaboration),
Physical Review Letters 90, 061101 (2003).
[22] K. Bays et al. (Super-Kamiokande Collaboration), Physical Review D 85, 052007 (2012).
[23] H. Zhang et al. (Super-Kamiokande Collaboration), Astroparticle Physics 60, 41 (2015).
[24] A. Gando et al. (KamLAND Collaboration), The Astrophysical Journal 745 (2012).
[25] S.
Abeetal.
(KamLAND
Collaboration),
arXiv:2108.08527 [astro-ph.HE] (2021).
[26] B. Aharmim et al. (SNO Collaboration), The Astrophysical Journal 653 (2006).
[27] B. Aharmim et al. (SNO Collaboration), Physical Review
D 102, 062006 (2020).
[28] M. Agostini et al. (Borexino Collaboration), Astroparticle Physics 125, 102509 (2021).
[29] S. Fukuda et al. (Super-Kamiokande Collaboration), Nuclear Instruments and Methods in Physics Research, Section A 501, 418 (2003).
[30] H. Kume et al., Nuclear Instruments and Methods inPhysics Research, Section A 205 (1983).
[31] K. Abe et al. (Super-Kamiokande Collaboration), Nuclear Instruments and Methods in Physics Research, Section A 737, 253 (2014).
[32] M. Nakahata et al. (Super-Kamiokande Collaboration),
Nuclear Instruments and Methods in Physics Research,
Section A 421, 113 (1999).

[33] E. Blaufuss et al. (Super-Kamiokande Collaboration),
Nuclear Instruments and Methods in Physics Research,
Section A 458, 638 (2001).
[34] H. Nishino et al., Nuclear Instruments and Methods inPhysics Research, Section A 610, 710 (2009).
[35] H. Watanabe et al. (Super-Kamiokande Collaboration),
Astroparticle Physics 31, 320 (2009).
[36] J. Hosaka et al. (Super-Kamiokande Collaboration),
Physical Review D 73, 112001 (2006).
[37] J. Cravens et al. (Super-Kamiokande Collaboration),
Physical Review D 78, 032002 (2008).
[38] K. Abe et al. (Super-Kamiokande Collaboration), Physical Review D 83, 052010 (2011).
[39] K. Abe et al. (Super-Kamiokande Collaboration), Physical Review D 94, 052010 (2016).
[40] M. Smy, Proceedings of the 30th International Cosmic
Ray Conference 5, 1279 (2008).
[41] E. E. Salpeter, The Astrophysical Journal 121, 161(1955).
[42] M. T. Keil et al., The Astrophysical Journal 590, 971(2003).
[43] A. Strumia and F. Vissani, Physics Letters B 564 (2003).
[44] R. Brun et al., Report No. CERN-W5013 (1994).
[45] A. Ankowski et al., Physical Review Letters 108, 052505(2012).
[46] K. Abe et al. (T2K Collaboration), Physical Review D
90, 072012 (2014).
[47] L. Wan et al. (Super-Kamiokande Collaboration), Physical Review D 99, 032005 (2019).
[48] K. Abe et al. (T2K Collaboration), Physical Review D
100, 112009 (2019).
[49] M. Honda et al., Physical Review D 75, 043006 (2007).
[50] M. Honda et al., Physical Review D 83, 123001 (2011).
[51] http://www.icrr.u-tokyo.ac.jp/~mhonda/nflx2011/index.html (HKKM2011 Flux).
[52] Y. Hayato, Acta Physica Polonica B 40, 2477 (2009).
[53] K. Abe et al. (Super-Kamiokande Collaboration),
arXiv:2012.03807 [hep-ex] (2020).
[54] M. Jiang et al. (Super-Kamiokande Collaboration),
Progress of Theoretical and Experimental Physics 2019,
053F01 (2019).
[55] G. Battistoni et al., AIP Conference Proceedings 896, 31(2007).
[56] S. W. Li and J. Beacom, Physical Review C 89, 045801(2014).
[57] S. W. Li and J. Beacom, Physical Review D 91, 105005(2015).
[58] S. W. Li and J. Beacom, Physical Review D 92, 105033(2015).
[59] Y. Zhang et al. (Super-Kamiokande Collaboration),
Physical Review D 93, 012004 (2016).
[60] https://github.com/Goldie643/SKReact (SKReact).
[61] M. Baldoncini et al., Physical Review D 91, 065002(2015).
[62] Technical Reports Series No. 428 (International Atomic
Energy Agency, Vienna, 2005).
[63] F. P. An et al. (Daya Bay), Chin. Phys. C 41, 013002(2017), arXiv:1607.05378 [hep-ex].
[64] S. M. Locke et al. (Super-Kamiokande Collaboration), toappear (2021).
[65] G. Carminati, in 32nd International Cosmic Ray Confer-

44
ence (2011).
[66] G. Carminati (Super-Kamiokande collaboration),
Physics Procedia 61, 666 (2015).
[67] M. Elnimr (Super-Kamiokande collaboration), Journal of
Physics Conference Series. 888, 012189 (2017).
[68] Z. Conner, Ph.D Thesis, Boston University (1997).
[69] S. Desai, Ph.D Thesis, Boston University (2004).
[70] K. Abe et al. (T2K Collaboration), Physical Review D
96, 092006 (2017).
[71] H. Zhang, Ph.D Thesis, Tsinghua University (2012).
[72] Y. Zhang, Ph.D Thesis, Tsinghua University (2015).
[73] Y. Nakano et al., Nuclear Instruments and Methods inPhysics Research, Section A 977, 164297 (2020).
[74] T. Chen and C. Guestrin, in Proceedings of the 22nd
ACM SIGKDD International Conference on Knowledge
Discovery and Data Mining, KDD ’16 (ACM, New York,NY, USA, 2016) pp. 785–794.
[75] R. Akutsu, Ph.D Thesis, The University of Tokyo (2019).
[76] K. Abe et al. (Super-Kamiokande Collaboration), to appear (2021).
[77] K. Abe et al. (Super-Kamiokande Collaboration), to appear (2021).
[78] W. A. Rolke et al., Nuclear Instruments and Methods inPhysics Research, Section A 551, 493 (2005).
[79] O. Benhar et al., Nuclear Physics A 579, 493 (1994).

[80] O. Benhar et al., Physical Review D 72, 053005 (2005).
[81] C. J. Horowitz and B. D. Serot, Nuclear Physics A 368,
503 (1981).
[82] C. Maieron et al., Physical Review C 68, 048501 (2003).
[83] J. A. Caballero et al., Physical Review Letters 95, 252502(2005).
[84] R. Gonzalez-Jimenez et al., Physical Review C 88,
025502 (2013).
[85] J. E. Amaro et al., Physical Review C 71, 015501 (2005).
[86] J. E. Amaro et al., Physical Review C 73, 035503 (2006).
[87] F. Capuzzi et al., Nuclear Physics A 524, 681 (1991).
[88] A. Meucci et al., Nuclear Physics A 739, 277 (2004).
[89] A. Meucci and C. Giusti, Physical Review D 89, 057302(2014).
[90] Y. Ashie et al. (Super-Kamiokande Collaboration), Physical Review D 71, 112005 (2005).
[91] C. Regis et al. (Super-Kamiokande), Physical Review D
86, 012006 (2012).
[92] W. T. Winter et al., Physical Review C 73, 025503(2006).
[93] B. Jegerlehner et al., Physical Review D 54, 1194 (1996).
[94] M. Sisti (JUNO Collaboration), Journal of Physics Conference Series 1468, 012150 (2020).
[95] K. Møller, A. M. Suliga, I. Tamborra, and P. B. Denton,
Journal of Cosmology and Astroparticle Physics 05, 066(2018).