Super-Kamiokande Collaboration

# Diffuse Supernova Neutrino Background Search at Super-Kamiokande

K. Abe  Kamioka Observatory, Institute for Cosmic Ray Research, University of
Tokyo, Kamioka, Gifu 506-1205, Japan  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  C.
Bronner  Kamioka Observatory, Institute for Cosmic Ray Research, University of
Tokyo, Kamioka, Gifu 506-1205, Japan  Y. Hayato  Kamioka Observatory,
Institute for Cosmic Ray Research, University of Tokyo, Kamioka, Gifu
506-1205, Japan  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  K. Hiraide  M. Ikeda  S.
Imaizumi  Kamioka Observatory, Institute for Cosmic Ray Research, University
of Tokyo, Kamioka, Gifu 506-1205, Japan  J. Kameda  Kamioka Observatory,
Institute for Cosmic Ray Research, University of Tokyo, Kamioka, Gifu
506-1205, Japan  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  Y. Kanemura  Y. Kataoka
S. Miki  Kamioka Observatory, Institute for Cosmic Ray Research, University of
Tokyo, Kamioka, Gifu 506-1205, Japan  M. Miura  S. Moriyama  Kamioka
Observatory, Institute for Cosmic Ray Research, University of Tokyo, Kamioka,
Gifu 506-1205, Japan  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  Y. Nagao  Kamioka
Observatory, Institute for Cosmic Ray Research, University of Tokyo, Kamioka,
Gifu 506-1205, Japan  M. Nakahata  Kamioka Observatory, Institute for Cosmic
Ray Research, University of Tokyo, Kamioka, Gifu 506-1205, Japan  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  S. Nakayama  Kamioka Observatory, Institute
for Cosmic Ray Research, University of Tokyo, Kamioka, Gifu 506-1205, Japan
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  T. Okada  K. Okamoto  A. Orii  G. Pronost
Kamioka Observatory, Institute for Cosmic Ray Research, University of Tokyo,
Kamioka, Gifu 506-1205, Japan  H. Sekiya  M. Shiozawa  Kamioka Observatory,
Institute for Cosmic Ray Research, University of Tokyo, Kamioka, Gifu
506-1205, Japan  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  Y. Sonoda  Y. Suzuki
Kamioka Observatory, Institute for Cosmic Ray Research, University of Tokyo,
Kamioka, Gifu 506-1205, Japan  A. Takeda  Kamioka Observatory, Institute for
Cosmic Ray Research, University of Tokyo, Kamioka, Gifu 506-1205, Japan  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  Y. Takemoto  A. Takenaka  Kamioka Observatory,
Institute for Cosmic Ray Research, University of Tokyo, Kamioka, Gifu
506-1205, Japan  H. Tanaka  Kamioka Observatory, Institute for Cosmic Ray
Research, University of Tokyo, Kamioka, Gifu 506-1205, Japan  S. Watanabe
Kamioka Observatory, Institute for Cosmic Ray Research, University of Tokyo,
Kamioka, Gifu 506-1205, Japan  T. Yano  Kamioka Observatory, Institute for
Cosmic Ray Research, University of Tokyo, Kamioka, Gifu 506-1205, Japan  S.
Han  Research Center for Cosmic Neutrinos, Institute for Cosmic Ray Research,
University of Tokyo, Kashiwa, Chiba 277-8582, Japan  T. Kajita  Research
Center for Cosmic Neutrinos, Institute for Cosmic Ray Research, University of
Tokyo, Kashiwa, Chiba 277-8582, Japan  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  K.
Okumura  Research Center for Cosmic Neutrinos, Institute for Cosmic Ray
Research, University of Tokyo, Kashiwa, Chiba 277-8582, Japan  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  T. Tashiro  J. Xia  Research Center for Cosmic Neutrinos, Institute for
Cosmic Ray Research, University of Tokyo, Kashiwa, Chiba 277-8582, Japan  G.
D. Megias  Institute for Cosmic Ray Research, University of Tokyo, Kashiwa,
Chiba 277-8582, Japan  D. Bravo-Berguño  L. Labarga  Ll. Marti  B. Zaldivar
Department of Theoretical Physics, University Autonoma Madrid, 28049 Madrid,
Spain  B. W. Pointon  Department of Physics, British Columbia Institute of
Technology, Burnaby, BC, V5G 3H2, Canada  TRIUMF, 4004 Wesbrook Mall,
Vancouver, BC, V6T2A3, Canada  F. d. M. Blaszczyk  Department of Physics,
Boston University, Boston, MA 02215, USA  E. Kearns  Department of Physics,
Boston University, Boston, MA 02215, USA  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  J. L.
Raaf  Department of Physics, Boston University, Boston, MA 02215, USA  J. L.
Stone  Department of Physics, Boston University, Boston, MA 02215, USA  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  L. Wan  Department of Physics, Boston
University, Boston, MA 02215, USA  T. Wester  Department of Physics, Boston
University, Boston, MA 02215, USA  J. Bian  N. J. Griskevich  W. R. Kropp  S.
Locke  S. Mine  Department of Physics and Astronomy, University of California,
Irvine, Irvine, CA 92697-4575, USA  M. B. Smy  H. W. Sobel  Department of
Physics and Astronomy, University of California, Irvine, Irvine, CA
92697-4575, USA  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  V. Takhistov  Department
of Physics and Astronomy, University of California, Irvine, Irvine, CA
92697-4575, USA  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  J. Hill  Department of
Physics, California State University, Dominguez Hills, Carson, CA 90747, USA
J. Y. Kim  I. T. Lim  R. G. Park  Institute for Universe and Elementary
Particles, Chonnam National University, Gwangju 61186, Korea  B. Bodur
Department of Physics, Duke University, Durham NC 27708, USA  K. Scholberg  C.
W. Walter  Department of Physics, Duke University, Durham NC 27708, USA  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  S. Cao  Institute For Interdisciplinary
Research in Science and Education, ICISE, Quy Nhon, 55121, Vietnam  L. Bernard
A. Coffani  O. Drapier  S. El Hedri  Corresponding author:
<EMAIL_ADDRESS>A. Giampaolo  Corresponding author:
<EMAIL_ADDRESS>M. Gonin  Th. A. Mueller  P. Paganini  B. Quilain
Ecole Polytechnique, IN2P3-CNRS, Laboratoire Leprince-Ringuet, F-91120
Palaiseau, France  T. Ishizuka  Junior College, Fukuoka Institute of
Technology, Fukuoka, Fukuoka 811-0295, Japan  T. Nakamura  Department of
Physics, Gifu University, Gifu, Gifu 501-1193, Japan  J. S. Jang  GIST
College, Gwangju Institute of Science and Technology, Gwangju 500-712, Korea
J. G. Learned  Department of Physics and Astronomy, University of Hawaii,
Honolulu, HI 96822, USA  L. H. V. Anthony  D. Martin  M. Scott  A. A. Sztuc
Y. Uchida  Department of Physics, Imperial College London , London, SW7 2AZ,
United Kingdom  V. Berardi  M. G. Catanesi  E. Radicioni  Dipartimento
Interuniversitario di Fisica, INFN Sezione di Bari and Università e
Politecnico di Bari, I-70125, Bari, Italy  N. F. Calabria  L. N. Machado  G.
De Rosa  Dipartimento di Fisica, INFN Sezione di Napoli and Università di
Napoli, I-80126, Napoli, Italy  G. Collazuol  F. Iacob  M. Lamoureux  M.
Mattiazzi  N. Ospina  Dipartimento di Fisica, INFN Sezione di Padova and
Università di Padova, I-35131, Padova, Italy  L. Ludovici  INFN Sezione di
Roma and Università di Roma “La Sapienza”, I-00185, Roma, Italy  Y. Maekawa
Y. Nishimura  Department of Physics, Keio University, Yokohama, Kanagawa,
223-8522, Japan  M. Friend  T. Hasegawa  T. Ishida  T. Kobayashi  M. Jakkapu
T. Matsubara  T. Nakadaira  High Energy Accelerator Research Organization
(KEK), Tsukuba, Ibaraki 305-0801, Japan  K. Nakamura  High Energy Accelerator
Research Organization (KEK), Tsukuba, Ibaraki 305-0801, Japan  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  Y. Oyama  K. Sakashita  T. Sekiguchi  T. Tsukamoto  High Energy
Accelerator Research Organization (KEK), Tsukuba, Ibaraki 305-0801, Japan  Y.
Kotsar  Y. Nakano  H. Ozaki  T. Shiozawa  Department of Physics, Kobe
University, Kobe, Hyogo 657-8501, Japan  A. T. Suzuki  Department of Physics,
Kobe University, Kobe, Hyogo 657-8501, Japan  Y. Takeuchi  Department of
Physics, Kobe University, Kobe, Hyogo 657-8501, Japan  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  S. Yamamoto  Department of Physics, Kobe University, Kobe, Hyogo
657-8501, Japan  A. Ali  Y. Ashida  Corresponding author:
<EMAIL_ADDRESS>J. Feng  S. Hirota  T. Kikawa  M. Mori  Department
of Physics, Kyoto University, Kyoto, Kyoto 606-8502, Japan  T. Nakaya
Department of Physics, Kyoto University, Kyoto, Kyoto 606-8502, Japan  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  R. A. Wendell  Department of Physics, Kyoto
University, Kyoto, Kyoto 606-8502, Japan  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  K.
Yasutome  Department of Physics, Kyoto University, Kyoto, Kyoto 606-8502,
Japan  P. Fernandez  N. McCauley  P. Mehta  K. M. Tsui  Department of Physics,
University of Liverpool, Liverpool, L69 7ZE, United Kingdom  Y. Fukuda
Department of Physics, Miyagi University of Education, Sendai, Miyagi
980-0845, Japan  Y. Itow  Institute for Space-Earth Environmental Research,
Nagoya University, Nagoya, Aichi 464-8602, Japan  Kobayashi-Maskawa Institute
for the Origin of Particles and the Universe, Nagoya University, Nagoya, Aichi
464-8602, Japan  H. Menjo  T. Niwa  K. Sato  Institute for Space-Earth
Environmental Research, Nagoya University, Nagoya, Aichi 464-8602, Japan  M.
Tsukada  Institute for Space-Earth Environmental Research, Nagoya University,
Nagoya, Aichi 464-8602, Japan  J. Lagda  S. M. Lakshmi  P. Mijakowski  J.
Zalipska  National Centre For Nuclear Research, 02-093 Warsaw, Poland  J.
Jiang  C. K. Jung  C. Vilela  M. J. Wilking  C. Yanagisawa  also at BMCC/CUNY,
Science Department, New York, New York, 1007, USA.  Department of Physics and
Astronomy, State University of New York at Stony Brook, NY 11794-3800, USA  K.
Hagiwara  M. Harada  T. Horai  H. Ishino  S. Ito  H. Kitagawa  Department of
Physics, Okayama University, Okayama, Okayama 700-8530, Japan  Y. Koshio
Department of Physics, Okayama University, Okayama, Okayama 700-8530, Japan
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  W. Ma  N. Piplani  S. Sakai  Department of
Physics, Okayama University, Okayama, Okayama 700-8530, Japan  G. Barr  D.
Barrow  Department of Physics, Oxford University, Oxford, OX1 3PU, United
Kingdom  L. Cook  Department of Physics, Oxford University, Oxford, OX1 3PU,
United Kingdom  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  A. Goldsack  Department
of Physics, Oxford University, Oxford, OX1 3PU, United Kingdom  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  S. Samani  Department of Physics, Oxford
University, Oxford, OX1 3PU, United Kingdom  D. Wark  Department of Physics,
Oxford University, Oxford, OX1 3PU, United Kingdom  STFC, Rutherford Appleton
Laboratory, Harwell Oxford, and Daresbury Laboratory, Warrington, OX11 0QX,
United Kingdom  F. Nova  Rutherford Appleton Laboratory, Harwell, Oxford, OX11
0QX, UK  T. Boschi  F. Di Lodovico  J. Gao  J. Migenda  M. Taani  S. Zsoldos
Department of Physics, King’s College London, London, WC2R 2LS, UK  J. Y. Yang
Department of Physics, Seoul National University, Seoul 151-742, Korea  S. J.
Jenkins  M. Malek  J. M. McElwee  O. Stone  M. D. Thiesse  L. F. Thompson
Department of Physics and Astronomy, University of Sheffield, S3 7RH,
Sheffield, United Kingdom  H. Okazawa  Department of Informatics in Social
Welfare, Shizuoka University of Welfare, Yaizu, Shizuoka, 425-8611, Japan  S.
B. Kim  J. W. Seo  I. Yu  Department of Physics, Sungkyunkwan University,
Suwon 440-746, Korea  K. Nishijima  Department of Physics, Tokai University,
Hiratsuka, Kanagawa 259-1292, Japan  M. Koshiba  Deceased.  The University of
Tokyo, Bunkyo, Tokyo 113-0033, Japan  K. Iwamoto  Department of Physics,
University of Tokyo, Bunkyo, Tokyo 113-0033, Japan  K. Nakagiri  Department of
Physics, University of Tokyo, Bunkyo, Tokyo 113-0033, Japan  Y. Nakajima
Department of Physics, University of Tokyo, Bunkyo, Tokyo 113-0033, Japan
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  N. Ogawa  Department of Physics, University of
Tokyo, Bunkyo, Tokyo 113-0033, Japan  M. Yokoyama  Department of Physics,
University of Tokyo, Bunkyo, Tokyo 113-0033, Japan  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  K. Martens  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  M. R. Vagins  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  Department of Physics and Astronomy,
University of California, Irvine, Irvine, CA 92697-4575, USA  M. Kuze  S.
Izumiyama  T. Yoshida  Department of Physics,Tokyo Institute of Technology,
Meguro, Tokyo 152-8551, Japan  M. Inomoto  M. Ishitsuka  H. Ito  T. Kinoshita
R. Matsumoto  K. Ohta  M. Shinoki  T. Suganuma  Department of Physics, Faculty
of Science and Technology, Tokyo University of Science, Noda, Chiba 278-8510,
Japan  A. K. Ichikawa  K. Nakamura  Department of Physics, Faculty of Science,
Tohoku University, Sendai, Miyagi, 980-8578, Japan  J. F. Martin  H. A. Tanaka
T. Towstego  Department of Physics, University of Toronto, ON, M5S 1A7, Canada
R. Akutsu  TRIUMF, 4004 Wesbrook Mall, Vancouver, BC, V6T2A3, Canada  V.
Gousy-Leblanc  also at University of Victoria, Department of Physics and
Astronomy, PO Box 1700 STN CSC, Victoria, BC V8W 2Y2, Canada.  M. Hartz
TRIUMF, 4004 Wesbrook Mall, Vancouver, BC, V6T2A3, Canada  A. Konaka  also at
University of Victoria, Department of Physics and Astronomy, PO Box 1700 STN
CSC, Victoria, BC V8W 2Y2, Canada.  P. de Perio  N. W. Prouse  TRIUMF, 4004
Wesbrook Mall, Vancouver, BC, V6T2A3, Canada  S. Chen  B. D. Xu  Y. Zhang
Department of Engineering Physics, Tsinghua University, Beijing, 100084, China
M. Posiadala-Zezula  Faculty of Physics, University of Warsaw, Warsaw, 02-093,
Poland  D. Hadley  M. O’Flaherty  B. Richards  Department of Physics,
University of Warwick, Coventry, CV4 7AL, UK  B. Jamieson  J. Walker
Department of Physics, University of Winnipeg, MB R3J 3L8, Canada  A. Minamino
K. Okamoto  G. Pintaudi  S. Sano  R. Sasaki  Department of Physics, Yokohama
National University, Yokohama, Kanagawa, 240-8501, Japan

###### Abstract

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

##  I Introduction

###  I.1 Diffuse supernova neutrino background

Core-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 core of the collapsing star.
Information about this core could be accessed by detecting neutrinos emitted
by supernova bursts, whose luminosity and energy spectra closely track the
different steps of the CCSN mechanism in a neutrino-heating scenario (see
Refs. Fukugita  _et al._ (1994); Kotake  _et al._ (2006); Nakazato  _et al._
(2013), for example). Existing neutrino experiments, however, are mostly
sensitive to supernova bursts occurring in our galaxy and its immediate
surroundings, that are extremely rare (a few times per century in our galaxy
Tammann  _et al._ (1994); Diehl  _et al._ (2006)). Learning about the
aggregate properties of supernovae in the Universe hence requires observing
the accumulation of neutrinos from all distant supernovae. The integrated flux
of these neutrinos forms the diffuse supernova neutrino background, or DSNB
(DSNB neutrinos are also referred to as supernova relic neutrinos (SRNs) in
various publications).

The DSNB is composed of neutrinos of all flavors whose energies have been
redshifted when propagating to Earth. Its spectrum therefore contains unique
information not only on the supernova neutrino emission process but also on
the star formation and Universe expansion history. Spectra predicted from
various models are shown in Fig. 1. Their overall normalization is mostly
determined by the supernova rate, related to the cosmic star formation rate.
While this redshift-dependent rate also impacts the DSNB spectral shape, the
latter 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 of neutron
stars, the fraction of black-hole-forming supernovae, the revival time of the
shock wave, etc. Possible combinations of these factors were systematically
studied in the Nakazato$+$15 model Nakazato  _et al._ (2015), the maximal and
minimal fluxes of which are shown in the figure. The minimal prediction from
this model and the Malaney97 gas infall model Malaney (1997) give the lower
bounds on the DSNB flux. Conversely the Totani$+$95 model Totani and Sato
(1995) can be considered an upper bound on the expected DSNB flux, on par with
the most optimistic predictions of the Kaplinghat$+$00 model Kaplinghat  _et
al._ (2000), also shown in the figure. Between these bounds, the DSNB flux can
vary over one order of magnitude and its spectral shape bears the imprint of a
wide range of physical effects. Black-hole-forming supernovae notably make the
DSNB spectrum harder, as can be seen with the Lunardini09 model Lunardini
(2009), that assumes that $17\%$ of core-collapse supernovae lead to black-
hole formation. Similarly, the Horiuchi$+$18 model Horiuchi  _et al._ (2018)
incorporates a fraction of black-hole-forming supernovae determined using the
stars’ compactness. Aside from accounting for black-hole formation, the
Kresse$+$21 model Kresse  _et al._ (2021), 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 Horiuchi  _et al._ (2021)) discusses
the impact of interactions in binary systems, such as mergers, and mass
transfer, on the DSNB flux. Detecting the DSNB would hence provide valuable
insights into a wide array of physical processes.

Figure 1: DSNB $\bar{\nu}_{e}$ flux predictions from various theoretical
models (Horiuchi$+$21 Horiuchi  _et al._ (2021), Tabrizi$+$21 Tabrizi and
Horiuchi (2021), Kresse$+$21 Kresse  _et al._ (2021), Horiuchi$+$18 Horiuchi
_et al._ (2018), Nakazato$+$15 Nakazato  _et al._ (2015), Galais$+$10 Galais
_et al._ (2010), Horiuchi$+$09 Horiuchi  _et al._ (2009), Lunardini09
Lunardini (2009), Ando$+$09 Ando  _et al._ (2003), Kaplinghat$+$00 Kaplinghat
_et al._ (2000), Malaney97 Malaney (1997), Hartmann$+$97 Hartmann and Woosley
(1997), and Totani$+$95 Totani and Sato (1995)). Refer to each publication for
the detailed descriptions of model. In the legend, “NH” and “IH” 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 bib (2005). The corresponding flux is larger
by a factor of 2.56 than in the original publication Ando  _et al._ (2003).

###  I.2 Experimental searches

While a significant fraction of DSNB neutrinos are expected to have energies
lower than 10 MeV —as shown in Fig. 1—, the $\mathcal{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 $\mathcal{O}(10)$
MeV. At these energies, the dominant detection channel in most experiments is
the inverse beta decay (IBD) of electron antineutrinos
($\bar{\nu}_{e}+p\rightarrow 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 a neutron (Fig. 2). Here the neutron does not emit
Cherenkov or scintillation light directly but its capture on a proton leads to
emission of a $2.2$ MeV photon. In pure water, the characteristic timescale
for neutron moderation and capture is of $204.8\pm 0.4~{}\mu$s Cokinos and
Melkonian (1977). Pairing the “prompt” positron signal with the “delayed”
photon emission from neutron capture is hence a key component of many
analyses.

Figure 2: Schematic illustration of an IBD process and the subsequent neutron
capture on another proton. The characteristic neutron capture time in water is
$\tau=204.8\pm 0.4~{}\mu$s Cokinos and Melkonian (1977).

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) Malek  _et al._ (2003); Bays  _et al._ (2012); Zhang  _et al._ (2015),
KamLAND Gando  _et al._ (2012), SNO Aharmim  _et al._ (2006, 2020) and
Borexino Agostini  _et al._ (2021). Note that, unlike other searches, the SNO
analysis is sensitive to electron neutrinos and, due to the irreducible solar
neutrino background, its effective energy threshold lies around the hep solar
neutrino flux endpoint, around 19 MeV. Among all past analyses, the SK and
KamLAND experiments placed the most stringent upper limits on the DSNB
$\bar{\nu}_{e}$ flux. At SK the first DSNB search was carried out in 2003
using a $22.5\times 1496$-kton$\cdot$day data set Malek  _et al._ (2003).
Using spectral shape fitting for signal and atmospheric neutrino backgrounds,
it placed an upper limit on the DSNB flux for a wide variety of models in the
19.3$-$83.3 MeV neutrino energy range. This analysis already allowed to
disfavor the most optimistic DSNB predictions, in particular the Totani$+$95
model Totani and Sato (1995), and constrain the parameter space of the
Kaplinghat$+$00 model Kaplinghat  _et al._ (2000). In 2012, an improved
analysis was performed at SK, using a $22.5\times 2853$-kton$\cdot$day
exposure, a lower neutrino energy threshold of $17.3$ MeV, and new event
selection cuts allowing a 50% increase in signal efficiency Bays  _et al._
(2012). Below 17.3 MeV, large backgrounds from cosmic muon spallation and
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 2012 Gando  _et al._
(2012) probed neutrino energies ranging from $8.3$ to $31.8$ MeV by
investigating coincident positron and neutron capture signals. While neutron
identification in a liquid scintillator detector such as KamLAND is
significantly easier than in pure water, KamLAND’s small fiducial volume only
allowed for an exposure of $4.31$ kton$\cdot$year. In 2015, a new SK analysis
including neutron identification led to stronger constraints down to $13.3$
MeV neutrino energies Zhang  _et al._ (2015). Since the SK triggers did not
allow to record the neutron capture signal until 2008, this analysis could
only be performed on a small part of the SK data set, with a total livetime of
$22.5\times 960$ kton$\cdot$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 Bays  _et al._ (2012) above 17.3 MeV.

In this study, we draw on the previous SK analyses to present two DSNB
searches for antineutrino energies ranging from $9.3$ to $81.3$ MeV, with
significantly improved background modeling and reduction techniques. In the
$9.3$ to $31.3$ MeV range, we derive differential upper limits on the
$\bar{\nu}_{e}$ flux independently from the DSNB model, following the strategy
outlined in Ref. Zhang  _et al._ (2015) using a $22.5\times
2970$-kton$\cdot$day data set. In the $17.3$ to $81.3$ MeV range we constrain
a wide variety of DSNB models using spectral fits analogous to the ones
described in Ref. Bays  _et al._ (2012). We then combine the results of this
analysis with the ones obtained in Ref. Bays  _et al._ (2012) for the former
SK phases, thus analyzing $22.5\times 5823$-kton$\cdot$days of data. This
unprecedented exposure will allow to probe the DSNB with an 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 Section II.
Then, details about modeling of the DSNB signal and the different backgrounds
are given in Sections III and IV, respectively. We then describe the data
reduction process in Section V. The procedures associated with the DSNB model-
independent and spectral analyses are described in Sections VI and VII,
respectively. Finally we discuss the current constraints on the DSNB flux and
future opportunities at SK with gadolinium and Hyper-Kamiokande in Section
VIII before concluding in Section IX.

##  II Super-Kamiokande

Super-Kamiokande is a 50-kton water Cherenkov detector located in the Kamioka
mine, Japan. It is structured by a cylindrical stainless steel tank with a
diameter of $39.3$ m and height of $41.4$ m and consists of two parts: an
outer detector (OD) that serves as a muon veto and an inner detector (ID)
where neutrino detection takes place. In order to reduce backgrounds due to
radioactivity near the detector wall, most analyses consider only events
reconstructed at least 2 m away from the ID wall, thus defining a 22.5-kton
fiducial volume (FV). It is located 1000 m underground, that allows reduction
of the cosmic ray muon flux by a factor of $10^{5}$. In order to ensure a
high-quality data taking, conditions inside SK are tightly controlled; water
is constantly recirculated and purified, and the ID wall is covered with
11,129 20-inch photomultipliers (PMTs) with a $3$-ns time resolution,
corresponding to a $40\%$ photocathode coverage. These features allow SK to
detect particles with energies ranging from a few MeV to a few TeVs. The OD
includes 1,885 8-inch PMTs, facing outwards, to detect the Cherenkov light
from muons. Further detailed descriptions of the SK detector and its
calibration can be found in Refs. Fukuda  _et al._ (2003); Kume  _et al._
(1983); Abe  _et al._ (2014a); Nakahata  _et al._ (1999); Blaufuss  _et al._
(2001).

SK is currently undergoing its sixth data taking phase since it started
functioning in 1996. The first phase lasted $1497$ days and ended for a
scheduled maintenance. Due to an accident following the maintenance, which
resulted in a loss 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, that lasted $562$ days. SK’s previous
spectral analysis of the DSNB Bays  _et al._ (2012) used data from all these
three phases. Since 2008, the front-end electronics has been replaced Nishino
_et al._ (2009) and a new trigger system that allows for neutron tagging has
been set up Watanabe  _et al._ (2009). SK operated with this new electronics
during phase 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 calibration and monitoring purposes (phase V) since January 2019, before
being loaded with gadolinium (phase VI) in July 2020. This study will
primarily focus on phase IV, and will present a combined analysis of the SK-I
to IV data. Additionally, we will briefly discuss the opportunities offered by
the gadolinium lodaded SK and by the future Hyper-Kamiokande experiments.

Data processing in SK makes use of multiple triggers, corresponding to
different thresholds on the number of PMT hits found in a $200$-ns window. In
SK-I to III, all PMT hits in a 1.3-$\mu$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-high-energy (“SHE”)
trigger, which requires 58 PMT hits (70 hits before September 2011) in $200$
ns. The data in a $[-5,+35]$-$\mu$s window around the main activity peak is
collected with this trigger. This window is too short to contain the delayed
neutron capture signal after IBDs. In order to identify this signal, a
500-$\mu$s ($350~{}\mu$s before November 2010) after trigger (“AFT”) window
follows each SHE trigger that is not associated with a pre-determined OD
trigger. The trigger conditions for the different periods in SK-IV are
summarized in Table 1. In the rest of this paper, we will describe the
analysis of hit patterns in these large SHE+AFT windows to identify the
coincident production of a positron and a neutron.

The prompt positron event is reconstructed using the dedicated solar neutrino
Hosaka  _et al._ (2006); Cravens  _et al._ (2008); Abe  _et al._ (2011, 2016)
and muon decay electron vertex and direction fitter, which is then used to
reconstruct the event energy Smy (2008). For this SK phase, we follow the
convention introduced in Ref. Abe  _et al._ (2016) and subtract the 0.511 MeV
electron mass from this energy to obtain the electron equivalent kinetic
energy $E_{\rm rec}$. This quantity can also be interpreted as the total
reconstructed positron energy, and will be used to present the results of this
study.

Table 1: Summary on the SHE and AFT trigger conditions and livetime for the different three periods in SK-IV. SHE-threshold | AFT-window | Livetime
---|---|---
70 hits | 350 $\mu$s | 25.0 days
70 hits | 500 $\mu$s | 869.8 days
58 hits | 500 $\mu$s | 2075.3 days

##  III Signal modeling

###  III.1 DSNB spectra and kinematics

In this analysis we consider both the DSNB models whose fluxes are shown in
Fig. 1 and models where we vary specific physical parameters. For these
parameterized models we compute the DSNB flux using the following formula:

$\displaystyle\Phi(E_{\nu})=$ $\displaystyle\frac{c}{H_{0}}\int\sum_{s}R_{\rm
SN}(z,s)\,\sum_{\nu_{i},\bar{\nu}_{i}}F_{i}\left(E_{\nu}(1+z),s\right)$
$\displaystyle\times\frac{dz}{\sqrt{\Omega_{M}(1+z)^{3}+\Omega_{\Lambda}}},$
(1)

where $E_{\nu}$ is the neutrino energy, $c$ is the speed of light in vacuum,
$H_{0}$ is the Hubble constant, $z$ is the redshift, $R_{\rm SN}(z)$ is the
redshift-dependent supernova rate, and $F_{i}$ is the supernova neutrino
emission spectrum for a given flavor $i$. The index $s$ represents the
different possible classes of supernovae —e.g. supernovae collapsing into
neutron stars or black-holes— associated with specific neutrino emission
spectra. The last factor accounts for the Universe expansion, with
$\Omega_{M}$ and $\Omega_{\Lambda}$ being the matter and dark energy
contributions to the energy density of the Universe, respectively. We evaluate
the supernova rate using a phenomenological model described in Ref. Horiuchi
_et al._ (2009), which extracts the redshift dependence of the cosmic star
formation history by fitting observations and uses the Salpeter initial mass
function Salpeter (1955) to obtain the fraction of supernova progenitors. We
then consider blackbody neutrino emission spectra, where the effective
temperature $T_{\nu}$ takes into account neutrino oscillation effects. For a
given neutrino flavor $i$ the spectrum for a model with effective temperature
$T_{\nu}$ is thus given by:

$\displaystyle F_{i}(E_{\nu})$ $\displaystyle=f_{\nu,i}E_{\nu}^{\rm
tot}\,\frac{120}{7\pi^{4}}\,\frac{E_{\nu,i}^{2}}{T_{\nu,i}^{4}}\left(e^{E_{\nu,i}/T_{\nu,i}}+1\right)^{-1},$
(2)

where $E_{\nu}^{\rm tot}$ is the total energy of the neutrinos emitted by the
supernovae and $f_{\nu,i}$ represents the fraction of neutrinos or
antineutrinos with flavor $i$ and can be roughly approximated by $1/6$ for a
given species. The blackbody model therefore only depends on the neutrino
temperatures and luminosities. In what follows, we will refer 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 Keil  _et al._ (2003);
however, the current exposure at SK does not allow to probe this pinching
effect. In this analysis we neglect the effects of degeneracies between the
pinching parameter and e.g. the neutrino temperature on the final constraints
on neutrino emission spectra.

This analysis exclusively focuses on the detection of electron antineutrinos
via IBD processes. In the rest of this paper, we model these processes using
the Strumia-Vissani IBD calculation Strumia and Vissani (2003), that
approximates the IBD cross section more precisely than the Beacom-Vogel
calculation used for previous SK analyses Bays  _et al._ (2012).

###  III.2 Detector simulation

SK-IV has been the longest data taking phase of the experiment, 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 the stability of the detector. A key ingredient of our
DSNB study is therefore a simulation of antineutrino IBD processes that
accounts for the evolution of the detector’s properties throughout the entire
SK-IV period.

In this study we simulate IBD processes uniformly distributed in the entire
inner detector, in order to account for 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 the wide diversity of DSNB models, we generate uniformly
distributed positron energies in 1$-$90 MeV and renormalize events later to
model specific spectra. This procedure will also allow to use this simulation
to model backgrounds associated to other IBD processes or $\beta$ decays, as
will be mentioned later. We then simulate detection signals using a dedicated
simulation, based on GEANT3 Brun  _et al._ (1994), that 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 to accurately model the time-
dependent noise contamination of the positron signal. Neutron capture,
however, produces a $2.2$ MeV $\gamma$ ray whose light yield is below the SK
trigger threshold and particularly difficult to distinguish from PMT dark
noise and low energy radioactive decays. In order to develop robust neutron
identification techniques we therefore collect noise samples from data using a
random wide trigger, and inject them into simulation results, after the
positron activity peak.

##  IV Background Sources

###  IV.1 Atmospheric neutrinos

Atmospheric neutrinos are an important background in the present DSNB
searches. Below about 15 MeV, neutral-current quasielastic (NCQE)
interactions, which induce nuclear $\gamma$ ray emission Ankowski  _et al._
(2012); Abe  _et al._ (2014b); Wan  _et al._ (2019); Abe  _et al._ (2019),
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
associated reconstructed energies will follow a Michel spectrum for
reconstructed energies of about 15 to 50 MeV. Atmospheric neutrino events are
simulated using the HKKM 2011 flux Honda  _et al._ (2007, 2011); bib (Flux) as
an input to the neutrino event generator NEUT 5.3.6 Hayato (2009) to model
interactions. Details about the model parameters used in this analysis can be
found in Ref. Abe  _et al._ (2020). The nuclear deexcitation $\gamma$’s are
simulated based on the spectroscopic factors for the $p_{\rm 1/2}$, $p_{\rm
3/2}$, and $s_{\rm 1/2}$ states calculated in Ref. Ankowski  _et al._ (2012).
More detailed descriptions of the excited states can be found in Refs. Abe
_et al._ (2014b, 2019) for the T2K NCQE measurement. One difference between
our procedure and Ref. Abe  _et al._ (2019) is treatment of the others state
—a state affected by short-range correlations or with a very high excitation
energy—. This state is treated as the highest excitation state ($s_{\rm 1/2}$)
in Ref. Abe  _et al._ (2019) while it is included in the ground state ($p_{\rm
1/2}$) in the present analysis as it is from the latest atmospheric neutrino
analysis in SK Jiang  _et al._ (2019). The systematic uncertainty regarding
this treatment is considered in the scaling factor used in Section VI.1.

The detector simulation for atmospheric neutrino events is performed using the
same GEANT3-based simulation tool as for the signal, but without corrections
for the PMT gain shift on the primary event. The associated systematic error
is estimated in Section VI.1 and found to be subdominant. Finally, random
trigger data are overlaid on top of neutron capture signals in order to
simulate background for neutron tagging, following the procedure described in
Section III.2. Note that for this last step the time evolution of the detector
properties is taken into account.

###  IV.2 Cosmic ray muon spallation

In spite of its 2700-m water-equivalent overburden, SK is still exposed to
cosmic ray muons with a rate of $\sim$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 then lead to the
production of radioactive isotopes, whose decays can be misidentified as IBD
events in the DSNB search window. Below about 20 MeV, the associated
background is $10^{6}$ times higher than the DSNB flux predictions, making
spallation reduction an essential aspect of this analysis. Since most
spallation isotope decays only produce $\beta$ and $\gamma$ rays, spallation
backgrounds can be significantly reduced using neutron tagging. Nonetheless,
accidental pairing between prompt $\beta$ or $\gamma$ 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
$\beta$+$n$ decay that mimics the IBD signal, and thus cannot be removed using
neutron tagging. Efficient spallation reduction therefore requires both
dedicated spallation reduction techniques tailored to the properties of the
different isotopes, and neutron tagging cuts.

Relevant spallation isotopes and their visibility in water have been studied
with the FLUKA simulation Battistoni  _et al._ (2007) in Refs. Li and Beacom
(2014, 2015a, 2015b). 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
${\rm{}^{14}B}$ and ${\rm{}^{11}Li}$. The spallation cuts used in this
analysis will therefore be applied up to $23.5$ MeV reconstructed kinetic
electron energy in order to account for energy resolution effects. Note also
that the isotope half-lives range from $\mathcal{O}(0.01)$ sec to a few tens
of seconds. Given the high rate of cosmic ray muons, pairing a given isotope
decay with its parent process is a particularly difficult endeavour.

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

Isotopes that undergo $\beta$$+$$n$ decays cannot be rejected using neutron
tagging and therefore need to be modeled with particular care. This is notably
the case for 8He, 11Li, 16C, and 9Li. Here 11Li is particularly short-lived
and therefore easy to eliminate using mild spallation cuts. Its production
yield is also expected to be particularly low (around
$10^{-9}~{}\mu^{-1}\cdot$g${}^{-1}\cdot$cm2 Li and Beacom (2014)). In
addition, 8He and 16C will remain subdominant due to their low production
yields (around 0.23 and $0.02\times 10^{-7}~{}\mu^{-1}\cdot$g${}^{-1}\cdot$cm2
Li and Beacom (2014)) and endpoint $\beta$ decay energies. On the other hand,
9Li has a non-negligible yield (around $1.9\times
10^{-7}~{}\mu^{-1}\cdot$g${}^{-1}\cdot$cm2 Li and Beacom (2014)), a medium
half-life of about $0.18$ sec, and decays into a $\beta$$+$$n$ pair with a
branching ratio of 50.8%. The associated $\beta$ spectrum is shown in Fig. 4,
and extends up to $13.5$ MeV reconstructed energy when accounting for
resolution effects. The $\beta$+$n$ decay of 9Li is hence similar to the IBD
of a DSNB neutrino, except that the $\beta$ from ${\rm{}^{9}Li}$ is an
electron and the neutron energy is higher than that from IBDs. Since SK does
not distinguish electrons from positrons and does not allow to measure the
neutron energy, we model 9Li decay using the IBD Monte-Carlo simulation,
renormalizing events to reproduce the 9Li $\beta$+$n$ spectrum. Finally, we
estimate the total expected number of 9Li $\beta$+$n$ decays in the DSNB
search window by using an earlier SK measurement Zhang  _et al._ (2016) that
estimated the 9Li production rate to be $0.86\pm 0.12({\rm stat.})\pm
0.15({\rm syst.})~{}{\rm kton^{-1}}$$\cdot$${\rm day^{-1}}$. Note that the
difference in neutron energies with respect to IBD processes leads to a
systematic uncertainty. However, this uncertainty is taken into account when
calibrating neutron tagging efficiencies, as neutrons from the calibration
source have similar energies to the 9Li decay products. More details about
this calibration procedure are given in Section V.4.

Figure 4: True electron kinetic energy (blue) and reconstructed kinetic energy
(orange) spectra from the 9Li $\beta$+$n$ decay.

###  IV.3 Reactor neutrinos

Antineutrinos emitted from nearby nuclear reactors can undergo IBDs in SK and
therefore be mistaken as the DSNB signal at low energies. We estimate the
reactor antineutrino flux using the SKReact code bib (eact) developed for the
SK reactor neutrino analysis. This code uses the reactor neutrino model from
Ref. Baldoncini  _et al._ (2015) based on IAEA records int (2005) to evaluate
the electron antineutrino fluxes for any reactor in the world, and simulates
oscillation effects. The resulting flux prediction also accounts for the time
dependence of the reactor antineutrino flux due to changes in reactor
activity. The expected total reactor neutrino flux at SK can thus be predicted
up to $E_{\nu}^{\text{max}}=9.2$ MeV. The expected reactor neutrino spectrum
is obtained from this predicted flux using the IBD simulation, as shown in
Fig. 5. Due to energy resolution effects, this spectrum extends up to $\sim$10
MeV reconstructed energy and dominates over DSNB predictions there. In the
7.5$-$9.5 MeV reconstructed kinetic electron energy range, for example, the
reactor antineutrino background is about 6 times higher than the DSNB signal
predicted by the Horiuchi$+$09 6 MeV model Horiuchi  _et al._ (2009). Reactor
antineutrino backgrounds thus effectively set a lower energy threshold for the
DSNB search.

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

###  IV.4 Solar neutrinos

Electron neutrinos from the Sun cause elastic interactions with electrons in
SK. In the DSNB search window, the 8B and hep solar neutrino fluxes represent
an important background up to $\sim$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 absence of
neutron tagging, a solar neutrino interaction can still be identified by
comparing the reconstructed directions of the scattered electron to the
direction of the Sun at the detection time. The angle between these two
directions, called the opening angle, is expected to 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 for the latest SK-IV solar neutrino analysis Abe
_et al._ (2016) 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 from the DSNB signature by spectral
shape studies. So far, only upper limits on the solar antineutrino flux have
been derived at SK, notably in Ref. Abe  _et al._ (2020).

##  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\times 2970.1$ kton$\cdot$days. The SHE
trigger efficiency is close to $100$% over the entire analysis window. The
following sections describe the four reduction steps used in this analysis:
noise reduction in Section V.1, spallation reduction in Section V.2, DSNB
positron candidate selection in Section V.3, and neutron tagging in Section
V.4.

###  V.1 Search energy range and noise reduction

The lower energy threshold for this analysis is the SHE trigger threshold,
which corresponds to $E_{\rm rec}>$ 7.5 or 9.5 MeV depending on the
observation time (see Table 1). The upper bound of the DSNB analysis window is
set to $29.5$ MeV for the model-independent analysis and $79.5$ MeV for the
spectral analysis.

We first select SHE-triggered events with $E_{\rm rec}$ below $79.5$ MeV and
apply a set of cuts aimed at removing PMT noise, 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 fiducial volume (FV) of $22.5$
kton. Additionally, we remove events associated with an OD trigger, or within
$50~{}\mu$s of another event with reconstructed electron kinetic energy larger
than $6$ MeV, in order to remove electrons produced by the decay of low energy
cosmic ray muons stopping in the detector. In addition, we apply fit quality
cuts to eliminate noise-like events. Inside the FV these noise reduction cuts
have a signal efficiency larger than $99\%$ as confirmed by the IBD Monte-
Carlo simulation. Finally, we require all SHE events to be followed by an AFT
trigger to allow for neutron tagging. Due to occasional deadtime induced by
trigger software failure, this step is associated with a $\sim$94% signal
efficiency. Events passing these criteria will be subsequently referred to as
“DSNB candidates”.

###  V.2 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 dedicated cuts that will allow for a substantial reduction of
these backgrounds, even in the absence of neutron tagging.

####  V.2.1 Spallation preselection

We apply two sets of preselection cuts to reduce contributions from the most
energetic muons and some long-lived spallation isotopes such as 16N, with
minimum harm to the signal efficiency.

First, since energetic muons could induce the production of multiple
radioactive isotopes, we remove all DSNB candidates observed within 60 sec and
4.9 m from at least one other low energy event with $E_{\rm rec}$ in the
5.5-24.5 MeV range Locke  _et al._ (2021). 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 a sample of low energy events from radioactive decays at SK ($4<E_{\rm
rec}<5$ MeV), whose reconstructed vertices have been replaced by random
vertices inside the ID. The cut parameters introduced above allow to remove
about $45\%$ of the background with a $98\%$ signal efficiency, as can be seen
in Fig. 6. Since this cut removes low energy events clustered in space and
time, it also removes low energy muons that are misclassified as electrons
with tens of MeV energies. Consequently, in this study, we will apply this cut
over the whole energy range for each of our analyses, including the
spallation-free $E_{\rm rec}>23.5$ MeV region.

Figure 6: Minimal distance between a given DSNB candidate (blue) or SK low
energy event with a random vertex (orange) and a well-reconstructed event with
$E_{\rm rec}>5.5$ MeV observed within $60$ sec of it. The DSNB candidates are
defined as all the SK-IV events with $E_{\rm rec}>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.

In addition to removing multiple spallation events, we locate muon-induced
showers by identifying neutron captures observed less than $500~{}\mu$s after
muons. These neutron clusters, called “neutron clouds”, are often produced in
muon-induced hadronic showers, that are the birthplaces of spallation
isotopes. Neutron capture events following muons are collected separately from
DSNB candidates, using a specific trigger called the Wideband Intelligent
Trigger (WIT) Carminati (2011, 2015); Elnimr (2017). This trigger, in addition
to setting a threshold on the number of PMT hits, applies cuts on the event
quality and distance of the reconstructed vertex to the ID wall. Events
triggered by WIT less than $500~{}\mu$s after a given muon are considered as
neutron captures if their reconstructed vertex is within $5$ m of the
reconstructed muon track. If multiple neutron captures are observed after a
muon event, we define a neutron cloud using the criteria from Ref. Locke  _et
al._ (2021). We then remove DSNB candidates that are found close in time and
space to these neutron clouds, following the procedure introduced in Ref.
Locke  _et al._ (2021). More details about the cut criteria are given in
Appendix A. To estimate the signal and background efficiencies of the neutron
cloud cut, we pair muons with DSNB candidates in data, defining two samples: a
“pre-sample”, including muons found up to $60$ sec _before_ a DSNB candidate,
and a “post-sample”, with muons found up to $30$ sec _after_ a candidate.
While the pre-sample contains a mixture of correlated and uncorrelated pairs,
the post-sample contains only uncorrelated pairs. This sample can hence be
used to both readily estimate signal efficiencies and characterize the
properties of the pairs formed by each isotope decay and their parent muon in
the pre-sample. An example of how to extract the distribution of the distances
between spallation events and the neutron clouds of their parent muons is
shown in Fig. 7. We optimize the neutron cloud cut by estimating the shape of
the neutron clouds from data and maximizing 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 trigger has only been available only during the last 388
live days of SK-IV. While the signal efficiency is larger than $99$% this cut
only removes about $10$% of the spallation background in this analysis. Over
the 388 days of the WIT trigger livetime, however, it removes about 40% of
spallation isotope decays.

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

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

####  V.2.2 Searching for parent muons

Since preselection cuts remove only about half of the spallation background, a
more in-depth study of the correlations between muons and DSNB candidates is
necessary. In particular, in order to identify the decays of spallation
isotopes, it is crucial to associate each isotope with its parent muon. Due to
the long half-lives of isotopes such as 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 to accommodate most decay products,
without processing a prohibitively high number of muons. As with neutron cloud
cut optimization, we define pre- and post-samples in order to extract the
observable distributions 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. Additionally, since the SK trigger windows can contain multiple 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 direction of the track, are then extracted
using a dedicated fitter Conner (1997); Desai (2004); Bays  _et al._ (2012).
This fitter classifies muons into five categories: misfits (1.0%), single
through-going muons (82.2%), stopping muons (4.9%), multiple-track muons
(7.6%), and corner clippers (4.3%); the 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 in order to reject non-muonic high energy
events. This cut does not affect removing spallation background because the
track length of these low charge muons through the ID is typically less than
50 cm.

After having paired DSNB candidates with the muons selected above, we extract
the following observables, related to the intrinsic properties of the muons
and their correlations with DSNB candidates.

  * •

$dt$: time difference between a muon and the DSNB candidate. For spallation
events, this observable reflects the half-lives of the produced isotopes.

  * •

$\ell_{t}$: transverse distance of the DSNB candidate to the muon track (see
Fig. 8). This observable is related to the path lengths of the secondary
particles from muon showers. For well-fitted single through-going muons paired
with their associated spallation isotope decays, $\ell_{t}$ is typically no
larger than a few meters.

  * •

$\ell_{l}$: longitudinal distance to the point of the muon track associated
with the maximal energy deposition (see Fig. 8). The energy deposition along
the muon track is determined by projecting back the light seen by each PMT to
a point on the muon track, following the procedure described in Ref. Bays  _et
al._ (2012). So $\ell_{l}$ provides an estimate of the distance between the
low energy event and the origin of the particle shower and, if the shower
position is correctly identified, should also not be larger than a few meters
for spallation pairs.

  * •

$Q_{\mu}$: total charge deposited by muons in the detector. Higher values of
this observable are expected when a shower is produced.

  * •

$Q_{\rm res}$: residual charge deposited by a muon compared to the value
expected from the minimum ionization. This observable can be expressed as:

$Q_{\rm res}=Q_{\mu}-q_{\rm MI}\cdot L,$

where $q_{\rm MI}$ and $L$ are the number of $p.e.$ per centimeter expected
from the minimum ionization and the track length, respectively. Here $q_{\rm
MI}$ depends on detector properties, such as water transparency, that vary
with time and is therefore recomputed for each run. For muons with multiple
tracks, $L$ is the length of the first track. This observable allows to
determine how likely a given muon is to induce a shower.

Figure 8: Schematic of the spallation observables $\ell_{t}$ and $\ell_{l}$.
The cylindrical shape represents the SK ID.

An example set of these spallation variables for single through-going muons is
shown in Appendix A. The distributions associated with corner clipping muons
show no evidence for spallation. Moreover, as discussed in Ref. Locke  _et
al._ (2021), spallation background events associated with corner 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 from these muons to spallation backgrounds.

Showers induced by muon spallation can be extremely energetic and involve up
to thousands of particles, notably neutrons and pions. Such neutrons cause
nuclear reactions and possibly produce deexcitation $\gamma$ rays, and are
later captured at a timescale up to $\sim$500 $\mu$s. The MeV-scale $\gamma$
rays or $\beta$ particles —produced by the initial neutron, secondary nuclear
reactions, or the decay product of a short-lived isotope, followed by a 2.2
MeV $\gamma$ ray from neutron capture can mimic an IBD signature. To prevent
these from leaking into the sample, we require DSNB candidates to be more than
$1$ ms 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 $\ell_{t}$ for
different energy ranges, exploiting the dependence of the isotopes’ half-lives
in their endpoint 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 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 (${\rm PDF}_{\rm
spall}^{i}$) and the random PDF (${\rm PDF}_{\rm random}^{i}$), for each
variable $i=dt,\ \ell_{t},\ \ell_{l},\ Q_{\mu},\ Q_{\rm res}$. We isolate
contributions from spallation events by subtracting the post-sample
distributions from the corresponding pre-sample distributions and obtain ${\rm
PDF}_{\rm spall}^{i}$ after area normalization of these subtracted
distributions. For ${\rm PDF}_{\rm random}^{i}$ we normalize post-sample
distributions by their areas. We repeat this procedure for each category of
muons. For single through-going and multiple track muons, that generate most
of the spallation background, we tune the PDFs in $dt$ and $\ell_{t}$ bins.
These bins, adjusted by considering the isotopes’ typical half-lives and the
$\ell_{t}$ distributions, account for the correlations between these
observables. For each set of PDFs, we then define log-likelihood ratios as:

$\displaystyle\mathcal{L}_{\text{spall}}=\log\left(\prod_{i}\frac{{\rm
PDF}_{\rm spall}^{i}}{{\rm PDF}_{\rm random}^{i}}\right).$ (3)

Separate likelihoods are defined for each muon category, except for corner
clippers, whose contributions to spallation backgrounds are negligible as
mentioned above. Note that for misfit muons only $dt$ is used to calculate
log-likelihood ratios, as the other observables are not reliable. An example
distribution of the log-likelihood ratio is given as Fig. 33 in Appendix A. We
finally determine cut conditions for each likelihood ratio, accounting for
complex correlations between spallation observables, 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-dominated data samples. These samples as well as our methodology to
estimate the spallation cut performance are described in detail in Appendix A.
The current cuts achieve significant improvement over the previous cuts used
in Refs. Bays  _et al._ (2012); Zhang  _et al._ (2015), especially at low
energies; for the same spallation remaining rate, the signal efficiency is
increased by up to 60% for $E_{\rm rec}<11.5$ MeV, 20% for $11.5<E_{\rm
rec}<13.5$ MeV, and comparable for higher energies.

###  V.3 DSNB positron candidate selection

In the region above $\sim$20 MeV, spallation backgrounds can be reduced to
negligible levels using a series of spallation cuts while keeping most of the
signal. However, significant backgrounds from atmospheric neutrino
interactions and radioactive decays remain. To identify them, we define the
following discriminating observables, aimed at characterizing the prompt
event.

####  V.3.1 Incoming event cut

Radioactivity near the detector wall as well as muon spallation in the rock
surrounding the detector, can lead to electrons or $\gamma$ rays entering FV.
Instead of tightening the FV cut, we consider the effective distance of each
event to the ID wall $d_{\text{eff}}$ Bays  _et al._ (2012); Hosaka  _et al._
(2006). This observable is computed by following the reconstructed direction
of each event backwards from its reconstructed vertex to the ID wall, and can
be 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
$d_{\text{eff}}$ distributions for the DSNB signal and for the events selected
using the cuts outlined in Section V.1 are shown in Fig. 10, for the
reconstructed energy ranges corresponding to the model-independent analysis,
$7.5$ to $29.5$ MeV, and the spectral analysis, $15.5$ to $79.5$ MeV. In this
study we impose lower thresholds on $d_{\text{eff}}$ ranging from $3$ to $5$ m
depending on the reconstructed energy:

$\displaystyle d_{\text{eff}}>\text{max}\left\\{300,\ 500-\frac{E_{\rm
rec}-15.5~{}\text{MeV}}{1~{}\text{MeV}}\times 50\right\\}~{}\text{cm}.$

Figure 9: Views of the ID, showing the effective distance of a reconstructed
event to the detector wall. The $d_{\text{eff}}$ is the distance that a
particle emitted near the ID wall —e. g. from by a radioactive decay— needs to
travel to give the observed signal.

Figure 10: The $d_{\text{eff}}$ distributions for the observed data and IBD
Monte-Carlo simulation, for events with $E_{\rm rec}\in[7.5,29.5]$ MeV (left)
and $E_{\rm rec}\in[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 $d_{\text{eff}}>1000$ cm. The lower panels show the relative
difference between the numbers of data and IBD simulation events. We consider
only events reconstructed in the FV and passing the noise reduction cuts. The
assumed DSNB spectrum is taken from the Horiuchi$+$09 6 MeV model Horiuchi
_et al._ (2009).

####  V.3.2 Pre- and post-activity cuts

Atmospheric neutrino interactions can produce both prompt signals from
photons, electrons, or energetic heavy particles, and delayed signals from
decay of muons and pions to electrons. These signals will be typically
separated by a few microseconds and can therefore share the same trigger
window. Depending on which particle deposits most light, unusually high
activity can be noticed either before or after the main peak. For pre-activity
we compute the maximal number of hits, $N_{\rm pre}^{\rm max}$, in a $15$-ns
time-of-flight subtracted window between $5~{}\mu$s and $12$ ns before the
main peak. We apply a similar procedure to post-activity, using an algorithm
developed for other analyses in SK and T2K that computes the number of
electrons from muon and pion decay $N_{{\rm decay}\mathchar 45\relax e}$ Abe
_et al._ (2017). The $N_{{\rm decay}\mathchar 45\relax e}$ distributions for
low and high energy events are shown in Fig. 11. In the analyses discussed
here we require $N_{\rm pre}^{\rm max}<12$ and $N_{{\rm decay}\mathchar
45\relax e}=0$.

Figure 11: The $N_{{\rm decay}\mathchar 45\relax e}$ distributions for the
atmospheric neutrino Monte-Carlo simulation, for events with $E_{\rm
rec}\in[7.5,29.5]$ MeV (left) and $E_{\rm rec}\in[15.5,79.5]$ MeV (right). We
consider only events reconstructed in the FV and passing the noise reduction
cuts.

####  V.3.3 Cherenkov angle

The opening angle of the Cherenkov light cone ($\theta_{\rm C}$) emitted by
highly relativistic particles, electrons and positrons in the current context,
in pure water is around $42^{\circ}$. Conversely, at the $\mathcal{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,
that will be reconstructed as a single cone with a particularly large opening
angle Abe  _et al._ (2014b, 2019). Consequently, $\theta_{\rm C}$ is one of
the most powerful observables for reducing both NCQE and $\mu/\pi$-producing
atmospheric backgrounds. The Cherenkov angle distributions for the signal and
atmospheric neutrino events are shown in Fig. 12 for the energy ranges of the
model-independent analysis and the spectral analysis. In what follows we will
require the Cherenkov angle in the signal regions to be in $[38,50]^{\circ}$.
In the spectral analysis, we will also use this angle to define sidebands to
evaluate atmospheric backgrounds.

Figure 12: The $\theta_{\rm C}$ distributions for the atmospheric neutrino and
IBD Monte-Carlo simulations, for events with $E_{\rm rec}\in[7.5,29.5]$ MeV
(left) and $E_{\rm rec}\in[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 Horiuchi  _et al._
(2009). The signal and background distributions are normalized to the same
area.

####  V.3.4 Ring clearness

Electrons and positrons do not follow a straight trajectory in SK due to
scattering and bremsstrahlung, which leads 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 compute the opening angles of the
cones formed by the directions of all possible 3-hit combinations. We then
identify the peak of this opening angle distribution $\theta_{0}$ and estimate
the “clearness” of the ring by computing:

$\displaystyle L_{\texttt{clear}}=\frac{N_{\text{triplets}}(\theta_{0}\pm
3^{\circ})}{N_{\text{triplets}}(\theta_{0}\pm 10^{\circ})}.$ (4)

The $L_{\texttt{clear}}$ distributions from the atmospheric and IBD Monte-
Carlo simulations are shown in Fig. 13. Events with lower $L_{\texttt{clear}}$
have fuzzier rings. For the analyses presented here we require
$L_{\texttt{clear}}<0.36$.

Figure 13: The $L_{\texttt{clear}}$ distributions for the atmospheric neutrino
and IBD Monte-Carlo simulations, for events with $E_{\rm rec}\in[7.5,29.5]$
MeV (left) and $E_{\rm rec}\in[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 Horiuchi
_et al._ (2009). The signal and background distributions are normalized to the
same area.

####  V.3.5 Average charge deposit

Energetic muons scatter less than electrons and hence often deposit more
charge in individual PMTs than $\mathcal{O}(10)$ MeV electrons and positrons.
We exploit this feature by considering a $50$-ns time-of-flight subtracted
window around the main activity peak and calculate the average charge
deposited per PMT hit ($q_{50}/n_{50}$) in this window. The distributions are
given in Fig. 14. We require $q_{50}/n_{50}$ to be less than $2$ in the
analysis.

Figure 14: The $q_{50}/n_{50}$ distributions for the atmospheric neutrino and
IBD Monte-Carlo simulations, for events with $E_{\rm rec}\in[7.5,29.5]$ MeV
(left) and $E_{\rm rec}\in[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 Horiuchi  _et al._
(2009). The signal and background distributions are normalized to the same
area.

####  V.3.6 Cut efficiencies and systematic uncertainties

We estimate the signal efficiencies of most of the cuts detailed above using
the IBD Monte-Carlo simulation. Hence, the main source of systematic
uncertainties on the cut efficiencies stems from the modeling of the SK
detector and particle propagation in water. We estimate these uncertainties
using calibration data taken by injecting a monoenergetic electron beam
produced by a LINAC at different positions inside the detector Nakahata  _et
al._ (1999). For a given cut, we compare the measured and predicted
efficiencies using a dedicated Monte-Carlo simulation, and define the
systematic uncertainty as the maximum discrepancy over all calibration tests.
LINAC calibration thus allows to estimate uncertainties for ring clearness,
the average charge, and Cherenkov angle cuts. Conversely, for the incoming
event cut —that requires considering uniformly distributed events— we use a
strategy developed for the SK solar neutrino analyses Hosaka  _et al._ (2006);
Cravens  _et al._ (2008); Abe  _et al._ (2011, 2016): we measure the variation
of the $d_{\text{eff}}$ cut efficiency after shifting the event directions and
vertices in the IBD Monte-Carlo simulation. The vertex and direction shifts
used for these estimates are obtained from calibration studies using a nickel
source Abe  _et al._ (2014a). Finally, efficiencies for the pre- and post-
activity cuts can be directly estimated using a spallation-dominated sample of
low energy events, with negligible systematic uncertainty. Overall, positron
candidate selection cuts allow 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 are of a few percent.

###  V.4 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 $\mu$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. Zhang (2012, 2015).

####  V.4.1 DSNB candidate selection and neutron preselection

As explained in Section II, most SHE triggered events are followed by a 500
$\mu$s AFT trigger window in order to record the neutron capture signal.
Before attempting to identify neutron candidates in this combined trigger
window, we apply a cut on the maximum hits inside a $200$-ns window:
$N_{200}<50$ to eliminate SHE+AFT trigger 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 trigger windows of the remaining events to
select hit clusters that could be associated with neutron captures. Here, we
take advantage of the fact that a typical neutron capture vertex is within a
few tens of centimeters of the well-reconstructed IBD vertex. Since the SK
resolution is about 50 cm for $\mathcal{O}(10)$ MeV events, this proximity
allows to reliably estimate the required photon emission time for each PMT hit
to be the result of the neutron capture. A neutron candidate is therefore
defined as a group of hits that is 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}\geq 7$ Zhang (2012). 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 this study,
we modify this preselection step as follows. First, we apply a cut to reduce
the continuous PMT dark noise, likely due to scintillation in the PMT glass,
that causes a single PMT to flash more than once in rapid succession with a
timescale of $\sim$10 $\mu$s. To do so, for each PMT we remove hit pairs
separated by less than 6 (12) $\mu$s for $N_{10}>6$ ($N_{10}\leq 6$). After
applying this cut, we define neutron candidates as hit clusters verifying
$N_{10}\geq 6$. This preselection cut will allow loosening of the neutron
tagging cut in energy regions where spallation no longer dominates, in
particular for the spectral analysis. Finally, since the continuous dark noise
cut relies on investigating PMT hit pairs, its efficiency sharply increases
near the edge of the SHE+AFT trigger window. We restrict the neutron search
window to [$14,523$] $\mu$s ([$14,373$] $\mu$s for events with a 350-$\mu$s
AFT trigger window). After these different steps, neutron candidate
preselection is associated with an efficiency of $44.7\%$.

####  V.4.2 Boosted Decision Tree

After preselection, we are still left with a large background from e. g. dark
noise fluctuations and radon decays Nakano  _et al._ (2020), totaling about 7
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 the neutron candidate, and its position with
respect to the primary event vertex. These $22$ observables are then used in a
Boosted Decision Tree (BDT) implemented using the xgboost Python library Chen
and Guestrin (2016). The BDT is trained using a data set of $2.8\times 10^{8}$
neutron candidates. In this sample, about $2\times 10^{6}$ candidates are
neutron captures simulated using the procedure described in Section III.2 and
the rest is composed of low energy accidental coincidences from random trigger
events. Table 2 summarizes the BDT training parameters. The final BDT
performance is shown in Fig. 15. Even accounting for preselection cuts, this
neutron tagging algorithm allows reduction the background rate by a factor of
$7$ compared to the analysis presented in Ref. Zhang (2015) for the same
signal efficiency.

Table 2: Parameters used for the training of the neutron tagging BDT in the xgboost Python library framework. Parameter | Value
---|---
learning rate | 0.02522
subsample | 0.97
max depth | 10
tree method | approx
max iterations | 1500
best iteration | 1170
Figure 15: Neutron identification efficiency and mistag rate of BDT used for
the current analysis. The red band depicts the uncertainty associated with the
time-dependence of the noise in SK, as described in Section VI.1.4. The
fraction of remaining IBD neutrons and the number of accidental coincidences
per event after preselection, averaged over the entire SK-IV period, are
$44.7\%$ and 7, respectively. These values depend linearly on the PMT noise
Akutsu (2019), and hence vary by about 15% over the SK-IV period.

####  V.4.3 Systematic uncertainty on neutron tagging efficiency

To estimate the systematic uncertainty on the BDT efficiency associated with
mismodeling of neutron propagation in water and the $2.2$ MeV photon emission,
we compare the efficiencies predicted by the IBD Monte-Carlo simulation to
measurements performed using an americium-beryllium (241Am/Be) radioactive
source embedded in a BGO scintillator. A detailed description of this
procedure is given in Ref. Watanabe  _et al._ (2009). To calculate the neutron
identification efficiency for data, we follow the methodology described in
Ref. Abe  _et al._ (2021a) and fit the timing distribution of tagged neutrons
by a decaying exponential plus a constant term as shown in Fig. 16. True
neutron capture event counts decay exponentially in time 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.

Figure 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\pm 6.1$ $\mu$s, is consistent with prior measurements.

###  V.5 Solar neutrino cut

Since even mild neutron tagging cuts can reduce solar neutrino backgrounds to
negligible levels, dedicated solar neutrino cuts are not implemented for the
DSNB model-independent analysis. For the spectral analysis, however, events
with no tagged neutrons are also considered in order to 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 $\theta_{\text{sun}}$ between the reconstructed direction of
candidate events and the direction of the Sun. The cosine of this angle is
peaked at 1 for solar neutrino interaction products, and is smeared by
kinematic and reconstruction effects. To estimate the impact of the latter, an
specific observable, the Multiple Scattering Goodness (MSG) $g$ has been
designed for the solar neutrino analysis Hosaka  _et al._ (2006); Cravens  _et
al._ (2008); Abe  _et al._ (2011, 2016). Simulated opening angle distributions
for solar neutrino events in different MSG bins are shown in Fig. 17.

Figure 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.

In this analysis we use the method described in Ref. Bays  _et al._ (2012) and
apply $\cos\theta_{\rm 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 $\theta_{\rm 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 to higher energies using the solar neutrino MC simulation from Ref. Abe
_et al._ (2021b), a procedure already used in the DSNB search at SK-I,II,III
Bays  _et al._ (2012). Here, we can safely assume the
$\cos\theta_{\text{sun}}$ distribution for non-solar event to be uniform Abe
_et al._ (2016), and therefore approximate the number of solar events in
13.5-15.5 MeV by:

$\displaystyle
N^{\text{solar}}_{\mathrm{low}}=N_{\mathrm{low}}(\cos\theta_{\text{sun}}>0)-N_{\mathrm{low}}(\cos\theta_{\text{sun}}<0),$
(5)

using the data passing the noise, spallation, and positron selection cuts
described in Secs V.1, 3, and V.3. Using the solar neutrino MC simulation from
Ref. Abe  _et al._ (2021b), the predicted number of solar neutrino events
above 15.5 MeV will then be

$\displaystyle N^{\text{solar}}\approx 0.12N^{\text{solar}}_{\mathrm{low}}.$
(6)

The impact of solar and neutron tagging cuts can then be assessed by rescaling
$N^{\text{solar}}$ by the corresponding efficiencies.

Finally, we evaluate the systematic uncertainties on the signal efficiency
with the solar event cut using the LINAC calibration results. Specifically, we
compute the scaling factor that allows to minimize the difference between the
predicted and observed MSG distributions for LINAC events. We then evaluate
the impact of this rescaling on the final efficiency to be of about 1%.

###  V.6 Cut optimization

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

####  V.6.1 DSNB positron candidate selection cuts

We determined the positron candidate selection cuts by comparing the IBD
Monte-Carlo simulation with the atmospheric neutrino simulation for the ring
clearness, average charge deposit, and Cherenkov angle, and the observed data
for the effective distance $d_{\text{eff}}$. Since the impact of the positron
candidate selection cuts depends weakly on the DSNB model considered, we
assumed the Horiuchi$+$09 spectrum Horiuchi  _et al._ (2009). As discussed in
Section V.3, we estimate systematic uncertainties using results from the solar
neutrino analysis Hosaka  _et al._ (2006); Cravens  _et al._ (2008); Abe  _et
al._ (2011, 2016) and the LINAC calibration Nakahata  _et al._ (1999). The
final cuts are the following:

$\displaystyle d_{\text{eff}}>\text{max}\left\\{300,\ 500-\frac{E_{\rm
rec}-15.5~{}\text{MeV}}{1~{}\text{MeV}}\times 50\right\\}~{}\text{cm},$
$\displaystyle\qquad\qquad\qquad\qquad N_{\rm pre}^{\rm max}<12,$
$\displaystyle\qquad\qquad\qquad\qquad N_{{\rm decay}\mathchar 45\relax e}=0,$
$\displaystyle\qquad\qquad\qquad\qquad L_{\texttt{clear}}<0.36,$
$\displaystyle\qquad\qquad\qquad\qquad q_{50}/n_{50}<2,$
$\displaystyle\qquad\qquad\qquad\qquad 38^{\circ}<\theta_{\rm C}<50^{\circ}.$

The efficiencies of these cuts are shown in Fig. 18 as a function of $E_{\rm
rec}$ and vary between 70 and 85%. Overall, positron candidate selection cuts
allow to remove $80\%$ of the atmospheric neutrino backgrounds, and bring
contamination from natural radioactivity to negligible levels.

####  V.6.2 Spallation and neutron tagging cuts

Unlike positron candidate selection cuts, spallation and neutron tagging cuts
need to be optimized in multiple energy regions to account for the important
variations of the spallation and atmospheric background rates in the analysis
windows. While only events with one tagged neutron are used 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 and neutron
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 $\neq$ 1 tagged neutron.

##### Events with one tagged neutron:

We devise a cut optimization scheme common to both the model-independent and
spectral analyses by simultaneously adjusting spallation and neutron tagging
cuts in $E_{\rm rec}$ 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 background-only hypothesis. We define the
sensitivity as the predicted number of signal events after cuts over the 90%
C.L. upper limit on the number of signal events, computed using the Rolke
method Rolke  _et al._ (2005). To simplify the interpretation of the results
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
the two analyses overlap. In particular, we require spallation and solar
neutrino backgrounds to be negligible above $15.5$ MeV. The remaining rates of
these backgrounds can be reliably evaluated by rescaling the number of events
observed after noise reduction, spallation, and positron candidate selection
cuts by the neutron mistag rate.

The signal efficiencies for our choice of cuts are shown in 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 Sections V.2,
V.3, and V.4, in Table 3.

Figure 18: DSNB signal efficiencies for the signal regions where one tagged neutron (top) and $\neq$ 1 tagged neutrons (bottom) are required, as started after the SHE$+$AFT trigger requirement and the noise reduction cuts. While the one tagged neutron region is common to both the model-independent and spectral analyses above 15.5 MeV, the $\neq$ 1 tagged neutron region is used only in the spectral analysis. Here, we applied noise reduction, spallation, third reduction, neutron tagging, and solar cuts sequentially and show cumulative efficiencies at each stage. The result in the final bin above 29.5 MeV corresponds to the efficiency for the range through 29.5 to 79.5 MeV, as calculated assuming the Horiuchi$+$09 6 MeV spectrum Horiuchi  _et al._ (2009). Note also that efficiencies for the solar neutrino cut in the bottom are shown in 1-MeV bins. Table 3: Systematic uncertainties for the reduction cuts, computed as described in Sections V.2, V.3, and V.4. The systematic uncertainties for the AFT trigger, noise reduction, pre- and post-activity cuts are negligible. The solar neutrino cut is used only for the spectral analysis, when events have $\neq$1 tagged neutrons. Cut | Relative systematic error
---|---
Spallation cut | $0.1\%$
$d_{\rm eff}$ cut | $0.1\%$
$L_{\texttt{clear}}$ cut | $0.2\%$
$q_{50}/n_{50}$ cut | $1.2\%$
$\theta_{\rm C}$ cut | $0.7\%$
Neutron tagging | $12.5\%$
Solar neutrino cut |  $1.0$%

##### Events with zero or more than one tagged neutrons:

As previously discussed, the events rejected by the neutron tagging cuts
obtained above can still be considered in the spectral analysis. In the
23.5$-$79.5 MeV region, that covers most of the spectral analysis window,
spallation and solar neutrino backgrounds are negligible and only noise
reduction and positron candidate selection cuts need to be applied to these
events. In the 19.5$-$23.5 MeV region, solar neutrino backgrounds are still
negligible and spallation backgrounds are mostly linked to short-lived
isotopes, that are easy to reject using the cuts from Section V.2. 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 and spallation backgrounds largely dominate and need to be reduced
using the targeted reduction strategies described in Sections V.2 and V.5. In
this analysis, we remove the solar neutrino backgrounds using the cuts in
Section V.5 with the cut points from Ref. Bays  _et al._ (2012). The cut
points and the corresponding signal efficiencies are shown in Table 4 for the
different multiple scattering goodness and reconstructed kinetic electron
energy regions. The remaining number of solar neutrino events after all cuts
is less than 2 for the entire energy region.

Table 4: Cuts on $\cos\theta_{\rm sun}$ (upper limits), and signal efficiencies for SK-IV. These cuts are only applied to events with $\neq 1$ tagged neutron, in the spectral analysis. They depend on the value of the multiple scattering goodness, $g$, and on the reconstructed kinetic electron energy. They are the same as the ones applied for SK-I and III in Ref. Bays  _et al._ (2012). The signal efficiency for SK-IV is within 1% of the one obtained in Ref. Bays  _et al._ (2012) for SK-I and III, which is similar to the associated systematic uncertainty. $E_{\rm rec}$ [MeV] | 15.5$-$16.5 | 16.5$-$17.5 | 17.5$-$18.5 | 18.5$-$19.5
---|---|---|---|---
$g<0.4$ | 0.05 | 0.35 | 0.45 | 0.93
$0.4<g<0.5$ | 0.39 | 0.61 | 0.77 | 0.93
$0.5<g<0.6$ | 0.59 | 0.73 | 0.81 | 0.93
$g>0.6$ | 0.73 | 0.79 | 0.91 | 0.93
Efficiency | 73.1% | 82.2% | 88.3% | 96.6%

In contrast, the impact of spallation cuts is difficult to evaluate 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 from
short-lived isotopes —with $\mathcal{O}(0.01)$ sec half-lives— to a negligible
level, using the cuts shown. in Appendix A. As can be seen in Fig. 3, the
remaining event sample will be dominated by decays from 8B and 8Li, that have
a half-life close to 1 sec. We then optimize the spallation likelihood cuts
using the method described in Appendix A.4 to compute the remaining amount of
8B and 8Li decays, and maximize $S/\sqrt{B}$ where $B$ and $S$ are the numbers
of spallation and non-spallation events respectively. The efficiencies of the
optimal cuts are of 65% and 63% in 15.5$-$17.5 and 17.5$-$19.5 MeV, for a
8B+8Li remaining fraction of about 2.5%. The corresponding remaining number of
spallation events, as well as the performance of this cut will be discussed
and compared to results from the previous SK DSNB analyses in Section VII.5.

###  V.7 Selected events in the final data samples

Here we discuss the behavior of the observed data after the cuts. In
particular, the impact of spallation, $d_{\mathrm{eff}}$, and the other
positron candidate selection cuts on data is shown in Fig. 19. This figure
allows to estimate the amount of spallation and radioactive backgrounds that
need to be removed, since these backgrounds cannot be currently simulated.

To validate our background modeling and neutron tagging procedure we also
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
the same neutron tagging cut, with a 20% signal efficiency, for all energies.
While most spallation cuts are removed in order to focus on neutron tagging,
multiple spallation cuts are conserved in order to avoid contaminating the
neutron detection window with isotope decays —an effect not accounted in the
mistag rate predictions. The predicted and observed event spectra are shown in
Fig. 20 and agree within uncertainties.

After the cuts derived in Section V.6 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 time constant of $\tau=154\pm 104~{}\mu$s,
which is compatible with the expected neutron capture time in water. In the
15.5$-$79.5 MeV region, that will be used for the spectral analysis, we also
find 248 events with zero neutrons and 9 events with more than one neutron.
The observation times and reconstructed vertices for these events are shown in
Figs. 22 and 23. No event time cluster has been observed and the time-
dependence of the event rate is consistent with the livetime variations over
the SK-IV period. Events are also uniformly distributed all over the FV,
irrespective of the number of neutrons. The space and time distributions of
the selected events are therefore consistent with what is expected for DSNB
candidates.

Figure 19: Data spectrum after noise reduction, spallation,
$d_{\mathrm{eff}}$, and all other positron candidate selection cuts. Here the
cuts are applied in the order shown above. Figure 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.
Figure 21: Time differences between reconstructed neutron captures and their
corresponding prompt events (positron candidates). Here we show the data for
one tagged neutron. The solid orange curve shows the fit result by an
exponential plus a constant, giving a time constant of $\tau=154\pm
104~{}\mu$s.

Figure 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 the data-taking rate over the SK-IV
period are taken into account.

Figure 23: Top (top) and transverse (bottom) views of reconstructed vertices
for events with one neutron (left), and with zero or more than 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.

##  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 of them, search for DSNB IBDs. In what follows, we describe how to
estimate the various backgrounds encountered in this analysis, model their
energy dependence, and estimate the corresponding uncertainties.

###  VI.1 Background estimation

####  VI.1.1 Atmospheric neutrinos

In the high end of the analysis window ($E_{\rm rec}>19.5$ MeV), the
atmospheric neutrino spectrum is dominated by decay of invisible muons or
pions from CC or NC interactions, forming the well-known Michel spectrum.
Systematic uncertainties in this region will therefore only affect the total
number of predicted atmospheric neutrino events. Here, we estimate this number
using the sideband region of $29.5<E_{\rm rec}<49.5$ MeV, that has a similar
background composition. 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 statistical fluctuations in the number of events in the sideband and is
around $19$%. The resulting spectrum is labeled as “non-NCQE” in Fig. 24.

Below $19.5$ MeV, the dominant atmospheric neutrino background arises from
NCQE interactions. The cross section measurements from T2K Abe  _et al._
(2019) 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 this measurement;
$f_{\rm\nu\mathchar 45\relax NCQE}=0.80\pm 0.08({\rm stat.})^{+{\rm
0.24}}_{-{\rm 0.18}}({\rm syst.})$ and $f_{\rm\bar{\nu}\mathchar 45\relax
NCQE}=1.11\pm 0.18({\rm stat.})^{+{\rm 0.29}}_{-{\rm 0.22}}({\rm syst.})$,
respectively. Note that this rescaling provides an inclusive estimation of the
NCQE and NC 2p2h interactions. The errors of these factors are taken as cross
section uncertainties.

Two additional flux uncertainties are considered for the NCQE background: the
atmospheric neutrino flux uncertainty and the uncertainty due to the flux
difference between the T2K beam and atmospheric neutrinos. For the former, a
$15$% uncertainty is taken from Refs. Honda  _et al._ (2007, 2011). 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 of cross section model uncertainties is different. To
estimate this uncertainty, the ratios of flux-averaged events in the T2K beam
and atmospheric neutrino fluxes for different cross section models are
calculated. Here six models are considered: Spectral Function (SF) Benhar  _et
al._ (1994, 2005), Relativistic Mean Field (RMF) Horowitz and Serot (1981);
Maieron  _et al._ (2003); Caballero  _et al._ (2005); Gonzalez-Jimenez  _et
al._ (2013), Superscaling (SuSA) Amaro  _et al._ (2005, 2006), Relativistic
Green’s Function (RGF) Capuzzi  _et al._ (1991); Gonzalez-Jimenez  _et al._
(2013); Meucci  _et al._ (2004); Meucci and Giusti (2014) with two functional
forms for the potential (EDAI, Democratic), and Relativistic Plane Wave
Impulse Approximation (RPWIA) Gonzalez-Jimenez  _et al._ (2013). 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 predicted efficiencies in the 241Am/Be
calibration. The neutron tagging efficiency will also depend on the distance
between the neutron emission and capture vertices. Indeed, neutrons produced
in atmospheric interactions can travel up to a few meters and are hence more
difficult to identify than neutrons produced in DSNB IBDs. The effect of the
uncertainties in modeling the neutron travel distance on the efficiency is of
around $7$% Wan  _et al._ (2019). 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. Here the recent measurement of the
neutron multiplicity after (mostly CC) interactions of neutrinos from the T2K
beam inside the SK tank is used Akutsu (2019). This measurement allows to
compare observations with Monte-Carlo simulation predictions for neutron
multiplicity as a function of the reconstructed squared momentum transferred
to the nucleus, $Q^{2}$, 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 $Q^{2}$ range, and renormalize
the present Monte-Carlo simulation prediction to accommodate this gap. This
procedure gives variances of approximately 40% in the number of events with
the tagged neutron being one, both for neutrinos and antineutrinos.

The NCQE spectral shape is sensitive to the $\gamma$ emission model, affecting
the current differential limit extraction. It is complicated to estimate the
associated systematic uncertainty, as $\gamma$ rays are produced at every
stage of the NCQE process, from the primary neutrino interaction to the
secondary nuclear reaction. The modeling of these $\gamma$ emission is in
particular a probable cause of the currently observed discrepancy between the
observed and predicted Cherenkov angle distributions using neutrinos from the
T2K beam Abe  _et al._ (2019). One notable effect of this mismodeling is a
smearing of the reconstructed energy analogous to a resolution effect. To
evaluate the associated systematic uncertainty we model this smearing by
convolving the original spectrum with a Gaussian distribution and using the
energy-dependent ratio between the smeared and unsmeared spectra to
renormalize the Monte-Carlo events. We then compute the number of renormalized
events in each energy bin after the reduction steps described in Section V.
The difference between this result and the nominal prediction is then taken as
the systematic uncertainty associated to the $\gamma$ emission modeling. For
this study, we set the smearing parameter $\sigma$ to 3 MeV for all energies;
this value allows to cover the discrepancy observed between the predicted and
observed $\theta_{\rm C}$ distributions in the T2K measurements mentioned
above. The resulting uncertainty ranges from 30% to 60% in the 7.5$-$29.5 MeV
range.

Finally, we combine all the systematic uncertainties described above by adding
them in quadrature. The total systematic uncertainty on the NCQE background
then varies from 60% to 80% for low to high energies. The predicted NCQE
spectrum with the T2K scaling is shown in Fig. 24.

####  VI.1.2 Spallation 9Li

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

$\displaystyle B_{{}^{9}\text{Li}}=R_{{}^{9}\text{Li}}\times
T_{\text{live}}\times 22.5~{}\text{kton}$
$\displaystyle\times\text{Br}[^{9}\text{Li}\rightarrow\beta+n]\times
f_{\text{window}}\times\epsilon_{\text{reduc}},$ (7)

where $R_{{}^{9}\text{Li}}$ is the production 9Li rate measured at SK
($0.86\pm 0.12({\rm stat.})\pm 0.15({\rm syst.})~{}{\rm
kton^{-1}}$$\cdot$${\rm day^{-1}}$) Zhang  _et al._ (2016), $T_{\text{live}}$
is the operational livetime (2075.3 and 2970.1 days below and above $E_{\rm
rec}=9.5$ MeV, respectively),
$\text{Br}[^{9}\text{Li}\rightarrow\beta+n]=0.508$ is the branching ratio for
$\beta$+$n$ decay, $f_{\text{window}}$ is the fraction of the 9Li decay energy
spectrum above the search energy threshold, and $\epsilon_{\text{reduc}}$ is
the reduction efficiency. We estimate $f_{\text{window}}$ using the Monte-
Carlo simulation elaborated for IBD events, renormalized to fit the 9Li
spectrum. This same simulation can be used to estimate the efficiencies for
the noise and positron reductions as well as the neutron tagging, that are the
same as for the DSNB signal. The spallation cut efficiency is estimated using
the procedure outlined in Section 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 in Appendix A.4, our
method assumes that the time difference between an isotope decay and its
parent muon is the only isotope-dependent spallation observable. In this
study, we account for the impact of other spallation observables by defining a
50% uncertainty on the number of 9Li events. The uncertainty on the event rate
is taken from the previous SK study Zhang  _et al._ (2016) as 22%. Other
uncertainties come from the reduction, 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%.

####  VI.1.3 Reactor neutrinos

Reactor neutrino backgrounds are estimated by renormalizing the IBD Monte-
Carlo simulation. They populates only in the lowest energy bin, as seen in
Fig. 24. In this analysis, a conservative 100% uncertainty is assigned to
their estimated rate.

####  VI.1.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 by computing the number of
events left after noise reduction, spallation, and positron reduction cuts
$N^{\text{data}}_{\text{pre-ntag}}$, and rescaling it by the neutron mistag
rate:

$\displaystyle B_{\text{acc}}=\epsilon_{\text{mis}}\times
N^{\text{data}}_{\text{pre-ntag}}.$ (8)

The associated systematic errors are highly cut-dependent and arise from the
time dependence of the PMT dark noise and natural radioactivity. In order to
evaluate its effect on the background rejection, we separate the random
trigger data mentioned in section III.2 over all the SK-IV period into $10$
time bins of about eight months each. Then we evaluate the BDT mistag rate for
each bin separately and, for a given signal efficiency, take the standard
deviation from the average mistag rate as systematic 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.

###  VI.2 Differential upper limit

The observed data and expected background spectra are shown in Fig. 24. The
corresponding p-values are calculated for each energy bin by performing pseudo
experiments as follows. We vary the number of expected background events
randomly based on Gaussian distributions with their widths being the 1$\sigma$
systematic uncertainties from each background source. We then calculate
p-values from the resulting distributions and the observed number of events in
each bin, as shown in Table 5. Here the most significant bins are found to be
at 2$\sigma$ 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 $\bar{\nu}_{e}$ flux using the
pseudo experiments above. Here we obtain the 90% C.L. upper bound on the
number of signal ($N_{90}^{\rm limit}$) as an excess of the observation over
the background expectation, by varying the number of observed events in each
energy bin by their statistical uncertainties as well as varying the number of
expected background events by their systematic uncertainties from each source.
The 90% C.L. upper limit on the $\bar{\nu}_{e}$ flux is then calculated as:

$\displaystyle\phi_{90}^{\rm limit}=\frac{N_{90}^{\rm limit}}{t\cdot
N_{p}\cdot\bar{\sigma}_{\rm IBD}\cdot\epsilon_{\rm sig}},$ (9)

where $t$ is the operational livetime [sec], $N_{p}$ is the number of free
protons in the SK’s FV, $\bar{\sigma}_{\rm IBD}$ is the IBD cross section
[${\rm 10^{-41}~{}cm^{2}}$] at a mean neutrino energy in the corresponding
region ($\bar{E}_{\nu}$), and $\epsilon_{\rm sig}$ is the signal efficiency.
Note that the neutrino energy is obtained as $E_{\nu}=E_{\rm rec}+1.8$ MeV in
IBD. For the expected sensitivity, the same procedure is applied while the
number of observed event is replaced with the number of nominal background
prediction. Here statistical uncertainties on the backgrounds are considered
instead. The expected sensitivities and observed upper limits from this work
are summarized in Table 5. 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 have been obtained using a
different methodology, more similar to the one used for the spectral analysis
in Section VII. In addition, the SK-I,II,III analysis not only used a higher
energy threshold but also did not involve neutron tagging Zhang (2012); Bays
_et al._ (2012). 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 the world’s tightest over the
search region 111The latest DSNB search result from KamLAND Abe  _et al._
(2021c), while still unpublished, serves tighter or comparable limits below
$E_{\nu}=13.3$ MeV.. The current limit disfavors the Totani$+$95 model Totani
and Sato (1995) and the most optimistic predictions of the Kaplinghat$+$00
model Kaplinghat  _et al._ (2000), and is reaching close to several other
model predictions. Due to not only a higher exposure but also higher cut
efficiencies and more precise background estimation, this new analysis
considerably improves on the previous SK-IV DSNB model-independent search
Zhang (2012).

Figure 24: Reconstructed energy spectra after data reductions for the 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 in the signal window representing the total systematic uncertainty. The red dashed line represents a DSNB signal expectation from the Horiuchi$+$09 model Horiuchi  _et al._ (2009) shown only for the signal window. Table 5: Summary on the 90% C.L. expected sensitivities and observed upper limits as well as the corresponding p-values in each electron antineutrino energy bin ($E_{\nu}=E_{\rm rec}+1.8$ MeV). $E_{\nu}$ [MeV] |  Expected [${\rm cm^{-2}sec^{-1}MeV^{-1}}$] |  Observed [${\rm cm^{-2}sec^{-1}MeV^{-1}}$] |  p-value
---|---|---|---
9.3$-$11.3 |  $4.44\times 10^{1}$ |  $3.71\times 10^{1}$ |  $0.346$
11.3$-$13.3 |  $1.14\times 10^{1}$ |  $2.04\times 10^{1}$ |  $0.886$
13.3$-$15.3 |  $4.17\times 10^{0}$ |  $9.34\times 10^{0}$ |  $0.938$
15.3$-$17.3 |  $1.87\times 10^{0}$ |  $3.29\times 10^{0}$ |  $0.830$
17.3$-$19.3 |  $8.48\times 10^{-1}$ |  $5.08\times 10^{-1}$ |  $0.243$
19.3$-$21.3 |  $4.64\times 10^{-1}$ |  $6.84\times 10^{-1}$ |  $0.686$
21.3$-$23.3 |  $3.28\times 10^{-1}$ |  $1.27\times 10^{-1}$ |  $0.073$
23.3$-$25.3 |  $2.11\times 10^{-1}$ |  $3.75\times 10^{-1}$ |  $0.597$
25.3$-$27.3 |  $2.13\times 10^{-1}$ |  $7.77\times 10^{-2}$ |  $0.051$
27.3$-$29.3 |  $1.98\times 10^{-1}$ |  $2.42\times 10^{-1}$ |  $0.605$
29.3$-$31.3 |  $1.50\times 10^{-1}$ |  $7.09\times 10^{-2}$ |  $0.126$
Figure 25: The 90% C.L. expected and observed upper limits on the
extraterrestrial electron antineutrino flux from the present work, in
comparison with previously published results from SK Bays  _et al._ (2012);
Zhang  _et al._ (2015) and KamLAND Gando  _et al._ (2012) and DSNB theoretical
predictions from Fig. 1 (in gray). The upper limit from Ref. Bays  _et al._
(2012) (blue) has been derived in Ref. Zhang  _et al._ (2015).

##  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 SK phases Bays
_et al._ (2012) to achieve a $22.5\times 5823$ kton$\cdot$day exposure. While
atmospheric neutrino backgrounds, notably from decays of invisible muons and
pions, will dominate over most of the analysis window, we also account for
possible residual spallation in 15.5$-$19.5 MeV.

###  VII.1 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 $\theta_{\rm C}\in[38,50]^{\circ}$, that will
contain most of the signal and the irreducible backgrounds, and two sidebands
with $\theta_{\rm C}\in[20,38]^{\circ}$ and $\theta_{\rm
C}\in[78,90]^{\circ}$. The low Cherenkov angle region will be populated with
mostly atmospheric backgrounds involving visible muons and pions while the
high angle region will be mostly populated by NCQE atmospheric neutrino events
with multiple $\gamma$ rays. Finally, we separate events with exactly one
identified neutron from the others, thus defining an “IBD-like” and a “non
IBD-like” region. Note that due to the low efficiencies of the neutron tagging
cuts the non IBD-like region is expected to contain a sizable 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 6.

Table 6: 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_{{\rm tagged\mathchar 45\relax}n}$ $\theta_{\rm C}$ | 20$-$38∘ | 38$-$50∘ | 78$-$90∘
---|---|---|---
1 | I | II | III
$\neq$1 | IV | V | VI

###  VII.2 Spectral shape fitting

We perform a simultaneous fitting of the signal and background spectra to the
observed data in all six regions of parameter space defined in Table 6 using
an extended maximum likelihood method. Performing this type of analysis
requires knowing the shapes of the signal and background spectra in each of
the regions. While the signal spectrum can be reliably predicted by the IBD
Monte-Carlo simulation for a given DSNB model, the treatment of the
atmospheric neutrino and spallation backgrounds is more complex. For this
study, we follow the method described in Ref. Bays  _et al._ (2012) and divide
these backgrounds into the following five categories.

##### 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 directly
from data, using a sample of electrons produced by cosmic ray muon decays (a
similar sample for atmospheric and accelerator neutrino oscillation and proton
decay analyses is described in Refs. Ashie  _et al._ (2005); Regis  _et al._
(2012)). We then use the atmospheric neutrino Monte-Carlo simulation to
compute the 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.

#####  $\nu_{e}$ CC interactions:

in this category we find backgrounds arising from CC interactions of electron
neutrinos and antineutrinos, with no visible muons and pions in the final
state. Their contributions will dominate in the signal regions II and V above
$50$ MeV. We estimate the associated spectral shapes in all regions using the
atmospheric neutrino Monte-Carlo simulation. Similarly to Michel electrons,
this background is negligible outside the regions II and V.

#####  $\mu/\pi$-producing interactions:

visible muons and pions will be associated with low Cherenkov angles, as these
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 the signal
regions. We extract the associated spectral shapes by considering an
atmospheric Monte-Carlo sample with only CC interactions, visible muons and
pions, and no electrons.

##### NCQE interactions:

unlike the other types of atmospheric neutrino events the NCQE spectrum peaks
at low energies, similarly to the DSNB spectrum. Here we define spectral
shapes using simulated NCQE events with no electrons. This background, often
involving multiple $\gamma$ rays, will dominate in the high Cherenkov angle
regions III and VI, but will also contribute to all other regions.

##### 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 in signal 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
the other two isotopes, it is expected to be subdominant. To model the
spallation spectrum and its variability we parameterize this spectrum and vary
its isotope composition as follows:

$\displaystyle S_{\rm spall}(E)=$ $\displaystyle f_{0}S_{{}^{9}{\rm C}}(E)$
(10) $\displaystyle+(1-f_{0})\left(f_{1}S_{{}^{8}{\rm
B}}(E)+(1-f_{1})S_{{}^{8}{\rm Li}}(E)\right),$ $\displaystyle\text{with
}f_{0}\in[$ $\displaystyle 0,0.5]\quad\text{and}\quad f_{1}\in[0,1]$

where $S_{i}(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, we constrain 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 the average
between the steepest and the least steep spectra and parameterize possible
deviations from this hypothesis as systematic uncertainties. We evaluate the
$\beta$ spectra using external measurements Winter  _et al._ (2006) for 8B
decays and GEANT4 for 8Li and 9C. As with 9Li, we then use these spectra to
renormalize events in the IBD Monte-Carlo simulation and thus model detector
resolution effects.

Although the first four background categories do not include all possible
types of atmospheric neutrino interactions, their associated spectra can be
linearly combined to model other background categories. Pion-producing NC
backgrounds can notably be treated as a superposition of NCQE and decay
electron backgrounds. The five categories described here hence allow to map
all the background configurations relevant to this analysis.

For each background category we define probability distribution functions
(PDFs) which model the corresponding energy spectra in each of the six
Cherenkov angle and neutron regions. These PDFs are normalized to $1$ across
all regions. We apply a similar procedure to obtain signal spectral PDFs
associated to different DSNB models, using the IBD Monte-Carlo simulation. For
$N_{\text{events}}$ observed events with energies
$\\{E_{1},...,E_{N_{\text{events}}}\\}$ we then extract the most probable
numbers of signal and background events by performing an extended unbinned
maximum likelihood fit, using the following likelihood function:

$\displaystyle\mathcal{L}(\\{N_{j}\\})$
$\displaystyle=e^{-\sum_{j=0}^{5}N_{j}}\,\prod_{i=1}^{N_{\text{events}}}\sum_{j=0}^{5}N_{j}\,\text{PDF}_{j}^{(r)}(E_{i}).$
(11)

Here, the index $j=0,...,5$ refers to the signal and the five background
categories described above. Then $\\{N_{j}\\}$ designates the numbers of
events for the signal and the different background types, across all six
Cherenkov angle and tagged neutron regions. For each event,
$\text{PDF}_{j}^{(r)}(E_{i})$ is the PDF for the signal or background category
$j$ in the region $r$ (for the regions defined in Table 6).

###  VII.3 Systematic uncertainties

One considerable advantage of using the likelihood defined in Equation 11 is
that the sizes of the different background contributions can be treated as
nuisance parameters. This likelihood, however, does not account for 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 account for systematic
uncertainties on the spallation and atmospheric neutrino background spectral
shapes, energy reconstruction, and reduction efficiencies by convolving the
likelihood introduced in Equation 11 with suitable probability distributions.
Here we neglect the systematic uncertainties associated with the $\mu/\pi$
background. The amounts of these background in the signal regions II and V are
indeed constrained to be negligible by observations from the low angle regions
I and IV, a result that is robust under large variations of the assumed
spectral shapes in these regions. We also neglect the systematic uncertainties
on the spectral shape of the background associated with invisible muon and
pion decays since this shape has been directly extracted from a large sample
of SK-IV data. We then consider five sources of systematic uncertainties:
energy scale and resolution, spectral shape predictions for the spallation,
NCQE, and $\nu_{e}$ CC backgrounds, and reduction cut efficiencies.

####  VII.3.1 Energy scale and energy resolution

To estimate the impact of energy scale and energy resolution uncertainties on
the signal and atmospheric background spectra, we follow the procedure
described in Ref. Bays  _et al._ (2012). The energy scale and energy
resolution uncertainties are obtained from the solar neutrino analysis and
increased to 2% and 3% respectively to account for the higher energies studied
here. We model the effect of these uncertainties by distorting the signal and
atmospheric background PDFs in Equation 11 using the energy resolution
functions from Refs. Hosaka  _et al._ (2006); Cravens  _et al._ (2008); Abe
_et al._ (2011, 2016). We parameterize the amount of distorsion applied to
these PDFs using the variable $\epsilon$ such that $\epsilon=1$ represents a
$1\sigma$ deviation in the energy scale and resolution. The final likelihood
can then be expressed as:

$\displaystyle\mathcal{L}^{\prime}=\int_{-\infty}^{+\infty}e^{-\epsilon^{2}/2}\mathcal{L}(\epsilon)\,\mathrm{d}\epsilon,$
(12)

where $\mathcal{L}(\epsilon)$ represents the likelihood from Equation 11 with
the PDFs distorted by an amount $\epsilon$. These uncertainties have a
negligible impact on the final result and are therefore not taken into account
in the final fit.

While energy scale and resolution uncertainties on signal and atmospheric
background spectra can be neglected in this analysis, the spallation
background spectrum is particularly steep and hence more sensitive to the
modeling of energy reconstruction. In this analysis, we will therefore account
for the associated modeling uncertainties when estimating spallation spectral
shape systematic uncertainty.

####  VII.3.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 the latter by varying the 8B, 8Li, and 9C fractions as shown in
Equation 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
$\pm$1$\sigma$ variations from the nominal spallation spectrum. Finally, we
parameterize the total spectral distorsion associated with the systematic
uncertainty by using a third order polynomial $\mathcal{P}_{3}$ such that:

$\displaystyle\text{PDF}_{\text{new}}$
$\displaystyle=\text{PDF}_{\text{old}}\times
N\left[1+\epsilon\mathcal{P}_{3}(E)\right],$ (13)

where $\epsilon$ is expressed in units of $\sigma$. More details about these
estimates are given in Appendix D. We then convolve the likelihood
$\mathcal{L}$ from Equation 11 with a Gaussian prior on $\epsilon$:

$\displaystyle\mathcal{L}^{\prime}=\int_{-2}^{2}\mathcal{L}(\epsilon)e^{-\frac{\epsilon^{2}}{2}}\,\mathrm{d}\epsilon.$
(14)

####  VII.3.3 Atmospheric $\nu_{e}$ CC spectral shapes

The shape of the $\nu_{e}$ CC background for electron energies below about
$100$ MeV has not been studied in-depth to date, although Monte-Carlo
simulations predict a linear increase in the number of events with the
electron energy. In our analysis, we model this spectral shape uncertainty by
allowing the slope of the $\nu_{e}$ CC spectrum to vary in the two signal
regions. We then integrate the distorted PDF using pseudo-Gaussian weights. A
detailed description of this procedure is provided in Ref. Bays  _et al._
(2012).

At this step of the analysis, we assume the spectral shape distortions in the
$=$1 and $\neq$1 tagged neutron regions to be fully correlated and model the
effect of neutron multiplicity uncertainties separately.

####  VII.3.4 Relative normalization for NCQE events

The Cherenkov angle distribution of atmospheric NCQE events is associated with
large uncertainties stemming from the modeling of $\gamma$ ray emission. We
incorporate this uncertainty in the current analysis by allowing the relative
normalization of the NCQE PDFs between the different Cherenkov angle regions
to vary. Following the procedure outlined in Ref. Bays  _et al._ (2012), we
define a 1$\sigma$ variation of this normalization to correspond to a $100\%$
change of the number of NCQE events in the intermediate Cherenkov angle
region. Additionally, in order for the overall normalization of the NCQE PDFs
across regions to remain constant, we reweight these PDFs as follows:

$\displaystyle\text{PDF'}^{\text{NCQE}}_{\text{med, n}}$
$\displaystyle=\text{PDF}^{\text{NCQE}}_{\text{med, n}}\times(1+\epsilon),$
(15) $\displaystyle\text{PDF'}^{\text{NCQE}}_{\text{high, n}}$
$\displaystyle=\text{PDF}^{\text{NCQE}}_{\text{high,
n}}\times(1-f_{\text{high}}\epsilon),$ (16)
$\displaystyle\text{PDF'}^{\text{NCQE}}_{\text{low, n}}$
$\displaystyle=\text{PDF}^{\text{NCQE}}_{\text{low,
n}}\times(1-f_{\text{low}}\epsilon),$ (17)

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 6 and
$n=0,1$ for the tagged neutron regions respectively. $f_{\rm low},f_{\rm
med},f_{\rm high}$ are scaling factors introduced to preserve the overall