PHYSICAL REVIEW D 104, 122002 (2021)

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,1S. 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,3D. Bravo-Berguñ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,9 I. 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,11 O. Drapier,11 S. El Hedri ,11,† A. Giampaolo ,11,‡ M. Gonin,11 Th. A. Mueller,11P. 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,¶ M. Hartz,47 A. Konaka,47,¶ 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, Massachusetts 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, California 92697-4575, USA
8
Department of Physics, California State University, Dominguez Hills, Carson, California 90747, USA
9
Institute for Universe and Elementary Particles, Chonnam National University, Gwangju 61186, Korea
10
Department of Physics, Duke University, Durham, North Carolina 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, Hawaii 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

2470-0010=2021=104(12)=122002(41)122002-1© 2021 American Physical Society

K. ABE et al.

PHYS. REV. D 104, 122002 (2021)

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, United Kingdom
23High Energy Accelerator Research Organization (KEK), Tsukuba, Ibaraki 305-0801, Japan
24
Department of Physics, Kobe University, Kobe, Hyogo 657-8501, Japan
25
Department of Physics, Kyoto University, Kyoto, Kyoto 606-8502, Japan
26
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,
New York 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, United Kingdom
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, and Daresbury 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 of Tokyo 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, 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, United Kingdom
51
Department of Physics, University of Winnipeg, 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
20(Received 24 September 2021; accepted 2 November 2021; published 10 December 2021)
We have conducted a new search for the diffuse supernova neutrino background (DSNB) flux at SuperKamiokande (SK), with a 22.5 × 2970-kton · day exposure from its fourth operational phase IV. With the
new analysis we improve on the existing background reduction techniques and systematic uncertainties andtake advantage of an improved neutron tagging algorithm to lower the energy threshold compared to theprevious 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*

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

122002-2

DIFFUSE SUPERNOVA NEUTRINO BACKGROUND SEARCH AT …

PHYS. REV. D 104, 122002 (2021)

three phases of SK to perform a joint analysis using 22.5 × 5823 kton · days of data. This analysis has theworld’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% CL upper limits on the DSNB ν̄e flux lie
around 2.7 cm−2 · sec−1 , strongly disfavoring the most optimistic predictions. Finally, potentialities of thegadolinium phase of SK and the future Hyper-Kamiokande experiment are discussed.
DOI: 10.1103/PhysRevD.104.122002

I. INTRODUCTION
A. Diffuse supernova neutrino backgroundCore-collapse supernovae (CCSNe) are among the most
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 of the coreof the collapsing star. Information about this core could beaccessed by detecting neutrinos emitted by supernova
bursts, whose luminosity and energy spectra closely trackthe different steps of the CCSN mechanism in a neutrinoheating scenario (see Refs. [1–3], for example). Existingneutrino experiments, however, are mostly sensitive tosupernova bursts occurring in our galaxy and its immediatesurroundings that are extremely rare (a few 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 distantsupernovae. The integrated flux of these neutrinos forms
the diffuse supernova neutrino background, or DSNB(DSNB neutrinos are also referred to as supernova relicneutrinos (SRNs) in various publications).
The DSNB is composed of neutrinos of all flavors whoseenergies have been redshifted when propagating to Earth.
Its spectrum therefore contains unique information not only
on the supernova neutrino emission process but also on thestar formation and Universe expansion history. Spectra
predicted from various models are shown in Fig. 1. Theiroverall 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 is mainly affected by the effective energies
of supernova neutrinos. Other factors shaping the DSNB
spectrum include the neutrino mass ordering, the initial
mass function of progenitors, the equation of state ofneutron stars, the fraction of black-hole-forming supernovae, the revival time of the shock wave, etc. Possiblecombinations of these factors were systematically studied
in the Nakazato þ 15 model [6], the maximal and minimal
fluxes of 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 upper bound on the expected DSNB flux, on par withthe most optimistic predictions of the Kaplinghat þ 00

model [9], also shown in the figure. Between these bounds,
the DSNB 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 corecollapse supernovae lead to black hole formation.
Similarly, the Horiuchi þ 18 model [11] incorporates afraction of black-hole-forming supernovae determinedusing the compactness of the stars. 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 theimpact 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 ofphysical 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 by current
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 these energies,
the dominant detection channel in most experiments is
the inverse beta decay (IBD) of electron antineutrinos(ν̄e þ p → eþ þ n), and contributions from subleading
channels such as electron neutrino elastic scattering or
charged-current interactions on oxygen can be neglected.
This process produces an easily identifiable positron and aneutron (Fig. 2). Here the neutron does not emit Cherenkov
or scintillation light directly, but its capture on a proton leads
to the emission of a 2.2 MeV photon. In pure water, thecharacteristic 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

122002-3

K. ABE et al.

PHYS. REV. D 104, 122002 (2021)

10
Horiuchi+21 (Extrapolated, α λ =0.1, NO)
Tabrizi+21 (NS+BH, NO)
Kresse+21 (W20-BH2.7- α 2.0, NO)
Horiuchi+18 (ξDSNB νe Flux [/cm2/sec/MeV]1= 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
νe Energy [MeV]7080

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 models. In thelegend, “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].

FIG. 2. Schematic illustration of an IBD process and thesubsequent neutron capture on another proton. The characteristic
neutron capture time in water is τ ¼ 204.8  0.4 μs [20].

analysis is sensitive to electron neutrinos and, due to theirreducible solar neutrino background, its effective energy
threshold lies around the hep solar neutrino flux end point,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 atlower energies. At SK the first DSNB search was carriedout 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 alreadyallowed us 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 performed at SK, using a22.5 × 2853-kton · day exposure, a lower neutrino energy

122002-4

DIFFUSE SUPERNOVA NEUTRINO BACKGROUND SEARCH AT …
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 spallationand solar neutrino interactions make it extremely difficult
to search for DSNB antineutrinos based on positron
identification alone. A search for extragalactic antineutrinos performed at KamLAND in 2021 [25] probed neutrinoenergies ranging from 8.3 to 30.8 MeV by investigating
coincident positron and neutron capture signals. Whileneutron identification in a liquid scintillator detector suchas KamLAND is significantly easier than in pure water,
KamLAND’s small fiducial volume 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 us to record the neutron capture
signal until 2008, this analysis could only be performed ona small part of the SK data set, with a total live time of
22.5 × 960 kton · days. Due to this low exposure and the
low neutron tagging efficiency in water, this search,however, yielded weaker limits than the SK-I,II,III analysis
22]] above 17.3 MeV.
In this study, we draw on the previous SK analyses topresent 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–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 a22.5 × 2970-kton · day data set. In the 17.3–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 the ones
obtained in Ref. [22] for the former SK phases, thus
analyzing 22.5 × 5823-kton · days of data. This unprecedented exposure will allow us to probe the DSNB withan unmatched sensitivity.
The rest of this article proceeds as follows. First, the SK
detector and the specific features of its data acquisition
system are described in Sec. II. Then, details about themodeling of the DSNB signal and the different backgrounds are given in Secs. III and IV, respectively. We thendescribe the data reduction process in Sec. V. The procedures associated with the DSNB model-independentand spectral analyses are described in Secs. VI and VII,
respectively. Finally we discuss the current constraints on
the DSNB flux and future opportunities at SK withgadolinium and Hyper-Kamiokande in Sec. VIII before
concluding in Sec. IX.
II. SUPER-KAMIOKANDE
Super-Kamiokande is a 50-kton water Cherenkov detector located in the Kamioka mine, Japan. It is structured by acylindrical stainless steel tank with a diameter of 39.3 m
and height of 41.4 m and consists of two parts: an outer

PHYS. REV. D 104, 122002 (2021)

detector (OD) that serves as a muon veto and an innerdetector (ID) where neutrino detection takes place. In orderto reduce backgrounds due to radioactivity near the detector
wall, most analyses consider only events reconstructed atleast 2 m away from the ID wall, thus defining a 22.5-kton
fiducial volume (FV). It is located 1000 m underground,which allows reduction of the cosmic ray muon flux by afactor of 105 . In order to ensure a high-quality data-taking,conditions inside SK are tightly controlled; water isconstantly recirculated and purified, and the ID wall iscovered with 11,129 20-inch photomultipliers (PMTs) witha 3-ns time resolution, corresponding to a 40% photocathode coverage. These features allow SK to detect
particles with energies ranging from a few MeV to afew TeVs. The OD includes 1,885 eight-inch PMTs, facing
outwards, to detect the Cherenkov light from muons.
Further detailed descriptions of the SK detector and its
calibration can be found in Refs. [29–33].
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. Due toan accident following the maintenance, which resulted in aloss of 60% of the ID PMTs, SK operated with a reduced
photocathode coverage for 794 days (phase II). The coverage was brought back to its nominal value for phase III,
which lasted 562 days. SK’s previous spectral analysis of
the DSNB [22] used data from all three phases. Since 2008,
the front-end electronics has been replaced [34] and a new
trigger system that allows for neutron tagging has been set
up [35]. SK operated with these new electronics duringphase IV, for 2970 days, until being stopped for refurbishment work in May 2018. After the maintenance, SK
operated for a short time with pure water for calibrationand monitoring purposes (phase V) beginning 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 the opportunities
offered by the gadolinium-loaded SK and by the future
Hyper-Kamiokande experiments.
Data processing in SK makes use of multiple triggers,
corresponding to different thresholds on the numberof 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 SK-IV,
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 super-highenergy (“SHE”) trigger, which requires 58 PMT hits
(70 hits before September 2011) in 200 ns. The data in
a ½−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

122002-5

K. ABE et al.

PHYS. REV. D 104, 122002 (2021)

TABLE I. Summary on the SHE and AFT trigger conditions
and live time for the three different periods in SK-IV.
SHE threshold
70 hits
70 hits
58 hits

AFT window

Live time

350 μs500 μs500 μs25.0 days
869.8 days
2075.3 days

each SHE trigger that is not associated with a predetermined OD trigger. The trigger conditions for the different
periods in SK-IV are summarized in Table I. In the rest of
this paper, we will describe the analysis of hit patterns inthese 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 to reconstruct
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 theelectron equivalent kinetic energy Erec. This quantity canalso be interpreted as the total reconstructed positronenergy, and will be used to present the results of this study.
III. SIGNAL MODELING
A. DSNB 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 parametrized
models we compute the DSNB flux using the following
formula:Z X
XcΦðEν Þ ¼
RSN ðz; sÞ Fi ðEν ð1 þ zÞ; sÞH0sνi ;ν̄idz
× pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ;ΩM ð1 þ zÞ3 þ ΩΛ

ð1Þ

where Eν is the neutrino energy, c is the speed of light invacuum, H 0 is the Hubble constant, z is the redshift, RSN ðzÞ
is the redshift-dependent supernova rate, and Fi is thesupernova neutrino emission spectrum for a given flavor i.
The index s represents the different possible classes ofsupernovae—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 the matter and dark energy
contributions to the energy density of the Universe,
respectively. We evaluate the supernova rate using aphenomenological model described in Ref. [16], which
extracts the redshift dependence of the cosmic star formation history by fitting observations and uses the Salpeter
initial mass function [41] to obtain the 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 a model with effective
temperature T ν is thus given by
Fi ðEν Þ ¼ f ν;i Etotν

120 E2ν;i Eν;i =T ν;iðeþ 1Þ−1 ;7π 4 T 4ν;ið2Þwhere Etot
ν is the total energy of the neutrinos emitted by thesupernovae and f ν;i represents the fraction of neutrinos orantineutrinos 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 andluminosities. In what follows, we will refer to blackbody 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 collapsing star,supernova neutrino emission is better modeled using a more
sharply peaked, “pinched,” Fermi-Dirac spectrum [42];however, the current exposure at SK does not allow us toprobe this pinching effect. In this analysis we neglect theeffects of degeneracies between the pinching parameter and,e.g., the neutrino temperature on the final constraints onneutrino emission spectra.
This analysis exclusively focuses on the detection ofelectron antineutrinos via IBD processes. In the rest of thispaper, we model these processes using the Strumia-Vissani
IBD calculation [43], which approximates the IBD crosssection 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 maintain thestability 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 to an
antineutrino, a positron, and a neutron. In a view of thewide diversity of DSNB models, we generate uniformly
distributed positron energies in 1–90 MeV and renormalized events later to model specific spectra. This procedure
will also allow us 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

122002-6

DIFFUSE SUPERNOVA NEUTRINO BACKGROUND SEARCH AT …
reproduces the properties of the water in SK, as well as
models the PMT properties and electronic response.
Detector properties such as water transparency and PMT
noise are being measured daily, allowing us to accurately
model the time-dependent noise contamination of thepositron 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 todevelop 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.
IV. BACKGROUND SOURCES
A. Atmospheric neutrinos
Atmospheric neutrinos are an important background in
the present DSNB searches. Below about 15 MeV, neutralcurrent 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 and pion production
dominate at higher energies. A typical visible signature
for these interactions is the Cherenkov light from electrons
produced by muon and pion decays. Hence, the associatedreconstructed energies will follow a Michel spectrum forreconstructed energies of about 15–50 MeV. Atmospheric
neutrino events are simulated using the HKKM 2011 flux
[49–51] as an input to the neutrino event generator NEUT
5.3.6 [52] to model interactions. Details about the model
parameters used in this analysis can be found in Ref. [53].
The nuclear de-excitation γs are simulated based on the
spectroscopic factors for the p1=2 , p3=2 , and s1=2 states
calculated in Ref. [45]. More detailed descriptions of theexcited states can be found in Refs. [46,48] for the T2K
NCQE measurement. One difference between our procedure and Ref. [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 is fromthe latest atmospheric neutrino analysis in SK [54]. The
systematic uncertainty regarding this treatment is considered in the scaling factor used in Sec. VI A.
The detector simulation for atmospheric neutrino eventsis 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 Sec. VI A and found to be subdominant. Finally, random trigger data are overlaid on top of
neutron capture signals in order to simulate a backgroundfor neutron tagging, following the procedure described inSec. III B. Note that for this last step the time evolution of
the detector properties is taken into account.

PHYS. REV. D 104, 122002 (2021)

B. Cosmic ray muon spallation
In spite of its 2700-m water-equivalent overburden, SK isstill 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 aspectof this analysis. Since most spallation isotope decays only
produce β and γ rays, spallation backgrounds can besignificantly reduced using neutron tagging. Nonetheless,accidental pairing between prompt β or γ events and PMT
hits due to dark noise or intrinsic radioactivity will allow
a sizable number of spallation events to pass as IBDs.
Furthermore, a few isotopes undergo a β þ n decay thatmimics 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, andneutron 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 14B and 11Li. The
spallation cuts used in this analysis will therefore be

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.

122002-7

K. ABE et al.

PHYS. REV. D 104, 122002 (2021)

applied up to 23.5 MeV reconstructed kinetic electronenergy in order to account for energy resolution effects.
Note also that the isotope half-lives range from Oð0.01Þ sec
to a few tens of seconds. Given the high rate of cosmic raymuons, pairing a given isotope decay with its parent process
is a particularly difficult endeavor.
Isotopes that undergo β þ n decays cannot be rejectedusing neutron tagging and therefore need to be modeledwith particular care. This is notably the case for 8He, 11Li,
16
C, and 9Li. Here, 11Li is particularly short-lived andtherefore easy to eliminate using mild spallation cuts. Its
production yield is also expected to be particularly low(around 10−9 μ−1 · g−1 · cm2 [56]). In addition, 8He and 16C
will remain subdominant due to their low production yields(around 0.23 and 0.02 × 10−7 μ−1 · g−1 · cm2 [56]) and end
point β decay energies. On the other hand, 9Li has a nonnegligible yield (around 1.9 × 10−7 μ−1 · g−1 · cm2 [56]), a
medium half-life of about 0.18 sec, and decays into a β þ npair with a branching ratio 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 decay of 9Li is hence similar to the
IBD of a DSNB neutrino, except that the β from 9Li
is an electron and the neutron energy is higher thanthat from IBDs. Since SK does not distinguish electrons
from positrons and does not allow us to measure theneutron energy, we model 9Li decay using the IBD
Monte-Carlo simulation, renormalizing events to reproduce
the 9Li β þ n spectrum. Finally, we estimate the total
expected number of 9Li β þ n decays in the DSNB searchwindow by using an earlier SK measurement [59] that
estimated the 9Li production rate to be 0.86  0.12ðstatÞ
0.15ðsystÞ kton−1 · day−1 . Note that the difference in neutron energies with respect to IBD processes leads to asystematic uncertainty. However, this uncertainty is taken
into account when calibrating neutron tagging efficiencies,

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

as neutrons from the calibration source have similarenergies to the 9Li decay products. More details about thiscalibration procedure are given in Sec. V D.
C. Reactor neutrinos
Antineutrinos emitted from nearby nuclear reactors canundergo IBDs in SK and therefore be mistaken as 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 code uses 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 thereactor antineutrino flux due to changes in reactor activity.
The expected total reactor neutrino flux at SK can thus bepredicted up to Emax
¼ 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 Emax
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 þ 096 MeV model [16]. Reactor antineutrino
backgrounds thus effectively set a lower energy thresholdfor the DSNB search.
D. Solar neutrinos
Electron neutrinos from the Sun cause elastic interactions with electrons in SK. In the DSNB search window,

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

122002-8

DIFFUSE SUPERNOVA NEUTRINO BACKGROUND SEARCH AT …
the 8B and hep solar neutrino fluxes represent an importantbackground up to ∼20 MeV. This background can be
drastically reduced by neutron tagging, although a small
fraction of electrons from solar neutrinos will be accidentally paired with background fluctuations or radioactive
decays, like for the spallation backgrounds. In the absenceof neutron tagging, a solar neutrino interaction can still beidentified by comparing the reconstructed directions of thescattered electron to the direction of the Sun at the detection
time. The angle between these two directions, called theopening angle, is expected to be close 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 forthe latest SK-IV solar neutrino analysis [39] to model solar
neutrino spectra and evaluate 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 fromthe 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 usedin this analysis: noise reduction in Sec. VA, spallation
reduction in Sec. V B, DSNB positron candidate selection
in Sec. V C, and neutron tagging in Sec. 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 (seeTable 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 PMTnoise, calibration events, radioactivity from the surrounding rock and the detector wall, and cosmic ray activity. In
particular, we require candidate events to be located at least
2 m away from the ID wall, thus defining a FVof 22.5 kton.
Additionally, we remove events associated with an OD
trigger, or within 50 μs of another event with reconstructed
electron kinetic energy larger than 6 MeV, in order toremove electrons produced by the decay of low energy
cosmic ray muons stopping in the detector. In addition, we

PHYS. REV. D 104, 122002 (2021)

apply fit quality cuts to eliminate noiselike events. Inside
the FV these noise reduction cuts have a signal efficiency
larger than 99% as confirmed by the IBD Monte-Carlosimulation. Finally, we require all SHE events to befollowed by an AFT to allow for neutron tagging. Due
to occasional deadtime induced by trigger software failure,
this step 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 somelong-lived spallation isotopes such as 16N, with minimum
harm to the signal efficiency.
First, since energetic muons could induce the productionof 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 such as
16N 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” using asample of low energy events from radioactive decays at SK(4 < Erec < 5 MeV), whose reconstructed vertices havebeen replaced by random vertices inside the ID. The cut
parameters introduced above allow us to remove about 45%of the background with a 98% signal efficiency, as can beseen in Fig. 6. Since this cut removes low energy events
clustered in space and time, it also removes low energymuons that are misclassified as electrons with tens of MeV
energies. These muons are not targeted by the otherspallation cuts applied in this study, which assume thatthe prompt event is an electron or a positron. Consequently,
in this study, we will apply this cut over the whole 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 inmuon-induced hadronic showers, which are the birthplacesof 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

122002-9

K. ABE et al.

PHYS. REV. D 104, 122002 (2021)

FIG. 6. Minimal distance between a given DSNB candidate(blue) or SK low energy event with a random vertex (orange) anda 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 this analysis with a red
vertical line.

threshold on the number of PMT hits, applies cuts on theevent quality and distance of the reconstructed vertex to the
ID wall. Events triggered by WIT less than 500 μs after agiven muon are considered as neutron captures if theirreconstructed vertex is within 5 m of the reconstructedmuon track. If multiple neutron captures are observed after
a muon event, we define a neutron cloud using the criteria
from Ref. [64]. We then remove DSNB candidates that arefound close in time and space to these neutron clouds,
following the procedure introduced in Ref. [64]. Moredetails about the cut criteria are given in Appendix A. To
estimate the signal and background efficiencies of theneutron cloud cut, we pair muons with DSNB candidates
in data, defining two samples: a “presample,” including
muons found up to 60 sec before a DSNB candidate, anda “postsample,” with muons found up to 30 sec after acandidate. While the presample contains a mixture ofcorrelated and uncorrelated pairs, the postsample contains
only uncorrelated pairs. This sample can hence be used toboth readily estimate signal efficiencies and characterize
the properties of the pairs formed by each isotope decay
and their parent muon in the presample. An example of how
to extract the distribution of the distances between spallation events and the neutron clouds of their parent muons isshown in Fig. 7. We optimize the neutron cloud cut by
estimating the shape of the neutron clouds from data andmaximizing the statistical significance of the signal over the
spallation background. The performance of this cut is,however, currently limited by the weakness of the neutron
capture signal and the fact that the WIT has only been
available during the last 388 live days of SK-IV. While thesignal efficiency is larger than 99% this cut only removes
about 10% of the spallation background in this analysis.

FIG. 7. Distributions of distance between DSNB candidates
and neutron clouds for pre- and postsamples, where the selected
neutron clouds have been observed up to 60 sec before and after
the DSNB candidates, respectively. The distribution for thepostsample 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 muons appear as a peak in
the presample below 10 m. The associated distribution can beextracted by subtracting the rescaled postsample from thepresample.

Over the 388 days of the WIT live time, however, it
removes about 40% of spallation isotope decays.
Finally, it should be noted that the multiple spallationand the neutron cloud cuts described above are expected tooverlap 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 us 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 thecorrelations between muons and DSNB candidates isnecessary. In particular, in order to identify the decaysof spallation isotopes, it is crucial to associate each isotope
with its parent muon. Due to the long half-lives of isotopes
such as 11Be and 16N, we need to investigate all possible
pairings between DSNB candidates and muons observed
up to 30 sec before them. This window size allows us toaccommodate most decay products, without processing aprohibitively high number of muons. As with neutron cloud
cut optimization, we define pre- and postsamples in order toextract the observable distributions associated with 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.

122002-10

DIFFUSE SUPERNOVA NEUTRINO BACKGROUND SEARCH AT …
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,
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 directionof 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%), andcorner clippers (4.3%), and the fractions of the different
categories in the SK-IV data are shown in parentheses. Misfit
muons with a charge lower than 1000 p:e: are removed inorder to reject nonmuonic high energy events. This cut does
not affect removing the spallation background because thetrack length of these low charge muons through the ID istypically 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.
(i) dt: time difference between a muon and the DSNB
candidate. For spallation events, this observablereflects the half-lives of the produced isotopes.
(ii) lt : transverse distance of the DSNB candidate to themuon track (see Fig. 8). This observable is related tothe path lengths of the secondary particles frommuon showers. For well-fitted single through-going
muons paired with their associated spallation isotope
decays, lt is typically no larger than a few meters.
(iii) ll : longitudinal distance to the point of the muontrack associated with the maximal energy deposition(see Fig. 8). The energy deposition along the muontrack is determined by projecting back the light seen
by each PMT to a point on the muon track, following
muon track

position wheredE/dx is largest

ll

DSNB candidate

lt

FIG. 8. Schematic of the spallation observables lt and ll . The
cylindrical shape represents the SK ID.

PHYS. REV. D 104, 122002 (2021)

the procedure described in Ref. [22]. So ll provides
an estimate of the distance between the low energy
event and the origin of the particle shower and, if theshower position is correctly identified, should also
not be larger than a few meters for spallation pairs.
(iv) Qμ : total charge deposited by muons in the detector.
Higher values of this observable are expected when ashower is produced. While this observable also
includes contributions from minimum ionization
along the muon track—and hence strongly depends
on the muon track length—it provides the mostrobust shower predictor for poorly fitted muonsand muon bundles. Conversely, for well-fitted single
through-going muons, the shower can be moreaccurately characterized by the Qres observable
described below.
(v) Qres : residual charge deposited by a muon comparedto the value expected from the minimum ionization.
This observable can be expressed as
Qres ¼ Qμ − qMI · L;
where qMI and L are the number of p:e: percentimeter expected from the minimum ionizationand the track length, respectively. Here qMI depends
on detector properties, such as water transparency,that vary with time and are therefore recomputed foreach run. For muons with multiple tracks, L is thelength of the first track. This observable allows us todetermine how likely a given muon is to induce ashower.
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 rest of
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 de-excitation γ rays, and are 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 ashort-lived isotope—followed by a 2.2 MeV γ ray fromneutron capture can mimic an IBD signature. To preventthese from leaking into the sample, we require DSNB
candidates to be more than 1 ms away from any muon. The
impact of this cut on the signal efficiency is negligible.
Additionally, we apply a set of cuts on dt and lt fordifferent energy ranges, exploiting the dependence of the

122002-11

K. ABE et al.

PHYS. REV. D 104, 122002 (2021)

isotopes’ half-lives in their end point energies shown in
Fig. 3. Detailed descriptions of these rectangular cuts are
given in Appendix A. These cuts are particularly efficient
above 15.5 MeV where short-lived isotopes dominate. In
15.5–19.5 MeV, notably, they allow us to remove about
85% of the spallation background while keeping 88% of
the signal events. Finally, in order to further eliminate
spallation backgrounds, we use distributions of the different observables considered here to define log-likelihood
ratios. First, we prepare two probability density functions(PDFs): the spallation PDF (PDFispall ) and the random PDF(PDFirandom ), for each variable i ¼ dt; lt ; ll ; Qμ ; Qres . We
isolate contributions from spallation events by subtracting
the postsample distributions from the corresponding
presample distributions and obtain PDFispall after area
normalization of these subtracted distributions. For
PDFirandom we normalize postsample distributions by theirareas. We repeat this procedure for each category of muons.
For single through-going and multiple-track muons,
which generate most of the spallation background, wetune the PDFs in dt and lt bins. These bins, adjusted byconsidering the typical half-lives of the isotopes and the lt
distributions, account for the correlations between theseobservables. For each set of PDFs, we then define loglikelihood ratios as
Lspall ¼ log

YPDFispall:
PDFirandomið3Þ

Separate likelihoods are defined for each muon category,
except for corner clippers, whose contributions to spallation backgrounds are negligible, as mentioned above. Notethat for misfit muons only dt is used to calculate loglikelihood 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 for complex correlations between spallation observables and geometrical and muon reconstruction effects.
The signal efficiencies and background rejection rates of
the resulting cuts are estimated using the random sample
introduced in this section, as well as spallation-dominateddata samples. These samples, as well as our methodology
to estimate the spallation cut performance, are describedin detail in Appendix A. Note that the estimates of thespallation remaining rate presented there are used only
for cut optimization and not for the final background
predictions shown in Secs. VI and VII. The current cuts
achieve significant improvement over the previous cuts
used in Refs. [22,23], especially at low energies; for thesame 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 for higher
energies.

C. DSNB positron candidate selection
In the region above ∼20 MeV, spallation backgrounds canbe reduced to negligible levels using a series of spallation cuts
while keeping most of the signal. However, significant
backgrounds from atmospheric neutrino interactions andradioactive decays remain. To identify them, we define thefollowing discriminating observables, aimed at characterizing
the prompt event.
1. Incoming event cut
Radioactivity near the detector wall, as well as muon
spallation in the rock surrounding the detector, can lead toelectrons or γ rays entering the FV. Instead of tightening theFV cut, we consider the effective distance of each event to
the ID wall deff [22,36]. This observable is computed byfollowing the reconstructed direction of each event backwards from its reconstructed vertex to the ID wall, and canbe interpreted as the minimal distance needed 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 Sec. VA are shown in Fig. 10, for thereconstructed energy ranges corresponding to the modelindependent analysis, 7.5–29.5 MeV, and the spectral
analysis, 15.5–79.5 MeV. Comparing the deff distributions
from the data and the Monte-Carlo simulation allows us toestimate contributions from radioactivity near the wall, anddetermine a suitable deff cut. In this study we impose lower
thresholds on deff ranging from 3 to 5 m depending on thereconstructed energy:
Erec − 15.5 MeV
deff > max 300; 500 −× 50 cm:
1 MeV
2. Pre- and postactivity cuts
Atmospheric neutrino interactions can produce bothprompt signals from photons, electrons, or energetic heavy
particles, and delayed signals from decay of muons andpions to electrons. These signals will typically be separated
by a few microseconds and can therefore share the same

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

122002-12

DIFFUSE SUPERNOVA NEUTRINO BACKGROUND SEARCH AT …

PHYS. REV. D 104, 122002 (2021)

350
50000
300
Events/100 cm

Events/100 cm

40000

30000

250
200
150

20000
100
10000
50

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

0

0.8
0.6
0.4
0.2
0
500

1000

1500

2000
2500
deff [cm]

3000

3500

0
0.8
0.6
0.4
0.2
0
−0.2
−0.4

4000

500

1000

1500

2000
2500
deff [cm]

3000

3500

4000

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 number of events
observed in the data for deff > 1000 cm. The lower panels show the relative difference between the numbers of data and IBD simulation
events. We only consider events reconstructed in the FV and passing the noise reduction cuts. The assumed DSNB spectrum is takenfrom the Horiuchi þ 09 6 MeV model [16].

trigger window. Depending on which particle deposits themost light, unusually high activity can be noticed either
before or after the main peak. For preactivity we computethe maximal number of hits, N max
pre , 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 postactivity, using analgorithm developed for other analyses in SK and T2K thatcomputes the number of electrons from muon and pion
decay N decay-e [70]. The N decay-e distributions for low andhigh energy events are shown in Fig. 11. In the analyses
discussed here, we require N max
pre < 12 and N decay-e ¼ 0.
3. Cherenkov angle
The opening angle of the Cherenkov light cone (θC )emitted by highly relativistic particles (electrons and
600

500

CCQE-like
CC1π

1400

CC other
NCQE

1200

NC1π
NC other

1000
EventsEvents400

positrons in the current context) in pure water is around
42°. Conversely, at the Oð10Þ MeV energies considered
here, heavier particles like muons and pions will be
typically observed near the Cherenkov threshold, leading
to cones with smaller angles in the current analysis range.
In addition, NCQE interactions with multiple photon
emission can produce multiple overlapping Cherenkov
cones, which 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

300

CCQE-likeCC1π
CC other
NCQE
NC1π
NC other

800
600

200
400
100

0
0

200
1

Ndecay-e

2

3

0
0

1

Ndecay-e23FIG. 11. The N decay-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.

122002-13

K. ABE et al.

PHYS. REV. D 104, 122002 (2021)

250

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

1000
Events/2.7-degree

300
Events/2.7-degree

1200

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

350

200
150

800
600
400

100
200

50
0

20

30

40

50

60

70

80

0

90

20

30

40

θC [degree]

50

60

70

80

90

θC [degree]

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 and background distributions are normalized
to the same area.

to be in ½38; 50°. In the spectral analysis, we will also use
this angle to define sidebands to evaluate atmosphericbackgrounds.

The Lclear distributions from the atmospheric and IBD
Monte-Carlo simulations are shown in Fig. 13. Events withlower Lclear have fuzzier rings. For the analyses presented
here we require Lclear < 0.36.

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 and computethe opening angles of the cones formed by the directions of
all possible 3-hit combinations. We then identify the peakof this opening angle distribution θ0 and estimate the“clearness” of the ring by computing

We estimate the signal efficiencies of most of the cutsdetailed above using the IBD Monte-Carlo simulation.
CCQE-likeCC1π
CC other
NCQE
NC1π
NC other
IBD signal

300
250
Events/0.005

Events/0.005

ð4Þ

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

80

Energetic muons scatter less than electrons and hence
often deposit more charge in individual PMTs thanOð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

N triplets ðθ0  3°Þ
Lclear ¼:
N triplets ðθ0  10°Þ100

5. Average charge deposit

60

40

200
150
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 and background
distributions are normalized to the same area.

122002-14

DIFFUSE SUPERNOVA NEUTRINO BACKGROUND SEARCH AT …
120

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

80

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

350
300
Events/0.5

Events/0.5

100PHYS. REV. D 104, 122002 (2021)

60

250
200
150

40
100
20
0

50
1

1.5

2

2.5
Q50/N50

3

3.5

0

4

1

1.5

2

2.5
Q50/N50

3

3.5

4FIG. 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.

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 amonoenergetic electron beam produced by a linear accelerator (LINAC) at different positions inside the detector[32]. For a given cut, we compare the measured andpredicted efficiencies using a dedicated Monte-Carlo simulation, and define the systematic uncertainty as the
maximum discrepancy over all calibration tests. LINAC
calibration thus allows us to estimate uncertainties for ringclearness, the average charge, and Cherenkov angle cuts.
Conversely, for the incoming event cut—which requiresconsidering uniformly distributed events—we use a strategy developed for the SK solar neutrino analyses [36–39]:
we measure 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 forthese estimates are obtained from calibration studies using
a nickel source [31]. Finally, efficiencies for the pre- andpostactivity cuts can be directly estimated using a spallation-dominated sample of low energy events, with negligible systematic uncertainty. Overall, positron candidate
selection cuts allow us to remove a large fraction of
atmospheric backgrounds while keeping up to 85% of
the signal. The total systematic uncertainty computed for
positron candidate selection cuts using these procedures areof 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 andneutron preselection
As explained in Sec. II, most SHE-triggered events arefollowed by a 500 μs AFT window in order to record theneutron capture signal. Before attempting to identify
neutron candidates in this combined trigger window, weapply a cut on the maximum hits inside a 200-ns window:
N 200 < 50 to eliminate SHE þ AFT windows containing
muons. This cut is associated with a signal efficiency
larger than 99%. After applying the N 200 cut we scan the
SHE þ AFT 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 is within a few tens of centimeters of thewell-reconstructed IBD vertex. Since the SK resolution isabout 50 cm for Oð10Þ MeV events, this proximity allows
us to reliably estimate the required photon emission timefor each PMT hit to be the result of the neutron capture. A
neutron candidate is therefore defined as a group of hits thatis clustered in emission time. We then look for clusters
maximizing N 10 , the number of hits in a 10-ns window.
In previous SK-IV searches, neutron candidates were
defined as clusters of hits satisfying N 10 ≥ 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, weapply a cut to reduce the continuous PMT dark noise, likely
due to scintillation in the PMT glass, that causes a single
PMT to flash more than once in rapid succession with atimescale of ∼10 μs. To do so, for each PMT we remove
hit pairs separated by less than 6 ð12Þ μs for N 10 > 6(N 10 ≤ 6). After applying this cut, we define neutron
candidates as hit clusters verifying N 10 ≥ 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 the122002-15

K. ABE et al.

PHYS. REV. D 104, 122002 (2021)

continuous dark noise cut relies on investigating PMT hit
pairs, its efficiency sharply increases near the edge of theSHE þ AFT window. We restrict the neutron search window to ½14; 523 μs (½14; 373 μs for events with a 350-μs
AFT window). After these different 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 seven neutron candidates per event. To
further reduce this background, for each neutron candidate
22 discriminating variables are calculated. These variables
are related to specific features of noise events (such as
clustered PMT hits), the Cherenkov light pattern of theneutron candidate, and its position with respect to theprimary 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 adata set of 2.8 × 108 neutron candidates. In this sample,
about 2 × 106 candidates are neutron captures simulatedusing the procedure described in Sec. III B and the rest arecomposed of low energy accidental coincidences fromrandom trigger events. Table II summarizes the BDT
training parameters. The final BDT performance is shownin Fig. 15. Even accounting for preselection cuts, thisneutron tagging algorithm allows reduction of the background rate by a factor of 7 compared to the analysis
presented in Ref. [72] for the same signal efficiency.
3. Systematic uncertainty on neutron tagging efficiency
To estimate the systematic uncertainty on the BDT
efficiency associated with the 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 (241Am=Be) radioactive source
embedded in a bismuth germanium oxide (BGO) scintillator. A detailed description of this procedure is given in
Ref. [35]. To calculate the neutron identification efficiency
for data, we follow the methodology described in Ref. [76]and fit the timing distribution of tagged neutrons by adecaying exponential plus a constant term as shown in

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 inSK, as described in Sec. 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, are44.7% and 7, respectively. These values depend linearly on the
PMT noise [75], and hence vary by about 15% over the SK-IV
period.

Fig. 16. True neutron capture event counts decay exponentially in time from the detection time of the prompt event,while uncorrelated backgrounds are constant in time. By
comparing the efficiency in data and simulation for BDT cut
points relevant in our analysis, as described in Appendix B,
we assign 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, dedicated solar

TABLE II. Parameters used for the training of the neutron
tagging BDT in the XGBoost Python library framework.
Parameter

Value

Learning rate
Subsample
Max depth
Tree method
Max iterations
Best iteration

0.025220.97
10approx1500
1170

FIG. 16. Exponential þ constant fit of neutron candidate timings in the 2016 241Am=Be sample collected at the center of the
SK tank. The resulting time constant, 195.5  6.1 μs, is consistent with prior measurements.

122002-16

DIFFUSE SUPERNOVA NEUTRINO BACKGROUND SEARCH AT …
neutrino cuts are not implemented for the DSNB modelindependent analysis. For the spectral analysis, however,
events with no tagged neutrons are also considered in orderto maximize the effective exposure. When considering these
events, we hence apply additional cuts specifically targeting
solar neutrino backgrounds. These cuts are based on the
opening angle θsun between the reconstructed direction ofcandidate events and the direction of the Sun. The cosine ofthis angle is peaked at 1 for solar neutrino interaction
products, and is smeared by kinematic and reconstruction
effects. To estimate the impact of the latter, a specificobservable, the multiple scattering goodness (MSG) g has
been designed for the solar 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 in Ref. [22]and apply cos θsun cuts in the four MSG bins shown in
Fig. 17. We estimate the impact of these cuts on DSNB
signal events by assuming that these events have a uniformly distributed θsun , and we obtain the associated MSG
distributions using the IBD Monte-Carlo simulation.
Conversely, to estimate the number of solar neutrino events
above 15.5 MeV, we evaluate the number of solar events in
13.5–15.5 MeV from data and extrapolate the result tohigher energies using the solar neutrino Monte-Carlosimulation from Ref. [77], a procedure already used in
the DSNB search at SK-I,II,III [22]. Here, we can safely
assume the cos θsun distribution for a nonsolar event to beuniform [39], and therefore approximate the number ofsolar events in 13.5–15.5 MeV by
N solarlow ¼ N low ðcos θ sun > 0Þ − N low ðcos θsun < 0Þ;

ð5Þ

using the data passing the noise, spallation, and positron
selection cuts described in Secs. VA, V B, and V C,

PHYS. REV. D 104, 122002 (2021)

respectively. Using the solar neutrino Monte-Carlo simulation from Ref. [77], the predicted number of solar
neutrino events above 15.5 MeV will then be
N solar ≈ 0.12N solar
low :

ð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 us 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 cuts bycomparing 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 Sec. V C,
we estimate systematic uncertainties using results fromthe solar neutrino analysis [36–39] and the LINAC calibration [32]. The final cuts are the following:
Erec − 15.5 MeV
deff > max 300; 500 −× 50 cm;
1 MeV
N max
pre < 12;
N decay-e ¼ 0;
Lclear < 0.36;
q50 =n50 < 2;
38° < θC < 50°:
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 us to remove 80% of
the atmospheric neutrino backgrounds, and bring contamination from natural radioactivity to negligible levels.

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.

2. Spallation and neutron tagging cuts
Unlike positron candidate selection cuts, spallation andneutron tagging cuts need to be optimized in multiple

122002-17

K. ABE et al.

PHYS. REV. D 104, 122002 (2021)

Signal Efficiency [%]100

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

60

ring clearness cut
charge/hit cut
Cherenkov angle cut

40

neutron tagging (tagged-n =1)

20

0

10

15

20

25

30

Erec [MeV]

Signal Efficiency [%]

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 ≠ 1 tagged neutrons (bottom) are
required, as started after the SHE þ AFT requirement and thenoise reduction cuts. While the one tagged neutron region iscommon to both the model-independent and spectral analyses
above 15.5 MeV, the ≠ 1 tagged neutron region is used only inthe spectral analysis. Here, we applied noise reduction, spallation,
third reduction, neutron tagging, and solar cuts sequentially andshow 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–79.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.

energy regions to account for the important variations of thespallation and atmospheric background rates in the analysis
windows. While only events with one tagged neutron areused in the DSNB model-independent analysis, the spectral
analysis also considers events with zero or more than one
tagged neutron. Here, we first optimize spallation andneutron tagging cuts for the signal region where exactly one
tagged neutron is required. Using the thus optimized
neutron tagging cuts, we then derive optimal spallation
cuts for events with ≠ 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 Erec bins.
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 backgroundonly hypothesis. We define the sensitivity as the predicted
number of signal events after cuts over the 90% CL upper
limit on the number of signal events, computed using theRolke method [78]. To simplify the interpretation of theresults presented in this paper, we apply the same cuts
for both the model-independent and spectral analyses in the
15.5–29.5 MeV reconstructed energy range, where thetwo analyses overlap. In particular, we require spallationand solar neutrino backgrounds to be negligible above
15.5 MeV. The remaining rates of these backgrounds canbe reliably evaluated by rescaling the number of events
observed after noise reduction, spallation, and positroncandidate selection cuts by the neutron 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
uncertainties, computed as described in Secs. V B, V C,
and V D, in Table III.
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, which
covers most of the spectral analysis window, spallation andsolar 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 are still negligible and spallation backgrounds are mostly linked to short-lived isotopes
that are easy to reject using the cuts from Sec. V B.
For these energies, we therefore apply the spallation cut
derived above, whose efficiency is close to 90%. In the
15.5–19.5 MeV range, however, both solar neutrino andspallation backgrounds largely dominate and need to be
reduced using the targeted reduction strategies described inSecs. V B and V E. In this analysis, we remove the solar
TABLE III. Systematic uncertainties for the reduction cuts,
computed as described in Secs. V B, V C, and V D. The
systematic uncertainties for the AFT, noise reduction, and preand postactivity cuts are negligible. The solar neutrino cut is used
only for the spectral analysis, when events have ≠ 1 tagged
neutrons.

122002-18

Cut
Spallation cut
deff cutLclear cut
q50 =n50 cut
θC cutNeutron tagging
Solar neutrino cut

Relative systematic error
0.1%
0.1%
0.2%
1.2%
0.7%
12.5%
1.0%

DIFFUSE SUPERNOVA NEUTRINO BACKGROUND SEARCH AT …

Erec (MeV)

15.5–16.5

16.5–17.5

17.5–18.5

18.5–19.5

g < 0.4
0.4 < g < 0.5
0.5 < g < 0.6
g > 0.6
Efficiency

0.05
0.39
0.59
0.73
73.1%

0.35
0.61
0.73
0.79
82.2%

0.45
0.77
0.81
0.91
88.3%

0.93
0.93
0.93
0.93
96.6%

neutrino backgrounds using the cuts in Sec. V E with thecut points from Ref. [22]. The cut points and the corresponding signal efficiencies are shown in Table IV for thedifferent MSG 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.
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 8B and 8Li, that have half-lives
close to 1 sec. We then optimize the spallation likelihood
cuts using the method described in Appendix A 4 tocompute the remaining
amount of 8B and 8Li decays,
pﬃﬃﬃﬃ
and maximize S= B where B and S are the numbers ofspallation and nonspallation events, respectively. The
efficiencies of the optimal cuts are of 65% and 63% in
15.5–17.5 MeV and 17.5–19.5 MeV, for a 8B þ 8Li
remaining fraction of about 2.5%. The corresponding
remaining number of spallation events, as well as theperformance of this cut, will be discussed and compared toresults from the previous SK DSNB analyses in Sec. VII E.

6

10

Noise reduction
Spallation
deff
Positron ID

5

Events/2-MeV

10

4

10

3

10

2

10

10
10

15

20

25

30

Erec [MeV]

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

the same neutron tagging cut, with a 20% signal efficiency,
for all energies. While most spallation cuts are removed inorder to focus on neutron tagging, multiple spallation cuts
are conserved in order to avoid contaminating the neutron
detection window with isotope decays—an effect notaccounted 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 Sec. V F are applied, we find
102 events with exactly one neutron in 7.5–79.5 MeV.
For these events, the distribution of the time 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 241Am=Be calibration, yields a
700
SK-IV Data
Atmospheric-ν (non-NCQE)
Atmospheric-ν (NCQE)
9
Spallation Li
Reactor-ν
Accidental coincidence

600
500
Events/2-MeV

TABLE IV. Cuts on cos θsun (upper limits), and signal efficiencies for SK-IV. These cuts are only applied to events with ≠ 1
tagged neutron, in the spectral analysis. They depend on the valueof the multiple scattering goodness, g, and on the reconstructedkinetic 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.

PHYS. REV. D 104, 122002 (2021)

400
300

G. Selected events in the final data samples

200

Here we discuss the behavior of the observed data after
the cuts. In particular, the impact of spallation, deff , and theother positron candidate selection cuts on data are shown in
Fig. 19. As shown in this figure, for reconstructed energies
below 20 MeV, contributions from spallation and radioactivity near the wall are significant. To validate ourbackground modeling and neutron tagging procedure wealso compare the observed data to the background expectations after the noise reduction, multiple spallation cuts,
DSNB positron candidate selection, and neutron tagging
cuts. To simplify the interpretation of the data, we impose

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 a uniform
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 9Li decays. The systematic
uncertainties are shown by the hatched lines.

122002-19

K. ABE et al.

PHYS. REV. D 104, 122002 (2021)

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 solidorange curve shows the fit result by an exponential plus a
constant, giving a time constant of τ ¼ 154  104 μs.

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, which will be used for the spectral
analysis, we also find 248 events with zero neutrons andnine events with more than one neutron. The observation
times and reconstructed vertices for these events are shownin Figs. 22 and 23. No event time cluster has been observed
and the time dependence of the event rate is consistent withthe live time variations over the SK-IV period. Events arealso uniformly distributed all over the FV, irrespective of
the number of neutrons. The space and time distributions of
the selected events are therefore consistent with what is
expected for DSNB candidates.

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

VI. 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 inthis analysis, model their energy dependence, and estimate
the corresponding uncertainties.
A. Background estimation
1. 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 chargedcurrent (CC) or neutral-current (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, which has a similar backgroundcomposition. The number of events in this sideband is
88.8% of the atmospheric neutrino Monte-Carlo prediction,
and we therefore use this factor to rescale the non-NCQE
backgrounds. The associated uncertainty is due to statisticalfluctuations in the number of events in the sideband andis 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 Monte-Carlo simulation
predictions of neutrino and antineutrino NCQE interactions
are renormalized using the scaling factors from thismeasurement: f ν-NCQE ¼ 0.80  0.08ðstatÞþ0.24−0.18 ðsystÞ andf ν̄-NCQE ¼ 1.11  0.18ðstatÞþ0.29
ðsystÞ,respectively.
Note
−0.22
that this rescaling provides an inclusive estimation of the
NCQE and NC 2p2h interactions. The errors of thesefactors are taken as cross section uncertainties.

122002-20

DIFFUSE SUPERNOVA NEUTRINO BACKGROUND SEARCH AT …

PHYS. REV. D 104, 122002 (2021)

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

Two additional flux uncertainties are considered for
the NCQE background: the atmospheric neutrino flux
uncertainty and the uncertainty due to the flux difference
between the T2K beam and atmospheric neutrinos. For theformer, 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 and then the effect ofcross section model uncertainties is different. To estimate
this uncertainty, the ratios of flux-averaged events in theT2K beam and atmospheric neutrino fluxes for different
cross section models are calculated. Here, six models are
considered: Spectral Function [79,80], Relativistic Mean
Field [81–84], Superscaling [85,86], Relativistic Green’s
Function [84,87–89] with two functional forms for the
potential (EDAI, Democratic), and Relativistic Plane Wave
Impulse Approximation [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 efficiency
and neutron multiplicity. For the uncertainty on the tagging efficiency, 12.5% is employed in this analysis as the
maximum difference between the measured and predictedefficiencies in the 241Am/Be calibration. The neutron tagging
efficiency will also depend on the distance between theneutron 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 inmodeling the neutron travel distance on the efficiency is ofaround 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 us tocompare observations with Monte-Carlo simulation predictions for neutron multiplicity as a function of the reconstructed squared momentum transferred to 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 eventswith the tagged neutron being one, both for neutrinos andantineutrinos.

122002-21

K. ABE et al.

PHYS. REV. D 104, 122002 (2021)

25

Sec. V. The difference between this result and the nominal
prediction is then taken as the systematic uncertainty
associated with the γ emission modeling. For this study,
we set the smearing parameter σ to 3 MeV for all energies;
this value allows us to cover the discrepancy observed
between the predicted and observed θC distributions in theT2K 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 thenvaries from 60% to 80% for low to high energies. The
predicted NCQE spectrum with the T2K scaling is shownin Fig. 24.

SK-IV Data
Atmospheric-ν (non-NCQE)
Atmospheric-ν (NCQE)
Spallation Li

20
Events/2-MeV

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

15

10

5

0

10

20

30

40

50

60

70

Erec [MeV]

2. Spallation 9Li

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

The background rate of spallation 9Li can be estimatedby the following formula:

Spallation Li
Reactor-ν

10
Events/2-MeV

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

B9Li ¼ R9Li × T live × 22.5 kton
× Br½9Li → β þ n × f window × ϵreduc ;
1

10

20

30

40

50

60

70

Erec [MeV]

FIG. 24. Reconstructed energy spectra after data reductions forthe background expectation and observation, including the signal
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 hatched lines inthe signal window representing the total systematic uncertainty.
The red dashed line represents a DSNB signal expectation fromthe Horiuchi þ 09 model [16] shown only for the signal window.

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 thesecondary nuclear reaction. The modeling of these γ
emissions is in particular a probable cause of the currently
observed discrepancy between the observed and predicted
Cherenkov angle distributions using neutrinos from theT2K beam [48]. One notable effect of this mismodeling is asmearing of the reconstructed energy analogous to a
resolution effect. To evaluate the associated systematic
uncertainty we model this smearing by convolving theoriginal spectrum with a Gaussian distribution and using
the energy-dependent ratio between the smeared andunsmeared spectra to renormalize the Monte-Carlo events.
We then compute the number of renormalized events ineach energy bin after the reduction steps described inð7Þ

where R9Li is the production 9Li rate measured at SK(0.86  0.12ðstatÞ  0.15ðsystÞ kton−1 · day−1 ) [59], T liveis the operational live time (2075.3 and 2970.1 days
below and above Erec ¼ 9.5 MeV, respectively), Br½9Li →β þ n ¼ 0.508 is the branching ratio for β þ n decay,
f window is the fraction of the 9Li decay energy spectrum
above the search energy threshold, and ϵreduc is the
reduction efficiency. We estimate f window using theMonte-Carlo simulation elaborated for IBD events, renormalized to fit the 9Li spectrum. This same simulation can beused 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 Sec. V. The
dominant systematic uncertainty on the number of 9Li
events in each energy bin after cuts arises from estimating
the impact of spallation cuts. Indeed, as discussed inAppendix A 4, our method assumes that the time difference
between an isotope decay and its parent muon is the onlyisotope-dependent spallation observable. In this study, we
account for the impact of other spallation observables bydefining a 50% uncertainty on the number of 9Li events.
The uncertainty on the event rate is taken from the previous
SK study [59] as 22%. Other uncertainties come from thereduction, especially from the spallation cut and neutron
tagging, that is 20% in total. The total uncertainty assigned
for the 9Li background estimation is therefore 60%.
3. Reactor neutrinos
Reactor neutrino backgrounds are estimated by renormalizing the IBD Monte-Carlo simulation. They populate

122002-22

DIFFUSE SUPERNOVA NEUTRINO BACKGROUND SEARCH AT …
only in the lowest energy bin, as seen in Fig. 24. In thisanalysis, a conservative 100% uncertainty is assigned totheir 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 reduction,
spallation, and positron reduction cuts N data
pre-ntag , and rescaling it by the neutron mistag rate:
Bacc ¼ ϵmis × N datapre−ntag :each bin, as shown in Table V. Here the most significant
bins are found to be at the 2σ level. We therefore conclude
that no significant excess is observed in the data over the
background prediction in any energy bin, and place upper
limits on the extraterrestrial ν̄e flux using the pseudoexperiments above. Here we obtain the 90% CL upper bound on
the number of signal (N limit
90 ), as an excess of the observation over the background expectation, by varying thenumber of observed events in each energy bin by theirstatistical uncertainties, as well as varying the number of
expected background events by their systematic uncertainties from each source. The 90% CL upper limit on the ν̄eflux is then calculated as
ϕlimit90 ¼

ð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 Sec. III B over theentire SK-IV period into ten 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. Differential upper limit
The observed data and expected background spectra areshown in Fig. 24. The corresponding p-values are calculated for each energy bin by performing pseudoexperiments
as follows. We vary the number of expected background
events randomly based on Gaussian distributions with theirwidths being the 1σ systematic uncertainties from eachbackground source. We then calculate p-values from theresulting distributions and the observed number of events in

PHYS. REV. D 104, 122002 (2021)

N limit90;t · N p · σ̄ IBD · ϵsig

ð9Þwhere t is the operational live time [sec], N p is the numberof free protons in the FV of the SK, σ̄ IBD is the IBD crosssection [10−41 cm2 ] at a mean neutrino energy in thecorresponding region (Ēν ), and ϵsig is the signal efficiency.
Note that the neutrino energy is obtained as Eν ¼ Erec þ
1.8 MeV in IBD. For the expected sensitivity, the sameprocedure is applied while the number of observed events is
replaced with the number of nominal background predictions. Here statistical uncertainties on the backgrounds are
considered instead. The expected sensitivities and observed
upper limits from this work are summarized in Table V.
These results are also compared with the previously
published results and theoretical predictions in Fig. 25.
Note that the SK-I,II,III limits presented in this figure havebeen obtained using a different methodology, more similar
to the one used for the spectral analysis in Sec. 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 for the model-independent differential
limit. The sensitivities obtained with this analysis are theworld’s tightest above neutrino energies of 11.3 MeV.

TABLE V. Summary on the 90% CL expected sensitivities and observed upper limits as well as the corresponding
p-values in each electron antineutrino energy bin (Eν ¼ Erec þ 1.8 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

Expected (cm−2 sec−1 MeV−1 )1Observed (cm−2 sec−1 MeV−1 )1

3.71 × 10
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−2

4.44 × 10
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

122002-23

p-value
0.346
0.886
0.938
0.830
0.243
0.686
0.073
0.597
0.051
0.605
0.126

K. ABE et al.

PHYS. REV. D 104, 122002 (2021)

SK-IV 2970 days, Observed 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

2

νe Flux Upper Limit [/cm /sec/MeV]

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

10

1

10−1
10

15

20
25
νe Energy [MeV]30

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

The current limit disfavors the Totani þ 95 model [8] andthe most optimistic predictions of the Kaplinghat þ 00
model [9], and is reaching close to several other model
predictions. Due to not only a higher exposure but also
higher cut efficiencies and more precise backgroundestimation, this new analysis considerably improves on
the previous SK-IV DSNB model-independent search [71].
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 mostof the signal and the irreducible backgrounds, and two
sidebands with θC ∈ ½20; 38° and θC ∈ ½78; 90°. The low
Cherenkov angle region will be populated with mostly
atmospheric backgrounds involving visible muons andpions while the high angle region will be mostly populatedby 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 efficienciesof the neutron tagging cuts the non-IBD-like region is
expected to contain a sizable amount of signal. Our analysis
will hence involve six regions of parameter space: two
signal regions with intermediate values of the Cherenkov
angle and four sidebands with low and high Cherenkov
angle values, as summarized in Table VI.
B. Spectral shape fitting
We perform a simultaneous fitting of the signal andbackground spectra to the observed data in all six regions ofparameter 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

122002-24

DIFFUSE SUPERNOVA NEUTRINO BACKGROUND SEARCH AT …
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 dominated by
the interactions producing visible muons and pions and by NCQE
interactions, respectively. We assign numbers for each region.
N tagged-n
θC

20–38°

38–50°

78–90°

1
≠1

I
IV

II
V

III
VI

spectrum can be reliably predicted by the IBD Monte-Carlosimulation for a given DSNB model, the treatment of theatmospheric neutrino and spallation backgrounds is more
complex. For this study, we follow the method described inRef. [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 of invisible
muons and pions. The energy distribution of these electrons
follows a Michel spectrum, whose shape is 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 the Michel spectrum at SK directlyfrom the data, using a sample of electrons produced bycosmic ray muon decays (a similar sample for atmospheric
and accelerator neutrino oscillation and proton decay
analyses is described in Refs. [90,91]). We then use theatmospheric neutrino Monte-Carlo simulation to computethe fractions of background events in the different signal
and background regions. Since electrons cannot be distinguished from positrons at SK, this background dominates in the signal 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 associatedspectral shapes in all regions using the atmospheric neutrinoMonte-Carlo simulation. Similarly to Michel electrons, thisbackground 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 shapes byconsidering an atmospheric Monte-Carlo sample with only
CC interactions, visible muons and pions, and no electrons.
d. NCQE interactions: Unlike the other types ofatmospheric neutrino events, the NCQE spectrum peaks

PHYS. REV. D 104, 122002 (2021)

at low energies, similarly to the DSNB spectrum. Here
we define spectral shapes using simulated NCQE eventswith no electrons. This background, often involving multiple γ rays, will dominate in the high Cherenkov angle
regions III and VI, but will also contribute to all otherregions.
e. Spallation backgrounds: These backgrounds are
negligible in the one-tagged-neutron regions I, II, and III
after the cuts. However, in the other regions, especially insignal region V, residual spallation backgrounds from 8Li,
8B, and 9C decays could remain. Since 9C has a shorter
half-life and is produced closer to the muon track than theother two isotopes, it is expected to be subdominant.
To model the spallation spectrum and its variability weparametrize this spectrum and vary its isotope composition
as follows:
Sspall ðEÞ ¼ f 0 S9C ðEÞþ ð1 − f 0 Þðf 1 S8B ðEÞ þ ð1 − f 1 ÞS8Li ðEÞÞ;with f 0 ∈ ½0; 0.5 and f 1 ∈ ½0; 1;

ð10Þ

where Si ðEÞ refers to the decay spectrum for a given
isotope i, normalized to one in the spectral analysis energy
window, and f 0 and f 1 are the isotope fractions. Here, weconstrain 9C to always be less abundant than the other two
isotopes after spallation cuts are applied. In what follows,
we will take the nominal spallation spectrum to be theaverage between the steepest and the least steep spectra andparametrize possible deviations from this hypothesis assystematic uncertainties. We evaluate the β spectra using
external measurements [92] for 8B decays and GEANT4 for
8Li and 9C. As with 9Li, we then use these spectra torenormalize events in the IBD Monte-Carlo simulation andthus model detector resolution effects.
Although the first four background categories do
not include all possible types of atmospheric neutrino
interactions, their associated spectra can be linearly
combined to model other background categories. Pionproducing NC backgrounds can notably be treatedas a superposition of NCQE and decay electron backgrounds. The five categories described here hence allow
us to map all the background configurations relevant tothis analysis.
For each background category we define PDFs whichmodel the corresponding energy spectra in each of the six
Cherenkov angle and neutron regions. These PDFs arenormalized to 1 across all regions. We apply a similar
procedure to obtain signal spectral PDFs associated todifferent DSNB models, using the IBD Monte-Carlosimulation. For N events observed events with energies
fE1 ; …; EN events g, we then extract the most probable numbers of signal and background events by performing an
extended unbinned maximum likelihood fit, using thefollowing likelihood function:

122002-25

K. ABE et al.

PHYS. REV. D 104, 122002 (2021)−LðfN j gÞ ¼ e

P5

N
j¼0 j

NY
5
events X

ðrÞ

N j PDFj ðEi Þ:

L0 ¼

ð11Þ

Z þ∞

i¼1 j¼0Here, the index j ¼ 0; …; 5 refers to the signal and the fivebackground categories described above. Then fN j g designates the numbers of events for the signal and thedifferent background types, across all six Cherenkov angleðrÞand tagged neutron regions. For each event, PDFj ðEi Þ isthe 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 Eq. (11) is that the sizes of the different background contributions can be treated as nuisance parameters.
This likelihood, however, does not account for the considerable uncertainties on the spallation and atmospheric background spectral shapes in the energy window of 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 accountfor systematic uncertainties on the spallation and atmospheric neutrino background spectral shapes, energy
reconstruction, and reduction efficiencies by convolving
the likelihood introduced in Eq. (11) with suitable probability distributions. Here we neglect the systematic uncertainties associated with the μ=π background. The amounts ofthese backgrounds 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
shape of the background associated with invisible muon andpion decays since this shape has been directly extracted froma large sample of SK-IV data. We then consider five sourcesof systematic uncertainties: energy scale and resolution,
spectral shape predictions for the spallation, NCQE, νe CC
backgrounds, and reduction cut efficiencies.

−∞2e−ϵ =2 LðϵÞdϵ;ð12Þ

where LðϵÞ represents the likelihood from Eq. (11) with the
PDFs distorted by an amount ϵ. These uncertainties have anegligible impact on the final result and are therefore nottaken into account in the final fit.
While energy scale and resolution uncertainties on signal
and atmospheric background spectra can be neglected inthis analysis, the spallation background spectrum is particularly steep and hence more sensitive to the modeling ofenergy reconstruction. In this analysis, we will therefore
account for the associated modeling uncertainties whenestimating 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 8B, 8Li, and 9C fractions as shown in
Eq. (3). Then, we distort the steepest and least steep spectra
found using this method in order to account for the 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, weparametrize the total spectral distortion associated with thesystematic uncertainty by using a third order polynomial
P 3 such that
PDFnew ¼ PDFold × N½1 þ ϵP 3 ðEÞ;

ð13Þ

where ϵ is expressed in units of σ. More details about theseestimates are given in Appendix D. We then convolve thelikelihood L from Eq. (11) with a Gaussian prior on ϵ:
0

L ¼

Z 2
−2

ϵ2

LðϵÞe− 2 dϵ:

ð14Þ

1. Energy scale and energy resolution

3. Atmospheric νe CC spectral shapes

To estimate the impact of energy scale and energy
resolution uncertainties on the signal and atmosphericbackground spectra, we follow the procedure describedin Ref. [22]. The energy scale and energy resolution
uncertainties are obtained from the solar neutrino analysis
and increased 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 Eq. (11) using the energy resolutionfunctions from Refs. [36–39]. We parametrize the amountof distortion applied to these PDFs using the variable ϵ suchthat ϵ ¼ 1 represents a 1σ deviation in the energy scale andresolution. The final likelihood can then be expressed as

The shape of the νe CC background for electron energies
below about 100 MeV has not been studied in-depth todate, although Monte-Carlo simulations predict a linear
increase in the number of events with the electron energy.
In our analysis, we model this spectral shape uncertainty byallowing the slope of the νe CC spectrum to vary in the two
signal regions. We then integrate the distorted PDF using
pseudo-Gaussian weights. A detailed description of thisprocedure is provided in Ref. [22].
At this step of the analysis, we assume the spectral shape
distortions in the ¼1 and ≠1 tagged neutron regions to be
fully correlated and model the effect of neutron multiplicity
uncertainties separately.

122002-26

DIFFUSE SUPERNOVA NEUTRINO BACKGROUND SEARCH AT …
4. Relative normalization for NCQE events
The Cherenkov angle distribution of atmospheric NCQE
events is associated with large uncertainties stemmingfrom the modeling of γ ray emission. We incorporate thisuncertainty 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 thisnormalization to correspond to a 100% change of thenumber of NCQE events in the intermediate Cherenkov
angle region. Additionally, in order for the overall 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 − f high ϵÞ;

ð16Þ

NCQE
PDF’NCQE
low;n ¼ PDFlow;n × ð1 − f low ϵÞ;ð17Þ

PHYS. REV. D 104, 122002 (2021)

conservative value of 40%. We then distort the PDFs in
the ¼ 1 and ≠ 1 tagged neutron regions for each background category as follows:
IBD
PDFIBDj;new ¼ PDFj;old × ð1 þ αj ϵÞ;ð18Þ

¼ PDFnon-IBD
× ð1 − αj fϵÞ;PDFnon-IBD
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 thisuncertainty into the final likelihood by integrating theinitial likelihood similarly to Eq. (14).
6. Combining the background uncertainties

where the “low,” “med,” and “high” indices refer to the low(20–38°), medium (38–50°), and high (78–90°) Cherenkov
angle ranges from Table VI and n ¼ 0, 1 for the tagged
neutron regions, respectively. f low ; f med ; f high are scaling
factors introduced to preserve the overall normalization of
the PDFs. The ϵ parametrizes the PDF distortion in units ofσ. We finally inject these distorted PDFs into Eq. (11) andconvolve the new likelihood with a weighted Gaussianprior, as with νe CC background uncertainties.

The three categories of systematic uncertainties on backgrounds described above—spallation and νe CC spectral
shapes, NCQE Cherenkov angle distribution, and neutron
multiplicity—are incorporated into Eq. (11) using themodified likelihood defined in Eq. (14) and the modified
PDFs defined in Ref. [22] as well as in Eqs. (17) and (19)with distortion parameters ϵνe CC , ϵNCQE , and ϵn , respectively.
The final likelihood Lsyst can then be written as
Lsyst ðfN j gÞ ¼Z 3−1dϵνe CCZ 3−1

dϵNCQE

Z 3
−2

dϵn

Z 2
−2

dϵspall× LðfN j g; ϵνe CC ; ϵNCQE ; ϵn ; ϵspall Þ
ϵ2n

ϵ2
spall

e− 2 e− 2× Gðϵνe CC Þ × GðϵNCQE Þ × pﬃﬃﬃ × pﬃﬃﬃ ;ππ

5. Uncertainty on the neutron multiplicity
The fractions of signal and background events observedin the ¼1 and ≠1 tagged neutron regions depend on boththe neutron tagging efficiency and the true neutron multiplicity of the observed events. As shown in the 241Am/Be
calibration study, the efficiency of the neutron tagging cut
for the DSNB spectral analysis is associated with asystematic uncertainty of around 6%. Conversely, theactual neutron multiplicity for atmospheric neutrino interactions is poorly known. As discussed in Sec. 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.
We model the effect of neutron multiplicity uncertainties
by varying the relative normalization of the PDFs in the ¼1
and ≠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 the systematic
uncertainty associated with μ=π backgrounds, which areassociated with negligible contributions in the signal
region. As discussed above, we set the 1σ variation of
the fraction of events with one tagged neutron to a

ð20Þ
where LðfN j g; ϵνe CC ; ϵNCQE ; ϵn ; ϵspall Þ is the likelihood
defined in Eq. (11) with distorted PDFs and G is theweighted Gaussian introduced in Ref. [22].
7. Energy-independent efficiency systematic error
Maximizing the likelihood Lsyst defined in Eq. (20) overthe numbers of events in all five background categories
results in a marginalized likelihood Lopt ðN s Þ where N s isthe number of signal events. Using this likelihood toconstrain the total number of DSNB events that passed
through the detector R requires rescaling N s by the averagevalue of the signal efficiency over all analysis regions andtheir associated energy windows, ϵsig . We therefore need toincorporate the systematic uncertainty on ϵsig into Lopt ,following a procedure similar to the treatment of 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 summarized in

122002-27

K. ABE et al.

PHYS. REV. D 104, 122002 (2021)

TABLE VII. Uncertainties associated with the fiducial volume
cut, the Strumia-Vissani IBD cross section calculation, and thelive time 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 phase
Fiducial volume
Cross section
Live time

I

II

III

IV

1.3%
1.0%
0.1%

1.1%
1.0%
0.1%

1.0%
1.0%
0.1%

1.5%
1.0%
0.1%

Table III. We also incorporate the uncertainties associatedto the fiducial volume cut, IBD cross section, and live
time calculation, which have been taken from the solar
neutrino analyses and are summarized in Table VII. Adding
all these uncertainties in quadrature, we then rewrite thelikelihood as

simulations used in Ref. [22]—with events generated inEeþ 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, forwhich DSNB fluxes are available only up to antineutrinoenergies of 50 MeV, we assume the predicted flux to be
zero in the 50–81.3 MeV range. This assumption is notexpected 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 inRef. [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 thespectra and associated systematic uncertainty using thesame procedure as at SK-IV. We then combine the four fit
results by summing the corresponding likelihoods maximized over the nuisance parameters.
E. Results

L0 ðRÞ ¼

Z 1
0−

Lopt ðϵRÞeðϵ−ϵsig Þ2
2σ 2
sig

dϵ;

ð21Þ

where σ sig is the 1σ systematic uncertainty on ϵsig .
Systematic uncertainties associated with reduction efficiencies weaken the final limits by around 1%. Similarly,
while the uncertainties associated with spallation andatmospheric neutrino backgrounds are of order 100%, theiroverall effect on the final limits is of about 5%. The limited
impact of spectral shape uncertainties is due to the large
dominance of backgrounds associated with invisible muonand pion decays, whose spectrum is well-known. The final
analysis will hence be largely dominated by statisticaluncertainties, which will allow us to easily combine SK-IV
observations with results from previous SK phases.
D. Combination with SK-I,II,III
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],which performed a DSNB spectral analysis over the22.5 × 2853 kton · days of data corresponding to the first
three SK phases. While this previous analysis did notinclude neutron tagging and used slightly different background categories and reduction cuts, it is also dominatedby statistical uncertainties. Hence, the SK-I,II,III results canbe readily combined with the SK-IV observations.
In order to study a wide variety of the discrete andparametrized models introduced in Sec. 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 theWe 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
shown in Figs. 26 and 27 for the DSNB signal and for the
spallation and atmospheric neutrino background categories
introduced in Sec. 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, which strongly depends on
the background composition. As expected, the dominantbackgrounds in the low and high Cherenkov angle sidebands are backgrounds with visible muons and pions and
NCQE backgrounds, respectively. While most of the signal
region is dominated by backgrounds from invisible muonsand pions, νe CC contributions become sizable above about
50 MeV. For SK-IV, the background rates are extremely
low in the IBD-like region, which will hence be particularly
sensitive to IBD signals in spite of the low efficiency of theneutron 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 numberof spallation events for each phase never exceed two 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, thespallation cut effectiveness at SK-IV should remain on parwith the ones observed at previous SK phases. Yet, an
important amount of spallation backgrounds (∼8 events) ispredicted below 19.5 MeV at SK-IV. The origin of thisimportant spallation background will be further studied infuture SK analyses with an enhanced performance by

122002-28

DIFFUSE SUPERNOVA NEUTRINO BACKGROUND SEARCH AT …

PHYS. REV. D 104, 122002 (2021)

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 the two signal regions and four sidebands in Table VI.

gadolinium. Note that this effect was not observed in otherSK searches probing a similar 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 can be
explained by statistical fluctuation, it led us to investigate
possible nonspallation sources for the excess observed overatmospheric neutrino background predictions in 15.5–
19.5 MeV. We first ruled out the possibility of a strong
DSNB signal, as it would be associated with a steeply
falling spectrum and hence to an unrealistically highsupernova rate. The hypothesis of an isolated astrophysicalevent is also strongly disfavored since the observed eventsare uniformly distributed over the SK-IV period, as discussed in Sec. V G. A radioactive origin for the observed
events is also highly unlikely, given their high energies andtheir uniform spatial distribution. The event positions arenotably not correlated with the locations of calibration

sources. Residual spallation therefore remains the mostplausible 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, with a total
exposure of 22.5 × 5823 kton · days. The associated limits
on the DSNB flux are shown in Tables VIII and IX. 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% CL 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 SK-I,
II,III,IV analysis. The observed upper limits on the DSNB
flux for SK-IV and for the combined analysis are of

122002-29

K. ABE et al.

PHYS. REV. D 104, 122002 (2021)

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.

the background prediction is higher than the 0.9σ excessfrom the previous combined analysis [22]. While notstatistically significant, this result is compatible with awide range of DSNB predictions with a flux comparable tothe one of the Horiuchi þ 09 model [16].
VIII. DISCUSSION
A. Constraints on the DSNB parameter space

FIG. 27. Best-fit signal and background spectra for SK-I (top), II(middle), and III (bottom), assuming the DSNB spectrum predictedby the Horiuchi þ 09 model [16]. The three regions presented here
correspond to different ranges 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.

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

Over the two search strategies presented in this paper, theanalysis described in Sec. VI gives a model-independent
differential limit on the ν̄e flux while the model-dependent
spectral analysis uses data from all phases of SK andyields the tightest constraints on DSNB models. In thissection, we present the results of this analysis using both
discrete DSNB models—most of them based on supernova
simulations—and simplified parametrizations.
The 90% CL expected and observed upper limits on the
DSNB rate and flux for theoretical models are shown in
Fig. 29 for the combined SK-I,II,III,IV analysis, with thecorresponding results being tabulated in Tables VIII and
IX. With an exposure of 22.5 × 5823 kton · days, the sensitivity of the combined analysis at 90% CL 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 starformation history allowed by observations. Additionally,
the Kresse þ 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σ hasbeen observed in the data and the resulting 90% CL limits
on the DSNB flux improve on the constraints obtained in

122002-30

DIFFUSE SUPERNOVA NEUTRINO BACKGROUND SEARCH AT …

PHYS. REV. D 104, 122002 (2021)

TABLE VIII. Best-fit values and the 90% CL upper limits on the DSNB fluxes (in cm−2 · sec−1 ) for the theoreticalmodels for phases SK-I to IV as well as for the combined analysis. Here the upper limits are given forEν > 17.3 MeV. For the Kresse þ 21 models, the “High,” “Fid,” and “Low” predictions correspond to the “W20BH2.7-α2.0,” “W18-BH2.7-α2.0,” and “W18-BH2.7-α2.0-He33” models from Ref. [12], respectively.
Best-fit
Model
Totani þ 95 ConstantKaplinghat þ 00 HMA (max)Horiuchi þ 09 6 MeV, max
Ando þ 03 (updated 05)
Kresse þ 21 (High, NO)
Galais þ 09 (NO)Galais þ 09 (IO)
Horiuchi þ 18 ξ2.5 ¼ 0.1
Kresse þ 21 (High, IO)
Kresse þ 21 (Fid, NO)
Kresse þ 21 (Fid, IO)
Kresse þ 21 (Low, NO)
Tabrizi þ 21 (NO)Kresse þ 21 (Low, IO)
Lunardini09 Failed SN
Hartmann þ 97 CE
Nakazato þ 15 (max, IO)
Horiuchi þ 18 ξ2.5 ¼ 0.5
Horiuchi þ 21
Malaney97 CGINakazato þ 15 (min, NO)

90% CL limit

SK4

All

SK1

SK2

SK3

SK4

All

Pred.

2.5þ1.4−1.3
2.6þ1.5−1.3
2.6þ1.4−1.3
2.7þ1.5−1.4
2.7þ1.5−1.3
2.5þ1.4−1.3
2.6þ1.4−1.3
2.6þ1.4−1.3
2.7þ1.5−1.3
2.7þ1.5−1.3
2.7þ1.5−1.3
2.7þ1.5−1.4
2.7þ1.5−1.3
2.7þ1.5−1.4
2.8þ1.5−1.4
2.6þ1.4−1.3
2.7þ1.5−1.4
2.7þ1.5−1.4
2.1þ1.3−1.2
2.7þ1.5−1.3
2.8þ1.5−1.4

1.3þ0.9−0.9
1.3þ0.9−0.9
1.3þ0.9−0.9
1.4þ0.9−0.9
1.4þ0.9−0.9
1.3þ0.9−0.9
1.3þ0.9−0.9
1.4þ0.9−0.9
1.4þ0.9−0.9
1.4þ0.9−0.9
1.4þ0.9−0.9
1.4þ0.9−0.9
1.4þ0.9−0.9
1.4þ0.9−0.9
1.4þ0.9−0.9
1.3þ0.9−0.9
1.4þ1.0−0.9
1.3þ0.9−0.9
1.2þ0.9−0.9
1.3þ0.9−0.9
1.4þ1.0
−0.9

2.3
2.3
2.4
2.3
2.3
2.3
2.3
2.4
2.3
2.3
2.3
2.3
2.4
2.3
2.4
2.3
2.4
2.2
3.4
2.3
2.3

6.3
6.7
6.0
6.6
6.7
6.3
6.4
6.1
6.7
6.8
6.8
6.8
6.6
6.8
6.8
6.5
6.5
7.1
4.3
6.8
6.8

7.0
7.1
7.0
7.2
7.2
7.0
7.0
7.1
7.1
7.2
7.2
7.2
7.1
7.2
7.3
7.1
7.2
7.1
5.9
7.1
7.2

4.5
4.7
4.6
4.7
4.7
4.5
4.5
4.6
4.7
4.7
4.7
4.8
4.7
4.8
4.8
4.6
4.8
4.8
3.9
4.7
4.8

2.6
2.6
2.6
2.7
2.7
2.6
2.6
2.7
2.7
2.7
2.7
2.7
2.7
2.7
2.8
2.6
2.7
2.6
2.5
2.6
2.7

4.67
3.00
1.94
1.74
1.57
1.56
1.50
1.23
1.21
1.20
1.02
0.96
0.92
0.84
0.73
0.63
0.53
0.55
0.28
0.26
0.19

TABLE IX. Best-fit values and the 90% CL upper limits on the DSNB rates (in events · year−1 ) for the theoreticalmodels for phases SK-I to IV as well as for the combined analysis. Here the upper limits are given forEν > 17.3 MeV. For the Kresse þ 21 models, the “High,” “Fid,” and “Low” predictions correspond to the “W20BH2.7-α2.0,” “W18-BH2.7-α2.0,” and “W18-BH2.7-α2.0-He33” models from Ref. [12], respectively.
Best-fit
Model
Totani þ 95 ConstantKaplinghat þ 00 HMA (max)Horiuchi þ 09 6 MeV, max
Ando þ 03 (updated 05)
Kresse þ 21 (High, NO)
Galais þ 09 (NO)Galais þ 09 (IO)
Horiuchi þ 18 ξ2.5 ¼ 0.1
Kresse þ 21 (High, IO)
Kresse þ 21 (Fid, NO)
Kresse þ 21 (Fid, IO)
Kresse þ 21 (Low, NO)
Tabrizi þ 21 (NO)

90% CL limit

SK4

All

SK1

SK2

SK3

SK4

All

Pred.

4.5þ2.5−2.2
4.4þ2.4−2.2
4.6þ2.5−2.3
4.7þ2.5−2.4
4.5þ2.5−2.3
4.4þ2.4−2.2
4.4þ2.4−2.2
4.8þ2.6−2.4
4.6þ2.5−2.3
4.5þ2.5−2.2
4.5þ2.5−2.2
4.5þ2.5−2.2
4.6þ2.5−2.3

2.3þ1.6−1.5
2.2þ1.5−1.4
2.4þ1.7−1.6
2.4þ1.6−1.6
2.4þ1.6−1.5
2.2þ1.6−1.4
2.2þ1.6−1.4
2.5þ1.7−1.6
2.4þ1.6−1.5
2.3þ1.5−1.5
2.3þ1.5−1.5
2.2þ1.6−1.4
2.4þ1.6
−1.5

4.2
3.9
4.4
4.2
4.1
4.0
4.0
4.5
4.1
4.0
4.0
3.9
4.2

11.2
11.2
11.7
11.8
11.7
11.0
11.0
12.1
11.7
11.7
11.6
11.6
11.8

12.2
11.7
12.5
12.4
12.2
11.9
11.9
12.8
12.2
11.9
11.9
11.8
12.3

7.9
7.7
8.2
8.2
8.1
7.8
7.7
8.4
8.1
7.9
7.9
7.9
8.2

4.6
4.4
4.8
4.7
4.6
4.5
4.5
4.9
4.6
4.5
4.5
4.5
4.7

8.35
5.09
3.54
3.09
2.76
2.74
2.62
2.29
2.14
2.06
1.75
1.65
1.64

(Table continued)

122002-31

K. ABE et al.

PHYS. REV. D 104, 122002 (2021)

TABLE IX. (Continued)
Best-fit
Model
Kresse þ 21 (Low, IO)
Lunardini09 Failed SN
Hartmann þ 97 CE
Nakazato þ 15 (max, IO)
Horiuchi þ 18 ξ2.5 ¼ 0.5
Horiuchi þ 21
Malaney97 CGINakazato þ 15 (min, NO)

90% CL limit

SK4

All

SK1

SK2

SK3

SK4

All

Pred.

4.5þ2.5−2.2
5.0þ2.7−2.5
4.4þ2.4−2.2
5.0þ2.7−2.5
4.2þ2.3−2.1
3.7þ2.2−2.1
4.3þ2.4−2.1
4.7þ2.5−2.4

2.2þ1.6−1.4
2.6þ1.7−1.7
2.2þ1.6−1.4
2.5þ1.8−1.6
2.1þ1.5−1.4
2.0þ1.6−1.6
2.2þ1.5−1.4
2.4þ1.6
−1.6

4.0
4.5
4.0
4.5
3.5
6.0
3.8
4.0

11.6
12.7
11.2
12.1
11.4
7.5
11.2
11.8

11.9
13.1
11.9
13.0
10.9
10.2
11.3
12.1

7.9
8.8
7.7
8.7
7.4
6.8
7.6
8.2

4.5
5.1
4.5
5.0
4.1
4.2
4.3
4.6

1.43
1.36
1.09
0.99
0.87
0.49
0.44
0.34

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 thestability of the observed and expected limits on the
DSNB flux over all models. This stability results fromthe 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 Sec. III, as inRef. [22]. The associated two-dimensional 90% CL limiton the neutrino emission temperature and on the DSNB rate
for Eν > 17.3 MeV is shown in Fig. 30 for the combined
SK-I,II,III,IV analysis. This figure also shows the expected
rates for different effective neutrino emission temperatures.
Here, the band accounts for the current uncertainty on the

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

122002-32

DIFFUSE SUPERNOVA NEUTRINO BACKGROUND SEARCH AT …

FIG. 30. The 90% CL excluded region for the DSNB rate forEν > 17.3 MeV and the effective neutrino temperature, for theblackbody models described in Sec. III (red surface). To guide theeye, we also show the one-dimensional exclusion limit forindividual neutrino temperatures. The blue band shows thepredictions for blackbody models, for the low, fiducial, and highstar formation rates used in Ref. [16].

cosmic star formation rate, using the low, fiducial, and highestimates from Ref. [16]. Since we used two-dimensional
exclusion contours instead of 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 complicated by the high SK analysis
threshold, we also show 90% CL exclusion contours forthe total ν̄e energy and the effective 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 this supernova is particularly low, future
experiments should be able to probe at least the associated
IMB region.
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 us toconsiderably 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 overthe DSNB up to about 15.5 MeV. Reactor antineutrinos,
in particular, provide an irreducible background up to
9.5 MeV. Up to about 11.5 MeV, 9Li decays might then

PHYS. REV. D 104, 122002 (2021)

FIG. 31. The 90% CL excluded region for the total ν̄e energy
emitted during a supernova and the effective neutrino temperature
for the blackbody models described in Sec. III (red surface). To
guide the eye, we also show the one-dimensional exclusion limit
for individual neutrino temperatures. We also show the 90% CL
contours corresponding to the range of parameters allowed forSN1987A.

dominate, and spallation reduction will remain a key
component of future analyses. Locating muon-induced
showers using neutrons will become especially powerfulat SK-Gd. After eliminating spallation, atmospheric NCQE
interactions will dominate over a wide energy range and willcall for new dedicated reduction techniques. At SK-Gd,
notably, the high visibility of the neutron capture signal willallow us to locate neutron capture without using the positronvertex, which could allow us 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 us to
improve the characterization of the NCQE rate, spectral
shape, and neutron multiplicity that are currently associatedwith large systematic uncertainties [48]. In the future, further
studies using the upcoming Intermediate Water Cherenkov
Detector will allow us 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].
IX. CONCLUSION
In this paper, we presented two analyses taking
advantage of the longest data-taking phase of SK, theimproved data reduction, and better background estimates.

122002-33

K. ABE et al.

PHYS. REV. D 104, 122002 (2021)

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 themodel-independent analysis, we achieved the world’s
tightest upper limit on the extraterrestrial electron antineutrino flux. For the latter region, we notably extended thespectral analysis developed for SK-I,II,III [22], and combined the results of the first four SK phases to reach an
exposure of 22.5 × 5823 kton · days. This exposure yields
SK a world-leading sensitivity to the DSNB flux at
90% CL, comparable to the fluxes of several realistic
models. A 1.5σ excess has been observed across all themodels probed and the final 90% CL 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 achieved at HK
will allow us to explore this parameter space, setting
meaningful constraints on astrophysical observables, withthe realistic prospect of a ground-breaking discovery.
ACKNOWLEDGMENTS
We gratefully acknowledge the cooperation of the
Kamioka Mining and Smelting Company. The SuperKamiokande experiment has been built and operated fromfunding by the Japanese Ministry of Education, Culture,
Sports, Science and Technology, the U.S. Department ofEnergy, and the U.S. National Science Foundation. Some ofus have been supported by funds from the National Research
Foundation of Korea NRF-2009-0083526 (KNRC) fundedby the Ministry of Science, ICT, and Future Planning and theMinistry of Education (No. 2018R1D1A3B07050696,
No. 2018R1D1A1B07049158), the Japan Society for the
Promotion of Science, the National Natural Science
Foundation of China under Grant No. 11620101004, theSpanish Ministry of Science, Universities and Innovation(Grant No. PGC2018-099388-B-I00), the Natural Sciencesand Engineering Research Council (NSERC) of Canada, theScinet and Westgrid consortia of Compute Canada, theNational Science Centre, Poland (2015/18/E/ST2/00758),
the Science and Technology Facilities Council (STFC) and
GridPPP, UK, the European Unions Horizon 2020 Researchand Innovation Programme under the Marie SklodowskaCurie Grant agreement No. 754496, H2020-MSCA-RISE2018 JENNIFER2 Grant agreement No. 822070, and H2020MSCA-RISE-2019 SK2HK Grant agreement No. 872549.
APPENDIX A: SPALLATION CUTTABLE 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)ll (cm)lt (cm)

0.1
1200
1200

1
800
800

30
383
219

60
548
268

60
603
379

60
712
490

60
766
548

vertex, ll is defined as the distance to this origin along themuon track while lt is the transverse distance to the track.
For a given cloud multiplicity, a DSNB candidate observed
less than dt after the associated muon, and within the ellipse
defined by ll and lt , is classified as spallation background.
Conditions for the neutron cloud cut are summarized inTable X for different multiplicities of the cloud.
2. Rectangular cuts
As mentioned in Sec. 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 eventswith 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 lt bins to account for the possible
correlations, 0–0.05, 0.05–0.5, and 0.5–30 sec for dt(determined by the isotope half-lives in Fig. 3), and
0–300, 300–1000, 1000–5000 cm for lt. Figure 32 shows
PDFs for each spallation variable from single through-going
muons. Here the results from the most spallation-rich region(0 < dt < 0.05 sec and 0 < lt < 300 cm) are shown. We
calculate the spallation log-likelihood ratio with the procedure described in Sec. V B. An example distribution isshown for single through-going muons in Fig. 33. The cut
criteria are tuned depending on the reconstructed energy as
well as the dt and lt bins. We use the common binning for
dt, 0 − 0.05, 0.05 − 0.5, and 0.5–30 sec, independent of themuon type. In contrast, for lt, we separate bins only for mostimpactful muon types, that is, single through-going andmultiple 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.

1. Neutron cloud cut
For each cloud, we define a specific coordinate system
whose origin is the projection of the cloud barycenter on
the muon track. For a given DSNB candidate reconstructed

4. Spallation cut performance estimation
In this section we present our methodology to estimate
the impact of spallation cuts on nonspallation “signal”

122002-34

DIFFUSE SUPERNOVA NEUTRINO BACKGROUND SEARCH AT …

PHYS. REV. D 104, 122002 (2021)

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 lt , at least one should be fulfilled.
Muon typeAll
MisfitWell-fitted single through-going
Stopping
Poorly-fitted stopping
All multiple

Energy range

dt requirement

lt requirement

Erec < 23.5 MeV
Erec < 23.5 MeV
15.5 < Erec < 19.5 MeV
15.5 < Erec < 17.5 MeV
15.5 < Erec < 19.5 MeV
15.5 < Erec < 17.5 MeV
15.5 < Erec < 19.5 MeV

dt > 0.001 sec
dt > 0.1 sec
dt > 1.5 sec
dt > 7 sec
dt > 0.05 sec
dt > 6 sec
dt > 0.05 sec

lt > 400 cm
lt > 150 cm




TABLE XII. Summary on the rectangular spallation cuts for events with tagged neutrons different from one. Whena given cut has requirements on both dt and lt , at least one should be fulfilled. Here, tighter rectangular cuts areapplied to the muon categories generating most of the spallation—single through-going and multiple muons—toremove contributions from short-lived isotopes.
Muon typeAll
MisfitWell-fitted single through-going
Stopping
All multiple
Well-fitted multiple

Energy range

dt requirement

lt requirement

Erec < 23.5 MeV
Erec < 23.5 MeV
15.5 < Erec < 19.5 MeV
15.5 < Erec < 17.5 MeV
17.5 < Erec < 19.5 MeV
15.5 < Erec < 19.5 MeV
15.5 < Erec < 19.5 MeV
15.5 < Erec < 19.5 MeV
15.5 < Erec < 19.5 MeV

dt > 0.001 sec
dt > 0.1 sec
dt > 1.5 sec
dt > 7 sec
dt > 7 sec
dt > 0.05 sec
dt > 0.05 sec
dt > 0.1 sec
dt > 1 sec

lt > 400 cm
lt > 150 cm
lt > 100 cm



lt > 400 cm

events (DSNB, solar, reactor, and atmospheric neutrino
interactions), as well as on spallation backgrounds. In what
follows, we treat 9Li decays separately as they provide thedominant β þ n signature associated with spallation.
As explained in Sec. 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 thepostsample, as described in the previous section. The
associated efficiency (ϵrandom ) can thus be expressed as
ϵrandom ¼ N post;after =N before ;ðA1Þ

where N before and N post;after are the numbers of events in thepostsample before and after spallation cuts. We validate thisefficiency by evaluating the impact of spallation cuts on
events from solar neutrino elastic interactions, using the
opening angle θsun defined in Sec. 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 wall and theCherenkov angle, which is defined in Sec. V C, and require
Erec < 15.5 MeV. The spallation remaining rate (ϵspall ) isthen estimated as the fraction of events remaining after

spallation cuts, after the expected number of remaining
atmospheric neutrino events from Monte-Carlo simulation
has been subtracted.
The procedure outlined above allows us 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 torenormalizing the spallation likelihood from Eq. (3) by
the isotope’s exponential decay time distribution. Then,
we apply the likelihood cuts described above to the pairs
with the modified dt and compute the fraction of remaining
events, ϵpre;iso . For each event, this fraction can be
expressed asðpairÞ

ϵpre;iso ¼ ϵiso;likeli × ϵrandom;iso ;ðpairÞ

ðA2Þ

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

122002-35

K. ABE et al.

PHYS. REV. D 104, 122002 (2021)
0.4

0.1

0.35
Spallation PDF

0.08

0.3

Random PDF
prob.

prob.

0.25
0.06

0.2
0.15

0.04

0.1
0.02
0.05
0
10−2

10−1

1
dt [sec]

0
0

10

500

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

0.16
−1

10

0.14
0.12

−2

10
prob.

prob.

0.1
0.08

−3

0.06

10

0.04
0.02

10−4

0
−5000 −4000 −3000 −2000 −1000

301000 2000 3000 4000 5000
ll [cm]

0

500

1000
1500
Qμ [p.e.]

2000

×10
2500

10−1

prob.

10−2

−3

10

−4

10

−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), lt (top right), ll (middle left), Qμ (middle right), and Qres (bottom), for thepre (blue) and post (red) samples. Results from the single through-going muons for the dt region of 0–0.05 sec and the lt regionof 0–300 cm are shown here.

muons to pass these cuts. This probability can be computed
using a postsample analogous to the one described inSec. V B 2 and, since we reassigned dt for all pairs in each
event, it will be artificially lower than the signal efficiency
ϵrandom . The probability for a given isotope decay to pass thelikelihood cuts (ϵiso;likeli ) is thus given by
ϵiso;likeli ¼ ϵpre;iso =ϵrandom;iso × ϵrandom :ðA3Þ

The fraction of remaining events after both preselection andlikelihood cuts (ϵiso;spall ) can then be estimated as

ϵiso;spall ¼ ϵpre;iso =ϵrandom;iso × ϵrandom × ϵpresel ;where ϵpresel ¼

ϵspall:ϵ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 andlikelihood cuts.

122002-36

DIFFUSE SUPERNOVA NEUTRINO BACKGROUND SEARCH AT …
underflow bin

overflow bin

104

Pre sample

EventsPost sample

3

10

102

-30

-20

-10

0
Lspall1020

30FIG. 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.

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 9Li background. Above
15.5 MeV, we will estimate the remaining fractions of 8Li
and 8B, two isotopes that are particularly difficult to remove
due to their long half-lives. This estimate will allow us tomodel the spallation background spectrum for the spectral
analysis.
In order to estimate the remaining amounts of individualisotopes without relying on a simulation, we made two
major assumptions. First, that the time difference between
spallation events and their parent muon is not correlatedwith other spallation observables. Second, that the impactof spallation cuts above 15.5 MeV can be estimated using asample of lower energy events. Both assumptions would beverified if the lt , ll , Qμ , and Qres spallation observables
described in Sec. V B 2 were similarly 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 2-MeV bins from 7.5 to
15.5 MeV that will have different isotope compositions.
For a given spallation cut, we then compute the spallation
remaining rates for the four subsamples and take the largest
deviation from the average value as our systematic uncertainty. This uncertainty depends on the isotope considered
but typically lies around 50%.
APPENDIX B: 241Am=Be CALIBRATION STUDYFOR NEUTRON TAGGINGAs described in Sec. V D 3, in order to validate theefficiency of the neutron tagging procedure on real data, we

PHYS. REV. D 104, 122002 (2021)

applied neutron tagging to data collected in the presence ofan 241Am=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 241Am=Be source is
not fully modeled by the simulation software; notably, no
modeling of the BGO scintillator is available. The MonteCarlo simulation that we use to model 241Am=Be events ishence similar to the IBD simulation. Because of thislimitation, 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 aneutron, thus artificially reducing the observed preselection
efficiency of the neutron candidates (see Sec. V D 1). To
account for this, we treated uncertainty associated with thepreselection separately from that of the BDT selection,
following the procedure described in Ref. [76].
The preselection efficiency is based on the numberof hits in a given time window and is hence particularly sensitive to the PMT photon detection efficiency(quantum efficiency, QE). Hence, we determine the associated systematic uncertainty by determining the value of
QE for which the predicted N 10 distribution fits the

FIG. 34. N 10 observed and expected distributions for neutrons.
The expected distribution corresponds to the PMT photon
detection efficiency for which the N 10 fits the observations best.
The observed distribution has been obtained by subtracting the
N 10 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
N 10 ¼ 6 is due to harsher continuous dark noise reduction cuts,
as mentioned in Sec. V D 1.

122002-37

K. ABE et al.

PHYS. REV. D 104, 122002 (2021)

TABLE XIII. Observed (“Data”) and predicted (“MC”) relative
efficiencies (in %) for different BDT cuts and for the 241Am=Be
calibration runs considered in this study. Data was taken in 2009
and 2016, at three different locations in the tank: in the center ofthe detector, near the barrel wall (Y12), and the top wall (Z15).
Year

Data

MC

ðData − MCÞ=Data %

0.40

2009
2016
2009
2016
2009
2016

69.4
80.0
74.6
83.8
87.6
80.2

76.2
72.3
77.4
77.6
77.1
74.1

−9.8
9.5
−3.7
7.2
12.0
−7.7

2009
2016
2009
2016
2009
2016

50.0
56.1
56.4
60.8
65.8
57.0

55.3
52.8
59.6
59.2
58.3
55.4

−10.6
5.8
−5.7
2.6
11.4
2.8

2009
2016
2009
2016
2009
2016

33.6
34.9
36.9
38.2
38.8
31.9

35.5
33.2
37.8
37.1
34.5
32.2

−5.8
4.8
−2.4
2.8
10.5
−0.1

Center
Y12
Z15

0.90

Center
Y12
Z15

0.99

Center
Y12
Z15

102
Events

BDT cut Position

CCQE-likeCC1π
CC other
NCQE
NC1π
NC other

3

10

10

1
10−1

0

1

2

3

4

5
6
Ntagged-n

7

8

9

10

FIG. 35. Predicted neutron multiplicity for atmospheric neutrino backgrounds. Here we use the atmospheric neutrino MonteCarlo simulation in the 7.5–79.5 MeV range, after the noise
reduction cuts.

APPENDIX C: NEUTRONS PRODUCED FROMATMOSPHERIC NEUTRINO INTERACTIONS

observations best. The N 10 predicted and observed distributions are shown in Fig. 34 for the best-fit QE, which isfound to be 2.1% larger than the QE used in the IBD
simulation from Sec. III B. The preselection efficiency
associated with this best-fit QE is 3.0% larger than theefficiency 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 fordifferent values of the BDT discriminant. We determine
these efficiencies by fitting the timing distributions ofneutron candidates by an exponential plus a constant.
Results for different BDT cuts are shown in Table XIII
for the different 241Am=Be datasets considered here. The
maximum discrepancy between data and Monte-Carlo lies
around 12% irrespective of the BDT cut.
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 includes possible discrepancies in neutron
kinetic energy distributions, like the ones between 9Li
decays and IBD events.

The neutron multiplicities predicted by the Monte-Carlosimulation for atmospheric neutrinos are shown in Fig. 35.
Note that a significant number of atmospheric neutrino
interactions produce one or more neutrons. However, theseneutrons are more energetic than for the IBDs of DSNB
neutrinos and sometimes travel further away from theirproduction vertex. Since the current neutron tagging
algorithm requires neutrons to be close to their associatedpositron vertex, it can also be used to reduce neutronproducing atmospheric backgrounds. A detailed study ofneutrons associated with atmospheric neutrino backgrounds has been performed in Ref. [76].
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 Sec. 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 inSec. VII, atmospheric background fluxes are treated as
nuisance parameters. In Table XIV, we compare thepredicted 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 Sec. VII, atmospheric events from, e.g., pionproducing NC interactions could be categorized as either
NCQE or μ=π, we quote ranges rather than numbers for

122002-38

DIFFUSE SUPERNOVA NEUTRINO BACKGROUND SEARCH AT …

PHYS. REV. D 104, 122002 (2021)

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

these two categories. For decay electrons and νe CCinteractions, the predicted and observed numbers of events
are within, at most, 23% of each other.
2. Modeling of spallation backgrounds
To model the spallation background PDF we vary theisotope composition of the background as described in
Eq. (3) and take the nominal spectrum as the average
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
Observed

Decay-e

νe CC

NCQE

μ=π

218.5
169.5

121.4
109.9

81.9–138.4
88.2

10.6–67.2
27.9

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.

122002-39

K. ABE et al.

PHYS. REV. D 104, 122002 (2021)ðErec þ 0.511 MeVÞα−βbetween the two extreme slopes. We parametrize thisaverage by the following function:ðErec þ 0.511 MeVÞα
Sspall ðErec Þ ¼ N exp −; ðD1Þ
β

S0spall ðErec Þ ¼ N 0 exp

where N is a normalization factor and α, β are determinedby a fit and depend on the SK phase considered. Then,
as described in Sec. VII C 2, we both vary the isotope
composition of the spallation spectrum and distort it toaccount for energy scale and resolution uncertainties. We
then fit the resulting extreme spectra by the function from
Eq. (D1) times a third order polynomial:and thus parametrize the 1σ systematic error distortions.
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 thespallation spectrum while avoiding unphysical results, we
impose S0spall ðErec Þ ¼ 0 for energies larger than the minimal
Erec for which the distorted PDF is zero.

[1] M. Fukugita et al., Physics and Astrophysics of Neutrinos(Springer, New York, 1994), ISBN 978-4-431-67029-2.
[2] K. Kotake, K. Sato, and K. Takahashi, Rep. Prog. Phys. 69,
971 (2006).
[3] K. Nakazato et al., Astrophys. J. Suppl. Ser. 205, 2 (2013).
[4] G. Tammann et al., Astrophys. J. Suppl. Ser. 92, 487 (1994).
[5] R. Diehl et al., Nature (London) 439, 45 (2006).
[6] K. Nakazato, E. Mochida, Y. Niino, and H. Suzuki,
Astrophys. J. 804, 75 (2015).
[7] R. A. Malaney, Astropart. Phys. 7, 125 (1997).
[8] T. Totani and K. Sato, Astropart. Phys. 3, 367 (1995).
[9] M. Kaplinghat, G. Steigman, and T. P. Walker, Phys. Rev. D
62, 043001 (2000).
[10] C. Lunardini, Phys. Rev. Lett. 102, 231101 (2009).
[11] S. Horiuchi, K. Sumiyoshi, K. Nakamura, T. Fischer, A.
Summa, T. Takiwaki, H.-T. Janka, and K. Kotake, Mon.
Not. R. Astron. Soc. 475, 1363 (2018).
[12] D. Kresse, T. Ertl, and H.-T. Janka, Astrophys. J. 909, 169(2021).
[13] S. Horiuchi, T. Kinugawa, T. Takiwaki, K. Takahashi, andK. Kotake, Phys. Rev. D 103, 043003 (2021).
[14] Z. Tabrizi and S. Horiuchi, J. Cosmol. Astropart. Phys. 05(2021) 011.
[15] S. Galais, J. Kneller, C. Volpe, and J. Gava, Phys. Rev. D 81,
053002 (2010).
[16] S. Horiuchi, J. F. Beacom, and E. Dwek, Phys. Rev. D 79,
083013 (2009).
[17] S. Ando, K. Sato, and T. Totani, Astropart. Phys. 18, 307(2003).
[18] D. H. Hartmann and S. E. Woosley, Astropart. Phys. 7, 137(1997).
[19] S. I. Ando, Proceedings of the Next Generation of Nucleon
Decay and Neutrino Detectors, Aussois, Savoie, France,
2005, edited by J. Dumarchez (2005).
[20] D. Cokinos and E. Melkonian, Phys. Rev. C 15, 1636 (1977).
[21] M. Malek et al. (Super-Kamiokande Collaboration), Phys.
Rev. Lett. 90, 061101 (2003).
[22] K. Bays et al. (Super-Kamiokande Collaboration), Phys.
Rev. D 85, 052007 (2012).

[23] H. Zhang et al. (Super-Kamiokande Collaboration), Astropart. Phys. 60, 41 (2015).
[24] A. Gando et al. (KamLAND Collaboration), Astrophys. J.
745, 193 (2012).
[25] S. Abe et al. (KamLAND Collaboration), arXiv:2108.08527.
[26] B. Aharmim et al. (SNO Collaboration), Astrophys. J. 653,
1545 (2006).
[27] B. Aharmim et al. (SNO Collaboration), Phys. Rev. D 102,062006 (2020).
[28] M. Agostini et al. (Borexino Collaboration), Astropart.
Phys. 125, 102509 (2021).
[29] S. Fukuda et al. (Super-Kamiokande Collaboration), Nucl.
Instrum. Methods Phys. Res., Sect. A 501, 418 (2003).
[30] H. Kume et al., Nucl. Instrum. Methods Phys. Res., Sect. A
205, 443 (1983).
[31] K. Abe et al. (Super-Kamiokande Collaboration), Nucl.
Instrum. Methods Phys. Res., Sect. A 737, 253 (2014).
[32] M. Nakahata et al. (Super-Kamiokande Collaboration),
Nucl. Instrum. Methods Phys. Res., Sect. A 421, 113(1999).
[33] E. Blaufuss et al. (Super-Kamiokande Collaboration), Nucl.
Instrum. Methods Phys. Res., Sect. A 458, 638 (2001).
[34] H. Nishino, K. Awai, Y. Hayato, S. Nakayama, K. Okumura,
M. Shiozawa, A. Takeda, K. Ishikawa, A. Minegishi, and Y.
Arai, Nucl. Instrum. Methods Phys. Res., Sect. A 610, 710(2009).
[35] H. Watanabe, H. Zhang, K. Abe, Y. Hayato, T. Iida, M.
Ikeda, J. Kameda, K. Kobayashi, Y. Koshio, and M. Miura(Super-Kamiokande Collaboration), Astropart. Phys. 31,
320 (2009).
[36] J. Hosaka et al. (Super-Kamiokande Collaboration), Phys.
Rev. D 73, 112001 (2006).
[37] J. Cravens et al. (Super-Kamiokande Collaboration), Phys.
Rev. D 78, 032002 (2008).
[38] K. Abe et al. (Super-Kamiokande Collaboration), Phys.
Rev. D 83, 052010 (2011).
[39] K. Abe et al. (Super-Kamiokande Collaboration), Phys.
Rev. D 94, 052010 (2016).

× ½1 þ P 3 ðErec Þ;

122002-40

ðD2ÞDIFFUSE SUPERNOVA NEUTRINO BACKGROUND SEARCH AT …
[40] M. Smy, Proc. Int. Cosmic Ray Conf. 5, 1279 (2008),
https://inspirehep.net/literature/1371616.
[41] E. E. Salpeter, Astrophys. J. 121, 161 (1955).
[42] M. T. Keil, G. G. Raffelt, and H.‐T. Janka, Astrophys. J. 590,
971 (2003).
[43] A. Strumia and F. Vissani, Phys. Lett. B 564, 42 (2003).
[44] R. Brun et al., CERN Report No. CERN-W5013, 1994.
[45] A. Ankowski, O. Benhar, T. Mori, R. Yamaguchi, and M.
Sakuda, Phys. Rev. Lett. 108, 052505 (2012).
[46] K. Abe et al. (T2K Collaboration), Phys. Rev. D 90, 072012(2014).
[47] L. Wan et al. (Super-Kamiokande Collaboration), Phys.
Rev. D 99, 032005 (2019).
[48] K. Abe et al. (T2K Collaboration), Phys. Rev. D 100,112009 (2019).
[49] M. Honda, T. Kajita, K. Kasahara, S. Midorikawa, and T.
Sanuki, Phys. Rev. D 75, 043006 (2007).
[50] M. Honda, T. Kajita, K. Kasahara, and S. Midorikawa, Phys.
Rev. D 83, 123001 (2011).
[51] http://www.icrr.u-tokyo.ac.jp/~mhonda/nflx2011/index.html (HKKM2011 Flux).
[52] Y. Hayato, Acta Phys. Pol. B 40, 2477 (2009), https://www.actaphys.uj.edu.pl/fulltext?series=Reg&vol=40&page=2477.
[53] K. Abe et al. (Super-Kamiokande Collaboration), arXiv:
2012.03807.
[54] M. Jiang et al. (Super-Kamiokande Collaboration), Prog.
Theor. Exp. Phys. 2019, 053F01 (2019).
[55] G. Battistoni, F. Cerutti, A. Fassò, A. Ferrari, S. Muraro, J.
Ranft, S. Roesler, and P. R. Sala, AIP Conf. Proc. 896, 31(2007).
[56] S. W. Li and J. Beacom, Phys. Rev. C 89, 045801 (2014).
[57] S. W. Li and J. Beacom, Phys. Rev. D 91, 105005 (2015).
[58] S. W. Li and J. Beacom, Phys. Rev. D 92, 105033 (2015).
[59] Y. Zhang et al. (Super-Kamiokande Collaboration), Phys.
Rev. D 93, 012004 (2016).
[60] https://github.com/Goldie643/SKReact (SKReact).
[61] M. Baldoncini, I. Callegari, G. Fiorentini, F. Mantovani, B.
Ricci, V. Strati, and G. Xhixha, Phys. Rev. D 91, 065002(2015).
[62] Technical Reports Series No. 428 International Atomic
Energy Agency, Vienna, 2005.
[63] F. P. An et al. (Daya Bay Collaboration), Chin. Phys. C 41,
013002 (2017).
[64] S. M. Locke et al. (Super-Kamiokande Collaboration),
arXiv:2112.00092.
[65] G. Carminati, in Proceedings of the 32nd International
Cosmic Ray Conference (2011), http://icrc2011.ihep.ac.cn/
paper/proc/v4.pdf.
[66] G. Carminati (Super-Kamiokande Collaboration), Phys.
Procedia 61, 666 (2015).
[67] M. Elnimr (Super-Kamiokande Collaboration), J. Phys.
Conf. Ser. 888, 012189 (2017).
[68] Z. Conner, Ph.D. Thesis, Boston University, 1997.
[69] S. Desai, Ph.D. Thesis, Boston University, 2004.

PHYS. REV. D 104, 122002 (2021)[70] K. Abe et al. (T2K Collaboration), Phys. Rev. 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, T. Hokama, M. Matsubara, M. Miwa, M.
Nakahata, T. Nakamura, H. Sekiya, Y. Takeuchi, S. Tasaka,
and R.A. Wendell, Nucl. Instrum. Methods Phys. Res., Sect.
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 be
published).
[77] K. Abe et al. (Super-Kamiokande Collaboration) (to be
published).
[78] W. A. Rolke, A. M. López, and J. Conrad, Nucl. Instrum.
Methods Phys. Res., Sect. A 551, 493 (2005).
[79] O. Benhar, A. Fabrocini, S. Fantoni, and I. Sick, Nucl. Phys.
A579, 493 (1994).
[80] O. Benhar, N. Farina, H. Nakamura, M. Sakuda, and R.
Seki, Phys. Rev. D 72, 053005 (2005).
[81] C. J. Horowitz and B. D. Serot, Nucl. Phys. A368, 503(1981).
[82] C. Maieron, M. C. Martínez, J. A. Caballero, and J. M.
Udías, Phys. Rev. C 68, 048501 (2003).
[83] J. A. Caballero, J. E. Amaro, M. B. Barbaro, T. W. Donnelly,
C. Maieron, and J. M. Udias, Phys. Rev. Lett. 95, 252502(2005).
[84] R. Gonzalez-Jimenez, J. A. Caballero, A. Meucci, C. Giusti,
M. B. Barbaro, M. V. Ivanov, and J. M. Udías, Phys. Rev. C
88, 025502 (2013).
[85] J. E. Amaro, M. B. Barbaro, J. A. Caballero, T. W. Donnelly,
A. Molinari, and I. Sick, Phys. Rev. C 71, 015501(2005).
[86] J. E. Amaro, M. B. Barbaro, J. A. Caballero, and T. W.
Donnelly, Phys. Rev. C 73, 035503 (2006).
[87] F. Capuzzi, C. Giusti, and F. D. Pacati, Nucl. Phys. A524,
681 (1991).
[88] A. Meucci, C. Giusti, and F. D. Pacati, Nucl. Phys. A739,
277 (2004).
[89] A. Meucci and C. Giusti, Phys. Rev. D 89, 057302 (2014).
[90] Y. Ashie et al. (Super-Kamiokande Collaboration), Phys.
Rev. D 71, 112005 (2005).
[91] C. Regis et al. (Super-Kamiokande Collaboration), Phys.
Rev. D 86, 012006 (2012).
[92] W. T. Winter, S. J. Freedman, K. E. Rehm, and J. P. Schiffer,Phys. Rev. C 73, 025503 (2006).
[93] B. Jegerlehner, F. Neubig, and G. Raffelt, Phys. Rev. D 54,
1194 (1996).
[94] M. Sisti (JUNO Collaboration), J. Phys. Conf. Ser. 1468,
012150 (2021).
[95] K. Møller, A. M. Suliga, I. Tamborra, and P. B. Denton,
J. Cosmol. Astropart. Phys. 05 (2018) 066.

122002-41