Indication of Electron Antineutrinos from the Diffuse Supernova Neutrino Background

Indication of Electron Antineutrinos from the Diffuse Supernova Neutrino Background

(Super-Kamiokande Collaboration) $^{18}$ Institute For Interdisciplinary Research in Science and Education, ICISE, Quy Nhon, 55121, Vietnam

$^{19}$ Department of Physics, Imperial College London, London, SW7 2AZ, United Kingdom

$^{20}$ Dipartimento Interuniversitario di Fisica, INFN Sezione di Bari and Università e Politecnico di Bari, I-70125, Bari, Italy

$^{21}$ Dipartimento di Fisica, INFN Sezione di Napoli and Università di Napoli, I-80126, Napoli, Italy

$^{22}$ Dipartimento di Fisica, INFN Sezione di Padova and Università di Padova, I-35131, Padova, Italy

$^{23}$ INFN Sezione di Roma and Università di Roma “La Sapienza”, I-00185, Roma, Italy

$^{24}$ ILANCE, CNRS - University of Tokyo International Research Laboratory, Kashiwa, Chiba 277-8582, Japan

$^{25}$ Department of Physics, Keio University, Yokohama, Kanagawa, 223-8522, Japan

$^{26}$ High Energy Accelerator Research Organization (KEK), Tsukuba, Ibaraki 305-0801, Japan

$^{27}$ Department of Physics, King’s College London, London, WC2R 2LS, UK

$^{28}$ Department of Physics, Kobe University, Kobe, Hyogo 657-8501, Japan

$^{29}$ Department of Physics, Kyoto University, Kyoto, Kyoto 606-8502, Japan

$^{30}$ Department of Physics, University of Liverpool, Liverpool, L69 7ZE, United Kingdom

$^{31}$ School of Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA

$^{32}$ Department of Physics, Miyagi University of Education, Sendai, Miyagi 980-0845, Japan

$^{33}$ Institute for Space-Earth Environmental Research, Nagoya University, Nagoya, Aichi 464-8602, Japan

$^{34}$ Kobayashi-Maskawa Institute for the Origin of Particles and the Universe, Nagoya University, Nagoya, Aichi 464-8602, Japan

$^{35}$ National Centre For Nuclear Research, 02-093 Warsaw, Poland

$^{36}$ National Institute of Technology, Niihama College, Niihama, Ehime 792-8580, Japan

$^{37}$ National Institute of Technology, Numazu College, Numazu, Shizuoka 410-8501, Japan

$^{38}$ Department of Physics and Astronomy, State University of New York at Stony Brook, NY 11794-3800, USA

$^{39}$ Department of Physics, Okayama University, Okayama, Okayama 700-8530, Japan

$^{40}$ Media Communication Center, Osaka Electro-Communication University, Neyagawa, Osaka, 572-8530, Japan

$^{41}$ Department of Physics, Oxford University, Oxford, OX1 3PU, United Kingdom

$^{42}$ Rutherford Appleton Laboratory, Harwell, Oxford, OX11 0QX, UK

$^{43}$ Department of Physics and Astronomy, Seoul National University, Seoul 151-742, Korea

$^{44}$ School of Mathematical and Physical Sciences, University of Sheffield, S3 7RH, Sheffield, United Kingdom

$^{45}$ Department of Informatics in Social Welfare, Shizuoka University of Welfare, Yaizu, Shizuoka, 425-8611, Japan

$^{46}$ August Chekhovski Institute of Physics, University of Silesia in Katowice, 75 Putku Piechoty 1, 41-500 Chorzów, Poland

$^{47}$ STFC, Rutherford Appleton Laboratory, Harwell Oxford, Oxford, UK

Daresbury Laboratory, Warrington, OX11 0QX, United Kingdom

$^{48}$ Department of Physics, Sungkyunkwan University, Suwon 440-746, Korea

$^{49}$ Department of Physics, Faculty of Science, Tohoku University, Sendai, Miyagi, 980-8578, Japan

$^{50}$ Department of Physics, University of Tokyo, Bunkyo, Tokyo 113-0033, Japan

$^{51}$ 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

$^{52}$ Department of Physics, Institute of Science Tokyo, Meguro, Tokyo 152-8551, Japan

$^{53}$ Department of Physics and Astronomy, Faculty of Science and Technology, Tokyo University of Science, Noda, Chiba 278-8510, Japan

$^{54}$ Faculty of Science, University of Toyama, Toyama City, Toyama 930-8555, Japan

$^{55}$ TRIUMF, 4004 Wesbrook Mall, Vancouver, BC, V6T2A3, Canada

$^{56}$ Department of Engineering Physics, Tsinghua University, Beijing, 100084, China

$^{57}$ Faculty of Physics, University of Warsaw, Warsaw, 02-093, Poland

$^{58}$ Department of Physics, University of Warwick, Coventry, CV4 7AL, UK

$^{59}$ Department of Physics, University of Winnipeg, MB R3J 3L8, Canada

$^{60}$ Department of Physics, Yokohama National University, Yokohama, Kanagawa, 240-8501, Japan

(Dated: October 27, 2026)

Introduction

Since the historic detection of 24 neutrinos from SN 1987A [1, 2, 3], a further detection of supernova neutrinos has been a central goal of neutrino astronomy. Core-collapse supernovae are expected to release on the order of $10^{53}$ erg of energy, approximately 99% of which is in the form of neutrinos (e.g., see [4, 5, 6, 7, 8, 9, 10]). Detecting bursts of neutrinos from nearby supernovae is important, but the burst rate is only a few per century [11]. Integrating over the observable universe, however, core collapses occur a few times per second [11, 12, 13]. The diffuse supernova neutrino background (DSNB) is the steady, isotropic flux of neutrinos and antineutrinos accumulated from all such collapses since the beginning of star formation [14, 15]. A measurement of its spectrum would offer unique insight into not only the core collapse mechanism, but also the fraction of collapses that form black holes and galactic chemical evolution. Many theoretical predictions have been made [16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37] and experimental searches have been conducted in water Cherenkov and liquid scintillator detectors [38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48], but no significant excess has been observed so far.

All neutrino flavors are emitted in a core collapse, but most DSNB searches target $\bar{\nu}_e$ via inverse $\beta$ decay (IBD) on free protons ($\bar{\nu}_e + p \to e^+ + n$, $E_e = E_\nu- 1.3$ MeV) [49], since in the relevant energy range for detection, roughly $E_e = 10$–$30$ MeV, the IBD cross section is 2 orders of magnitude larger than that of the subdominant $\nu$-electron elastic scattering [50]. Within this energy range, $\beta$ decays of radioisotopes produced by cosmic-ray muon spallation and atmospheric $\nu$ interactions are dominant backgrounds. The expected DSNB signal is only a few events per year [14], making the efficient and reliable tagging of the delayed neutron in coincidence with the prompt positron an important tool to suppress some backgrounds and to enhance the sensitivity of the search [51].

In this Letter, we present a search for DSNB $\bar{\nu}_e$ using neutron tagging in Super-Kamiokande (SK) with data collected between October 2008 and March 2026, yielding a total livetime of 5003.0 days. We first perform a DSNB model-independent analysis, comparing the number of observed events with the background expectation in the energy range $13.3 < E_\nu< 31.3$ MeV. We then perform a spectral fit analysis, in which the observed spectrum is fitted with spectra of the signal and six background components over the wider range $13.3 < E_\nu< 81.3$ MeV. From the latter, our primary analysis, we report the first indication (2.6$\sigma$) of the DSNB $\bar{\nu}_e$ signal.

SK is a water Cherenkov detector, with a 22.5-kton analysis fiducial volume, located in the Kamioka mine in Gifu, Japan [52]. The data reported here span six operational phases: SK-IV, V, VI, VII, VII.5, and VIII. Neutron tagging became possible starting in SK-IV (2970.1 days) following an electronics upgrade [53]. Between SK-IV and SK-V, detector refurbishment work was carried out in preparation for dissolving gadolinium sulfate into SK’s pure water. Then SK-V (379.2 days) operated as the last pure-water phase. Starting in SK-VI (552.2 days), 13 tons of Gd$_2$(SO$_4$)$_3 \cdot8$H$_2$O was dissolved [54] to enhance the neutron tagging efficiency [51]. In SK-VII (404.0 days), an additional 26.1 tons of the Gd compound was dissolved, raising the Gd mass concentration from 0.011% to 0.033% [55], and the neutron capture fraction on Gd from 50% to 75%. In October 2023, a failure occurred in the geomagnetic compensation coils arranged around the inner detector. Accordingly, in SK-VII.5 (191.8 days), photomultiplier tubes (PMTs) affected by the failure operated with reduced light collection efficiency. Following repairs to the compensation coil system, SK-VIII (505.7 days) was started and continues today.

Methodology

The IBD interactions of the DSNB signal are simulated using the supernova event generator SKSNSim [56]. Atmospheric $\nu$ backgrounds above 100 MeV are simulated using the HKKM 2011 flux [57, 58] as input to the neutrino event generator NEUT 5.9.2 [59, 60]. Quasielastic (QE) interactions, dominant in the relevant search energy range, are described by the spectral function model [61, 62, 63], and nuclear deexcitation by NucDeEx [64, 65, 66]. The detector response is modeled by SKG4, a Geant4-based SK detector simulation [67, 68, 69, 70], using the INCL cascade model coupled with NucDeEx [71]. All Monte Carlo (MC) samples and data are reconstructed using a dedicated low-energy event fitter [72] and run through the event selection steps described below. The total energy, $E_e$, of the prompt DSNB interaction candidates is reconstructed assuming all light is emitted from a single electron or positron. The reconstructed energy, $E_{\rm rec}$, is defined as the corresponding kinetic energy of that electron or positron. For IBD candidates the relation with neutrino energy is $E_\nu\simeq E_{\rm rec} + 1.8$ MeV.

Positron candidate events above $\sim$6 MeV are collected by a dedicated trigger, whose threshold (hit PMTs in 200 ns) was 70 until September 2011 and 58 thereafter in SK-IV, 60–70 in SK-V, and 60 from SK-VI onward. The prompt event window of 35 $\mu$s after the trigger time is followed by another 500 $\mu$s where all PMT hits are recorded. The full 535 $\mu$s is sufficient for the coincidence detection of the prompt positron and delayed neutron, since the expected neutron capture times range from $\sim$200 $\mu$s in pure water to $\sim$60 $\mu$s in 0.033% Gd water. The event selection then proceeds by application of noise cuts [43, 47], atmospheric $\nu$ and spallation background cuts, and neutron tagging.

Atmospheric $\nu$ interactions form major backgrounds to the DSNB signal. Neutral-current (NC) interactions resulting in nuclear deexcitation $\gamma$ rays, accompanied by a knocked-out neutron, mimic IBD [73, 74, 75, 76]. At higher energies, charged-current (CC) interactions of $\nu_\mu$ and $\bar{\nu}_\mu$ resulting in a decay electron together with a neutron can also mimic IBD [43]. These backgrounds are separated from signal by the pattern and timing of the prompt event’s PMT hits: deexcitation $\gamma$ rays produce hit patterns more isotropic than that of a single positron, while muons decay with a lifetime of 2.2 $\mu$s, producing an additional cluster of hits within the prompt event window. Two improvements to these rejection methods were made in the present analysis.

First, a machine learning algorithm was developed to target the NC interactions [77]. Its inputs are the reconstructed Cherenkov angle ($\theta_C$) [38], which quantifies the opening angle of the Cherenkov cone; multiple scattering goodness (MSG) [39, 72, 78], which quantifies the fuzziness of the Cherenkov ring from a low-energy electron; and the maximum number of hits in a sliding 50 ns window ($n_{50}$). The first two assume a single Cherenkov cone, which does not hold for NC interactions with multiple $\gamma$ rays, and hence multiple Compton electrons. If the $\gamma$ rays go in different directions, they reconstruct with unphysical opening angles ($\gg41^\circ$) and artificially small MSG values. Even if the $\gamma$ rays all align, MSG will be lower, since each Compton electron has lower energy and therefore increased multiple Coulomb scattering. Compared with previous cuts [47, 77], this technique increases the IBD selection efficiency from 79% to 89% and halves the remaining NCQE fraction, from 10% to 5%.

To reduce CC interactions, previous analyses required no decay electron candidates within the prompt event window [43, 46, 47]. In the Gd-water phases, this rejected IBD events whose neutron captures were misidentified as decay electrons, decreasing the signal efficiency by $\sim10\%$. We instead separate decay electrons from Gd neutron captures by finding the largest number of time-of-flight corrected hits in a sliding 15 ns window after the prompt peak, and using the ratio of this hit cluster’s time to its number of PMT hits. Decay electrons appear within a few microseconds of the prompt peak and produce more hits, while Gd neutron captures occur tens of microseconds later with fewer hits.

Cosmic-ray muon spallation with oxygen nuclei produces $\beta$-unstable isotopes at a rate up to $10^6$ times the predicted DSNB signal below about 20 MeV [43]. Most decay by $\beta$ emission alone and are removed by the requirement of a delayed neutron, entering the sample only in conjunction with a falsely tagged neutron. What remains is further reduced by correlating candidate events with the timing, reconstructed track, and energy deposition profile along the track of preceding muons [39, 79]. These spallation cuts and neutron tagging were optimized simultaneously to suppress the resulting accidental coincidence background. In contrast, ${}^{9}\mathrm{Li}$ undergoes $\beta+n$ decay [80, 81, 82, 83], exactly reproducing the IBD signature. Given its 0.178 s half-life, ${}^{9}\mathrm{Li}$ is difficult to efficiently remove with current spallation cuts. We therefore estimate its surviving contribution by normalizing the $\beta$ decay spectrum given in [43] to the ${}^{9}\mathrm{Li}$ production yield of 0.86 kton$^{-1}\!\cdot$day$^{-1}$ from SK’s spallation measurement [84]. Other $\beta+n$ emitters in the relevant energy range have small yields (e.g., ${}^{11}\mathrm{Li}$) and are negligible.

Delayed neutrons are identified using a neural network [47, 85] [86] trained separately for each phase using common input variables. The selection threshold is set per energy range in each phase to suppress the accidental coincidences. It is tightest at the lowest energies, where the number of events before neutron tagging is largest, and is gradually loosened above $E_{\rm rec}=15.5$ MeV, where accidental coincidences are subdominant. The resulting probability of falsely tagging a delayed neutron is $10^{-4}$ to $10^{-3}$ in all phases at the lowest energies, rising to $10^{-1}$ in the pure-water phases at the highest. The corresponding neutron tagging efficiencies are 10–30% in the pure-water phases, 40–50% in the 0.011% Gd-water phase, and 60–70% in the 0.033% Gd-water phases.

Systematic uncertainties for the atmospheric $\nu$ backgrounds arise from the interaction rate (flux times cross section) and the final-state kinematics ($\gamma$ ray and neutron multiplicities). For CC interactions, 20% is assigned to the flux [58] and 30% to the cross section [50, 87, 88, 89]. For NC interactions, the selection efficiency depends strongly on the $\gamma$ ray multiplicity, which is uncertain owing to final-state interactions inside the nucleus and secondary interactions in the detector [90, 91]. We therefore compared data and MC event rates after all cuts except neutron tagging using an independent T2K beam sample with the same final state topology [74, 77, 92], giving a 45% NC rate uncertainty, improved from the previous 70% [43, 47]. The neutron multiplicity uncertainty was estimated for CC and NC events by shaping the MC multiplicity distribution to match the mean value [43, 75] of SK’s measurement [93], reducing the uncertainty from the previous 30–40% [43, 47, 94] to 10–30%.

Systematic uncertainties for the spallation ${}^{9}\mathrm{Li}$ background are taken on the isotope yield and the remaining rate after cuts. For the yield, 22% is assigned based on [84]. Previous analyses assigned a conservative 50% error to the remaining rate, estimated from the ${}^{9}\mathrm{Li}$ half-life alone and assuming the same spatial correlation with the parent muon as other isotopes [43]. We instead evaluate this correlation directly by simulating ${}^{9}\mathrm{Li}$ production with SKG4 using the Bertini cascade model [69]. The simulated distance distribution to the parent muon is consistent with a ${}^{9}\mathrm{Li}$-dominant data sample selected with loose neutron tagging and spallation cuts applied to events below $E_{\rm rec}=15.5$ MeV. Comparing estimated and true ${}^{9}\mathrm{Li}$ efficiencies in simulation, the difference in spatial correlation was found to be subdominant to that in temporal correlation, reducing the remaining rate error to 30%.

The neutron tagging systematic uncertainty was evaluated with an americium-beryllium calibration source [95]. The data–MC difference in tagging efficiency was found to be 10% for the pure-water phases, 3% for the 0.011% Gd-water phase, and 5% for the 0.033% Gd-water phases. Table I summarizes all systematics.

Item	SK-IV	SK-V	SK-VI	SK-VII	SK-VII.5	SK-VIII
(A) Atmospheric $\nu$ flux	$\pm20\%$	$\pm20\%$	$\pm20\%$	$\pm20\%$	$\pm20\%$	$\pm20\%$
(B) Neutrino cross section	$\pm30\%$	$\pm30\%$	$\pm30\%$	$\pm30\%$	$\pm30\%$	$\pm30\%$
(C) Atmospheric $\nu$ NCQE rate	$\pm45\%$	$\pm45\%$	$\pm45\%$	$\pm45\%$	$\pm45\%$	$\pm45\%$
(D) Neutron multiplicity	$\begin{matrix}+27\\-30\end{matrix}\%$	$\begin{matrix}+27\\-30\end{matrix}\%$	$\begin{matrix}+18\\-31\end{matrix}\%$	$\begin{matrix}+12\\-32\end{matrix}\%$	$\begin{matrix}+12\\-32\end{matrix}\%$	$\begin{matrix}+12\\-32\end{matrix}\%$
(E) Neutron tagging efficiency	$\pm10\%$	$\pm10\%$	$\pm3\%$	$\pm5\%$	$\pm5\%$	$\pm5\%$
(F) Spallation ${}^{9}\mathrm{Li}$ yield	$\pm22\%$	$\pm22\%$	$\pm22\%$	$\pm22\%$	$\pm22\%$	$\pm22\%$
(G) Spallation ${}^{9}\mathrm{Li}$ remaining rate	$\pm30\%$	$\pm30\%$	$\pm30\%$	$\pm30\%$	$\pm30\%$	$\pm30\%$
(H) Accidental coincidence rate	$\pm5\%$	$\pm5\%$	$\pm5\%$	$\pm5\%$	$\pm5\%$	$\pm5\%$

Statistical Analysis—The data were divided into six analysis regions based on the reconstructed Cherenkov angle of the positron candidate ($\theta_{\rm C}\in[20^\circ,38^\circ]$, $\theta_{\rm C}\in[38^\circ,53^\circ]$, or $\theta_{\rm C}\in[78^\circ,90^\circ]$) and the number of tagged neutrons ($N_n\ne1$ or $N_n=1$). IBD positrons mostly populate the mid-$\theta_{\rm C}$ band ($\theta_{\rm C}\in[38^\circ,53^\circ]$, the signal region) in either $N_n=1$ or $N_n\ne1$ subregions, depending on the neutron tagging efficiencies of each phase. The remaining four regions are background-dominated sidebands.

In the DSNB model-independent analysis, events in the signal region with $N_n=1$ are compared with the expected background. Data at $E_{\mathrm{rec}} \in[29.5, 79.5]$ MeV were used to scale the atmospheric $\nu$ non-NCQE event component at lower energies. The resulting energy spectra are shown in Figure 1. The background-only $p$-value is obtained by comparing the observed count against a distribution built from pseudoexperiments, in which each background category is varied within its systematic uncertainty assuming Gaussian distributions [43, 47]. Applying this to the summed counts over $E_{\mathrm{rec}} \in[11.5, 29.5]$ MeV, we observe a $1.8\sigma$ tension with the background-only hypothesis ($p$-value 0.036).

Image not extracted

The spectral fit analysis uses all six analysis regions. To observe the DSNB over as wide an energy range as possible, the $N_n = 1$ threshold was lowered from the previous analysis’s $E_{\mathrm{rec}} = 15.5$ MeV to 11.5 MeV. This adds spallation ${}^{9}\mathrm{Li}$ (whose remaining rate is now constrained by simulation validated with data) and accidental coincidences as fit components. In the $N_n \ne1$ subregion, spallation $\beta$ decays without an accompanying neutron are not suppressed by neutron tagging. Previous searches [43, 47] modeled this contribution in $E_{\mathrm{rec}} = 15.5$–21.5 MeV with a single spectrum for a mixture of ${}^{8}\mathrm{Li}$, ${}^{8}\mathrm{B}$, and ${}^{9}\mathrm{C}$, the isotopes with high endpoint energies and yields. Because the isotope fractions are only loosely constrained by SK data [84], they were varied over wide ranges and the average of the steepest and least steep resulting spectra was taken as the nominal spectrum [43].

This treatment was adequate in previous searches, where the atmospheric $\nu$ backgrounds carried large uncertainties: 70% on the NCQE rate and 30–40% on the neutron multiplicity [43, 47]. With the improved NCQE rejection and the reduced uncertainties (45% and 10–30%, respectively), the fit in this energy range becomes more sensitive to the spallation modeling. This is most problematic in the pure-water phases, where most of the signal lies in the $N_n \ne1$ subregion alongside this background. A fit to these phases with the previous threshold indicated an additional contribution at $E_{\mathrm{rec}} = 15.5$–21.5 MeV, comparable in size to the expected DSNB signal there (see Appendix B). As this is most likely residual spallation background not fully described by the current modeling, we conservatively raise the $N_n \ne1$ threshold to $E_{\mathrm{rec}} = 21.5$ MeV in all phases, above which the spallation contribution is negligible. For the reference DSNB model [21], these threshold changes reduce the expected signal by 13.9% overall.

Atmospheric $\nu$ backgrounds are divided into four categories based on their final state kinematics. CC interactions of $\nu_e$ and $\bar{\nu}_e$ form a ramp spectrum above $\sim$55 MeV (“$\nu_e$ CC”). CC interactions of $\nu_\mu$ and $\bar{\nu}_\mu$ that produce muons below Cherenkov threshold yield Michel electrons, forming a bump spectrum between $\sim$20 and 55 MeV (“Decay $e$”). Events with visible muons or pions tend to be reconstructed at smaller $\theta_{\mathrm{C}}$ (“$\mu/\pi$”), while NCQE and other NC interactions producing nuclear de-excitation $\gamma$ rays tend to be reconstructed at larger $\theta_{\mathrm{C}}$

(“NC$\gamma$”). Background events from the atmospheric $\nu$ flux below 100 MeV, which is not included in the nominal simulation, were estimated by simulating the fluxes of [96, 97] with SKSNSim [56] and applying the full event selection. They contribute a few percent of the nominal atmospheric $\nu$ background at the lowest energies, rising to 10% at $E_{\mathrm{rec}} \sim30$ MeV. Since their spectral shape resembles that of the $\nu_e$ CC component, the normalization uncertainty assigned to that component covers this contribution.

An unbinned extended maximum-likelihood fit was performed on the multi-phase dataset, with a single DSNB rate shared across all phases. Parameters describing effects correlated across phases (the atmospheric and spallation background normalizations, along with their nuisance parameters) are shared. Parameters reflecting phase-specific detector conditions (accidental coincidence rates and signal nuisance parameters) are fit per phase.

Systematic uncertainties are incorporated as nuisance parameters of three classes: normalization, migration, and shape systematics. Each background normalization is constrained by a Gaussian penalty centered on its nominal MC prediction, with a width equal to its uncertainty. The DSNB normalization is free to float, but carries an uncertainty of 2.9% for SK-IV and 3.1% for SK-V–VIII, combining the IBD cross section, livetime, fiducial volume, and event selection uncertainties in quadrature.

Migration systematics move events between the $N_n$ regions while preserving the total. All components carry the neutron tagging efficiency uncertainty. Atmospheric backgrounds additionally carry the neutron multiplicity uncertainty, with separate parameters for the CC-like and NC categories, allowing them to migrate in different directions. Since spallation ${}^{9}\mathrm{Li}$ and accidental coincidences are fit only in the $N_{\mathrm{n}}=1$ region, the $N_{\mathrm{n}}$ migration systematic instead scales their normalization. Migration of NC$\gamma$ events from high to mid $\theta_{\mathrm{C}}$ due to the $\gamma$ multiplicity uncertainty lies well within the 45% normalization uncertainty and is not included.

The only shape uncertainty is that on the spallation ${}^{9}\mathrm{Li}$ spectrum. The 30% uncertainty on the ${}^{9}\mathrm{Li}$ remaining rate is implemented as an allowed change in the spectral shape above and below $E_{\mathrm{rec}}=13.5$ MeV, preserving the total event rate. The systematics implementation for each fit component is summarized in Table II.

Component	Fitted Region	Normalization	Migration	Shape
DSNB	Mid-$\theta_C$, both $N_n$	Free-floating	(E)	—
Atmospheric $\nu_e$ CC	Mid-$\theta_C$, both $N_n$	(A), (B)	(D), (E)	—
Atmospheric Decay-$e$	Mid-$\theta_C$, both $N_n$	(A), (B)	(D), (E)	—
Atmospheric $\mu/\pi$	Low-$\theta_C$, both $N_n$	(A), (B)	(D), (E)	—
Atmospheric NC$\gamma$	Mid/High-$\theta_C$, both $N_n$	(C)	(D), (E)	—
Spallation ${}^{9}\mathrm{Li}$	Mid-$\theta_C$, $N_n = 1$	(E), (F)	—	(G)
Accidental coincidence	Mid-$\theta_C$, $N_n = 1$	(E), (H)	—	—

Sensitivity studies show that the updated systematic uncertainties and spectral fitting method, together with the lowered $N_{\mathrm{n}}=1$ threshold, more than double the median discovery potential (see Appendix C).

The best-fit configuration from a joint fit of the full 5003.0-day dataset is shown in Figure 2. Using the Horiuchi+09 6 MeV model [21] for the signal shape, the best-fit DSNB flux is $3.6^{+1.6}_{-1.5}\ \mathrm{cm}^{-2}\mathrm{s}^{-1}$ ($4.8^{+2.2}_{-2.0}$ events yr$^{-1}$) for $E_{\nu}>13.3$ MeV. The profile likelihood ratio (Figure 5) has a single well-defined minimum and disfavors the null hypothesis at $2.6\sigma$ (99.5% C.L., one-sided). Using only the Gd water phases (pure water phases), the flux is determined as $3.9^{+2.2}_{-1.9}\ \mathrm{cm}^{-2}\mathrm{s}^{-1}$ ($3.3^{+2.6}_{-2.4}\ \mathrm{cm}^{-2}\mathrm{s}^{-1}$) from a fitted event rate of $5.2^{+2.9}_{-2.6}$ events yr$^{-1}$ ($4.4^{+3.5}_{-3.2}$ events yr$^{-1}$) disfavoring the null hypothesis by $2.2\sigma$ ($1.4\sigma$). Full fit results are provided in Table III.

Image not extracted

Image not extracted

Component	Parameter	Type	Joint	SK-IV	SK-V	SK-VI	SK-VII	SK-VII.5	SK-VIII
DSNB	$N_s$	rate per year	4.8	5.6	0.0	7.9	0.0	7.6	9.5
Atmospheric $\nu_e$ CC	$N_{\nu_e\mathrm{CC}}$	scale on MC nominal	1.1	1.3	1.3	0.8	0.8	1.0	0.8
Atmospheric Decay-$e$	$N_{\mathrm{Decay\text{-}e}}$	scale on MC nominal	1.2	1.2	1.2	1.2	1.3	1.2	1.0
Atmospheric $\mu/\pi$	$N_{\mu/\pi}$	scale on MC nominal	1.2	1.3	1.0	1.0	0.9	1.1	1.1
Atmospheric NC$\gamma$	$N_{\mathrm{NC}\gamma}$	scale on MC nominal	1.1	1.0	1.2	1.1	1.3	0.9	0.9
Spallation ${}^{9}\mathrm{Li}$	$N_{{}^{9}\mathrm{Li}}$	scale on MC nominal	1.1	1.0	1.0	1.0	1.0	1.0	1.0
Accidental coincidence	$N_{\mathrm{acc}}$	rate per phase	—	0.2	0.5	0.2	0.3	0.3	0.1
Atmospheric CC $N_n$ migration	$\epsilon_{\mathrm{ntag\ CC}}$	pull in units of $\sigma$	$-0.2$	0.0	$-1.1$	$-1.3$	$-0.9$	0.0	1.1
Atmospheric NC $N_n$ migration	$\epsilon_{\mathrm{ntag\ NC}}$	pull in units of $\sigma$	0.8	0.2	0.7	0.3	0.4	0.0	0.2
DSNB $N_n$ migration	$\epsilon_{\mathrm{ntag\ sig}}$	pull in units of $\sigma$	—	0.2	0.0	0.1	0.0	0.0	0.2
DSNB signal efficiency	$\epsilon_{\mathrm{eff\ sig}}$	pull in units of $\sigma$	—	0.0	0.0	0.0	0.0	0.0	0.0
${}^{9}\mathrm{Li}$ spectral shape	$\epsilon_{{}^{9}\mathrm{Li}}$	pull in units of $\sigma$	0.1	0.1	0.0	0.0	0.0	0.0	$-0.1$

Discussion

The best-fit DSNB flux lies near the upper end of the 2.1–3.9 $\mathrm{cm}^{-2}\mathrm{s}^{-1}$ range predicted by the Horiuchi+09 6 MeV model [21]. The observed flux using this best-fit value for normalization is shown against the predicted flux in Figure 3. Repeating the fit with the Kresse+21 [27] and Ashida+23 [32] signal spectral shapes gives similar significances ($2.7\sigma$ and $2.6\sigma$, respectively). This analysis combines data from SK-IV [43] and V (pure water) as well as VI, VII [47], VII.5, and VIII (Gd water). For $E_{\nu}>17.3$ MeV, this result corresponds to $1.8^{+0.8}_{-0.7}\ \mathrm{cm}^{-2}\mathrm{s}^{-1}$, consistent within uncertainties with the $1.4^{+1.5}_{-1.2}\ \mathrm{cm}^{-2}\mathrm{s}^{-1}$ obtained from [47].

Image not extracted

Conclusion

We report the first indication of the DSNB $\bar{\nu}_{e}$ signal, based on a spectral fit of 5003.0 days of pure water (after 2008) and Gd-water data in SK. The best-fit flux disfavors the background-only hypothesis at $2.6\sigma$ for $E_{\nu}>13.3$ MeV. Increased Gd-water exposure, improved atmospheric $\nu$ background rejection, and reduced systematic uncertainties allowed the $N_{\mathrm{n}}=1$ threshold to be lowered to $E_{\nu}>13.3$ MeV, with the $N_{\mathrm{n}}\ne1$ threshold conservatively raised to $E_{\nu}>23.3$ MeV. A complementary DSNB model-independent analysis gives a $1.8\sigma$ tension with the background-only hypothesis for the energy range $13.3<E_{\nu}<31.3$ MeV.

**Acknowledgements—**We gratefully acknowledge the cooperation of the Kamioka Mining and Smelting Company. The Super-Kamiokande experiment has been built and operated from funding by the Japanese Ministry of Education, Culture, Sports, Science and Technology; the U.S. Department of Energy; and the U.S. National Science Foundation. Some of us have been supported by funds from the National Research Foundation of Korea (NRF-2009-0083526, NRF-2022R1A5A1030700, NRF-2022R1A3B1078756, RS-2025-00514948, RS-2025-23525600, RS-2026-25489835) funded by the Ministry of Science, Information and Communication Technology (ICT); the Institute for Basic Science (IBS-R016-D1); and the Ministry of Education (2018R1D1A1B07049158, 2021R1I1A1A01042256, RS-2024-00442775); the Japan Society for the Promotion of Science; the National Natural Science Foundation of China (Grants No. 12375100 and 12521007); the Spanish Ministry of Science, Universities and Innovation (grant PID2024-157541NB-I00); the Natural Sciences and Engineering Research Council (NSERC) of Canada; the Scinet and Digital Research of Alliance Canada; the National Science Centre (UMO-2022/46/E/ST2/00336) and the Ministry of Science and Higher Education (2023/WK/04), Poland; the Science and Technology Facilities Council (STFC) and Grid for Particle Physics (GridPP), UK; the European Union’s Horizon Europe Research and Innovation Programme HORIZON-MSCA-SE-2023 JEN-NIFER3 project, grant agreement no.101183137, H2020-MSCA-RISE-2019 SK2HK grant agreement no. 872549; and European Union’s Next Generation EU/PRTR grant CA3/RSUE2021-00559; the National Institute for Nuclear Physics (INFN), Italy; Japan Society for the Promotion of Science (JSPS) Grant-in-Aid no. 25K24689, 26K00707, and 26K17172.

• Corresponding author: [email protected]

† Corresponding author: [email protected]

‡ Corresponding author: [email protected]

§ Corresponding author: [email protected]

References

End Matter

Data taking history

—Figure 4 shows the accumulation of events in the final data sample over the six phases covered in our analysis. The $N_{\mathrm{n}} = 1$ rate increases after each Gd loading owing to the higher neutron tagging efficiency.

Image not extracted

$N_{\mathrm{n}} \ne1$ threshold and signal impact

—For the pure-water phases, keeping the $N_{\mathrm{n}} \ne1$ threshold at the previous value of $E_{\mathrm{rec}} = 15.5$ MeV with a spallation model based on that of previous searches [43, 47] increases the fitted DSNB rate by about its $1\sigma$ uncertainty. Because the fits with the previous and current thresholds share all data except the $N_{\mathrm{n}} \ne1$ events in $E_{\mathrm{rec}} = 15.5$–$21.5$ MeV, a shift of this size is larger than expected from statistical fluctuations alone. Since spallation decays without a neutron enter the $N_{\mathrm{n}} \ne1$ subregion regardless of the Gd concentration, we apply the higher threshold to all phases. Excluding this energy range lowers the expected signal and is the conservative choice.

For the reference DSNB model [21], including the energy-dependent signal efficiencies, lowering the $N_{\mathrm{n}} = 1$ threshold increases the livetime-averaged expected signal in that subregion by 23.3%, while raising the $N_{\mathrm{n}} \ne1$ threshold reduces it in the other by 41.0%. The net effect is a signal loss of 24.0% in the pure-water phases and 3.2% in the 0.011% Gd-water phase, and a signal gain of 10.7% in the 0.033% Gd-water phases, corresponding to a loss of 13.9% across the entire dataset. The fraction of the expected signal in the higher-purity $N_{\mathrm{n}} = 1$ subregion increases from 42% with a common threshold of $E_{\mathrm{rec}} = 15.5$ MeV to 60% with the new thresholds. Other signal models give similar values.

Sensitivity study

—Sensitivity studies were performed with pseudo-datasets corresponding to 4040 days of 0.033% Gd-water livetime, ten times the SK-VII exposure, in which the Horiuchi+09 6 MeV signal prediction was injected along with the expected backgrounds. The combination of the normalization penalty terms, the updated NCQE, spallation ${}^{9}\mathrm{Li}$, and neutron multiplicity systematics, and lowering of the analysis threshold in the $N_{\mathrm{n}} = 1$ signal subregion from $E_{\mathrm{rec}} = 15.5$ MeV to 11.5 MeV more than doubled the expected median discovery potential, with the $N_{\mathrm{n}} \ne1$ threshold held at $E_{\mathrm{rec}} = 21.5$ MeV.

Profile likelihood

—Figure 5 shows $-2\Delta\ln\mathcal{L}$ versus DSNB flux, profiled over all nuisance parameters.

Supplemental Material for:

Indication of Electron Antineutrinos from the Diffuse Supernova Neutrino Background

This Supplemental Material shows the energy spectra of Figure 1 and Figure 2 of the main text separately for the pure-water phases (SK-IV and SK-V) and the Gd-loaded phases (SK-VI, SK-VII, SK-VII.5, and SK-VIII). Figure S1 and Figure S2 correspond to Figure 1, and Figure S3 and Figure S4 to Figure 2.

Image not extracted

Image not extracted

Image not extracted

Image not extracted