normalization of the PDFs. The $\epsilon$ parameterizes the PDF distortion, in
units of $\sigma$. We finally inject these distorted PDFs into Equation 11 and
convolve the new likelihood with a weighted Gaussian prior, as with $\nu_{e}$
CC background uncertainties.

####  VII.3.5 Uncertainty on the neutron multiplicity

The fractions of signal and background events observed in the $=$1 and $\neq$1
tagged neutron regions depend on both the neutron tagging efficiency and the
true neutron multiplicity of the observed events. As shown in the 241Am/Be
calibration study, the efficiency of the neutron tagging cut for the DSNB
spectral analysis is associated with a systematic uncertainty of around $6$%.
Conversely, the actual neutron multiplicity for atmospheric neutrino
interactions is poorly known. As discussed in Section VI.1, neutron
multiplicity uncertainties of up to $\sim$40% have been measured for
reconstructed energies above $100$ MeV using the T2K neutrino beam Akutsu
(2019). In what follows, we therefore consider a $40\%$ systematic uncertainty
on the neutron tagging.

We model the effect of neutron multiplicity uncertainties by varying the
relative normalization of the PDFs in the $=$1 and $\neq$1 tagged neutron
regions for the $\nu_{e}$ CC, NCQE, and invisible muon and pion backgrounds.
As mentioned at the beginning of this section we neglect the systematic
uncertainty associated with $\mu/\pi$ backgrounds, which are associated with
negligible contributions in the signal region. As discussed above, we set the
1$\sigma$ variation of the fraction of events with one tagged neutron to a
conservative value of $40\%$. We then distort the PDFs in the $=$1 and $\neq$1
tagged neutron regions for each background category as follows:

$\displaystyle\text{PDF}_{j,\text{new}}^{\text{IBD}}$
$\displaystyle=\text{PDF}_{j,\text{old}}^{\text{IBD}}\times(1+\alpha_{j}\epsilon),$
(18) $\displaystyle\text{PDF}_{j,\text{new}}^{\text{non-IBD}}$
$\displaystyle=\text{PDF}_{j,\text{old}}^{\text{non-
IBD}}\times(1-\alpha_{j}f\epsilon),$ (19)

where $\alpha_{j}$ is the relative $1\sigma$ systematic uncertainty for
background category $j$, $\epsilon$ measures the magnitude of the distortion
in units of $\sigma$ and $f$ preserves the overall normalization of the PDFs.
Finally, we incorporate this uncertainty into the final likelihood by
integrating the initial likelihood similarly to Equation 14.

####  VII.3.6 Combining the background uncertainties

The three categories of systematic uncertainties on backgrounds described
above —spallation and $\nu_{e}$ CC spectral shapes, NCQE Cherenkov angle
distribution, and neutron multiplicity— are incorporated into Equation 11
using the modified likelihood defined in Equation 14 and the modified PDFs
defined in Ref. Bays  _et al._ (2012) as well as in Eqs 17 and 19 with
distorsion parameters $\epsilon_{\nu_{e}\text{CC}}$, $\epsilon_{\text{NCQE}}$,
and $\epsilon_{n}$ respectively. The final likelihood
$\mathcal{L}_{\text{syst}}$ can then be written as:

$\displaystyle\mathcal{L}_{\text{syst}}(\\{N_{j}\\})=$
$\displaystyle\int_{-1}^{3}\mathrm{d}\epsilon_{\nu_{e}\text{CC}}\int_{-1}^{3}\mathrm{d}\epsilon_{\text{NCQE}}\int_{-2}^{3}\mathrm{d}\epsilon_{n}\int_{-2}^{2}\mathrm{d}\epsilon_{\rm
spall}$
$\displaystyle\times\mathcal{L}(\\{N_{j}\\},\epsilon_{\nu_{e}\text{CC}},\epsilon_{\text{NCQE}},\epsilon_{n},\epsilon_{\rm
spall})$ (20) $\displaystyle\times G(\epsilon_{\nu_{e}\text{CC}})\times
G(\epsilon_{\text{NCQE}})\times\frac{e^{-\frac{\epsilon_{n}^{2}}{2}}}{\sqrt{\pi}}\times\frac{e^{-\frac{\epsilon_{\rm
spall}^{2}}{2}}}{\sqrt{\pi}},$

where
$\mathcal{L}(\\{N_{j}\\},\epsilon_{\nu_{e}\text{CC}},\epsilon_{\text{NCQE}},\epsilon_{n},\epsilon_{\rm
spall})$ is the likelihood defined in Equation 11 with distorted PDFs and $G$
is the weighted Gaussian introduced in Ref Bays  _et al._ (2012).

####  VII.3.7 Energy-independent efficiency systematic error

Maximizing the likelihood $\mathcal{L}_{\text{syst}}$ defined in Equation 20
over the numbers of events in all five background categories results in a
marginalized likelihood $\mathcal{L}_{\text{opt}}(N_{s})$ where $N_{s}$ is the
number of signal events. Using this likelihood to constrain the total number
of DSNB events that passed through the detector, $R$, requires rescaling
$N_{s}$ by the average value of the signal efficiency over all analysis
regions and their associated energy windows, $\epsilon_{\text{sig}}$. We
therefore need to incorporate the systematic uncertainty on
$\epsilon_{\text{sig}}$ into $\mathcal{L}_{\text{opt}}$, following a procedure
similar to the treatment of the background systematic uncertainty. Since the
neutron tagging efficiency only affects the relative normalization between
regions, the main sources of uncertainties on $\epsilon_{\text{sig}}$ will
arise from noise reduction and positron candidate selection cuts. These
uncertainties are summarized in Table 3. We also incorporate the uncertainties
associated to the fiducial volume cut, inverse beta decay cross section, and
livetime calculation, that have been taken from the solar neutrino analyses
and are summarized in Table 7. Adding all these uncertainties in quadrature,
we then rewrite the likelihood as:

$\displaystyle\mathcal{L}^{\prime}(R)=\int_{0}^{1}\mathcal{L}_{\text{opt}}(\epsilon
R)e^{-\frac{(\epsilon-\epsilon_{\text{sig}})^{2}}{2\sigma_{\text{sig}}^{2}}}\,\mathrm{d}\epsilon,$
(21)

where $\sigma_{\text{sig}}$ is the $1\sigma$ systematic uncertainty on
$\epsilon_{\text{sig}}$.

Table 7: Uncertainties associated with the fiducial volume cut, the Strumia-Vissani IBD cross section calculation, and the livetime calculation for SK-I to IV. The uncertainties for SK-I,II,III have been taken from the SK solar neutrino analyses Hosaka  _et al._ (2006); Cravens  _et al._ (2008); Abe  _et al._ (2011, 2016). For SK-IV the uncertainties have been obtained from the IBD Monte-Carlo simulation. SK phase | I | II | III | IV
---|---|---|---|---
Fiducial volume | 1.3% | 1.1% | 1.0% | 1.5%
Cross section | 1.0% | 1.0% | 1.0% | 1.0%
Livetime | 0.1% | 0.1% | 0.1% | 0.1%

Systematic uncertainties associated with reduction efficiencies weaken the
final limits by around $1\%$. Similarly, while the uncertainties associated
with spallation and atmospheric neutrino backgrounds are of order $100\%$,
their overall effect on the final limits is of about $5\%$. The limited impact
of spectral shape uncertainties is due to the large dominance of backgrounds
associated with invisible muon and pion decays, whose spectrum is well-known.
The final analysis will hence be largely dominated by statistical
uncertainties, which will allow to easily combine SK-IV observations with
results from previous SK phases.

###  VII.4 Combination with SK-I,II,III

The methodology described above, based on an extended maximum likelihood fit
over signal and sideband regions, is similar to the one described in Ref. Bays
_et al._ (2012), that performed a DSNB spectral analysis over the $22.5\times
2853$ kton$\cdot$days of data corresponding to the first three SK phases.
While this previous analysis did not include neutron tagging and used slightly
different background categories and reduction cuts, it is also dominated by
statistical uncertainties. Hence, the SK-I,II,III results can be readily
combined with the SK-IV observations.

In order to study a wide variety of the discrete and parameterized models
introduced in Section I, we build signal PDFs for all four phases of SK using
Monte-Carlo simulation to model detector resolution effects over most of the
energy range. For SK-IV, we renormalize the IBD Monte-Carlo simulation by the
predicted DSNB spectrum for all energies. For SK-I,II,III, we renormalize the
simulations used in Ref. Bays  _et al._ (2012) —with events generated in
$E_{e^{+}}$ in 10$-$90 MeV— and use the resolution functions from Refs. Hosaka
_et al._ (2006); Cravens  _et al._ (2008); Abe  _et al._ (2011) to extrapolate
them to lower energies. For the Galais+10 Galais  _et al._ (2010) and the
Nakazato+15 Nakazato  _et al._ (2015) models, for which DSNB fluxes are
available only up to antineutrino energies of 50 MeV, we assume the predicted
flux to be zero in the 50$-$81.3 MeV range. This assumption is not expected to
affect the spectral fit, as backgrounds from $\nu_{e}$ CC interactions largely
dominate at these energies. Finally, for each DSNB model, we fit the data
samples obtained in Ref. Bays  _et al._ (2012) for the first three SK phases.
Here, although we use the same cuts and search regions as in Ref. Bays  _et
al._ (2012), we account for possible residual spallation by modeling the
spectra and associated systematic uncertainty using the same procedure as at
SK-IV. We then combine the four fit results by summing the corresponding
likelihoods maximized over the nuisance parameters.

###  VII.5 Results

We apply the procedure described in the previous sections to the observed data
collected during phases I to IV of SK for the Horiuchi+09 DSNB model Horiuchi
_et al._ (2009). The associated spectral fits for SK-IV and for SK-I,II,III
are shown in Figs. 26 and 27 for the DSNB signal and for the spallation and
atmospheric neutrino background categories introduced in Section VII.2. Since
neutron tagging was not possible in the first three phases of SK, the
corresponding fits span only three regions, defined by the low, medium, and
high Cherenkov angle ranges defined in Table 6. For SK-IV, the low energy bump
in the predicted decay electron spectrum for IBD-like events is due to the
energy dependence of the neutron tagging cut, that strongly depends on the
background composition. As expected, the dominant backgrounds in the low and
high Cherenkov angle sidebands are backgrounds with visible muons and pions
and NCQE backgrounds, respectively. While most of the signal region is
dominated by backgrounds from invisible muons and pions, $\nu_{e}$ CC
contributions become sizable above about 50 MeV. For SK-IV, the background
rates are extremely low in the IBD-like region, that will hence be
particularly sensitive to IBD signals in spite of the low efficiency of the
neutron tagging cut.

In spite of introducing spallation backgrounds in the spectral fits, the
results shown in Fig. 27 for SK-I,II,III are highly similar to the ones shown
in Ref. Bays  _et al._ (2012). The number of spallation events for each phase
never exceed 2, within the upper limit obtained in Ref. Bays  _et al._ (2012).
Since the current analysis uses a similar procedure with more inclusive muon
selection criteria and 20$-$25% lower signal efficiencies, the spallation cut
effectiveness at SK-IV should remain on par with the ones observed at previous
SK phases. Yet, an important amount of spallation backgrounds ($\sim$8 events)
is predicted below 19.5 MeV at SK-IV. The origin of this important spallation
background will be further studied in future SK analyses with an enhanced
performance by gadolinium. Note that this effect was not observed in other SK
searches probing a similar energy range. These searches however used specific
cuts on e.g. the event timings or opening angles that we do not apply in this
analysis.

Even accounting for possible residual spallation, the number of events
observed in 17.5$-$19.5 MeV is higher than the best-fit predictions from Fig.
26. While this can be explained by statistical fluctuation, it led us to
investigating possible non-spallation sources for the excess observed over
atmospheric neutrino background predictions in 15.5$-$19.5 MeV. We first ruled
out the possibility of a strong DSNB signal, as it would be associated with a
steeply falling spectrum and hence to an unrealistically high supernova rate.
The hypothesis of an isolated astrophysical event is also strongly disfavored
since the observed events are uniformly distributed over the SK-IV period, as
discussed in Section V.7. A radioactive origin for the observed events is also
highly unlikely, given their high energies and their uniform spatial
distribution. The event positions are notably not correlated with the
locations of calibration sources. Residual spallation therefore remains the
most plausible explanation for the high number of events observed in
15.5$-$19.5 MeV at SK-IV.

The likelihoods associated to the SK-I to IV fits for the Horiuchi+09 model
are shown in Fig. 28 as a function of the DSNB rate for the reconstructed
equivalent electron kinetic energies $E_{\rm rec}>15.5$ MeV. This figure also
shows the likelihood for the combined analysis, with a total exposure of
$22.5\times 5823$ kton$\cdot$days. The associated limits on the DSNB flux are
shown in Table 8. The expected sensitivities of the SK-I,II,III analysis
change by only $\sim$3% with the inclusion of the spallation backgrounds while
the associated observed upper limits on the DSNB flux become tighter. The
expected sensitivities of the SK-IV and the combined analyses to the DSNB flux
at $90\%$ C.L. are $2.2~{}\bar{\nu}_{e}$ cm-2$\cdot$$\text{sec}^{-1}$ and
$1.5~{}\bar{\nu}_{e}$ cm-2$\cdot$$\text{sec}^{-1}$, respectively. Hence, the
predicted flux for this model, $1.9~{}\bar{\nu}_{e}$
cm-2$\cdot$$\text{sec}^{-1}$, is well within the reach of the SK-I,II,III,IV
analysis. The observed upper limits on the DSNB flux for SK-IV and for the
combined analysis are of $4.6~{}\bar{\nu}_{e}$ cm-2$\cdot$$\text{sec}^{-1}$
and $2.6~{}\bar{\nu}_{e}$ cm-2$\cdot$$\text{sec}^{-1}$, respectively. The
best-fit flux for the combined analysis is now of
$1.30^{+0.90}_{-0.85}~{}\bar{\nu}_{e}$ cm-2$\cdot$$\text{sec}^{-1}$. The
1.5$\sigma$ excess observed over the background prediction is higher than the
0.9$\sigma$ excess from the previous combined analysis Bays  _et al._ (2012).
While not statistically significant, this result is compatible with a wide
range of DSNB predictions with a flux comparable to the one of the Horiuchi+09
model Horiuchi  _et al._ (2009).

Figure 26: Best-fit signal and background spectra for SK-IV, assuming the DSNB
spectrum predicted by the Horiuchi+09 model Horiuchi  _et al._ (2009). The six
regions presented here include two signal regions and four sidebands in Table
6.

Figure 27: Best-fit signal and background spectra for SK-I (top), II (middle),
and III (bottom), assuming the DSNB spectrum predicted by the Horiuchi+09
model Horiuchi  _et al._ (2009). The three regions presented here correspond
to different range of the Cherenkov angle, as shown in Table 6. Since neutron
tagging is not possible in SK-I,II,III, only these three Cherenkov angle
regions are considered. Figure 28: Likelihoods associated with phases I to IV
of SK, as well as the combined likelihood. As mentioned in the main text, we
can neglect correlations between background systematic uncertainties when
combining all phases of SK.

##  VIII Discussion

###  VIII.1 Constraints on the DSNB parameter space

Over the two search strategies presented in this paper, the analysis described
in Section VI gives a model-independent differential limit on the
$\bar{\nu}_{e}$ flux while the model-dependent spectral analysis uses data
from all phases of SK and yields the tightest constraints on DSNB models. In
this section, we present the results of this analysis using both discrete DSNB
models —most of them based on supernova simulations— and simplified
parameterizations.

The $90\%$ C.L. expected and observed upper limits on the DSNB rate and flux
for theoretical models are shown in Fig. 29 for the combined SK-I,II,III,IV
analysis, with the corresponding results being tabulated in Table 8. With an
exposure of $22.5\times 5823$ kton$\cdot$days, the sensitivity of the combined
analysis at $90\%$ C.L. is the tightest to date and is on par with predictions
from the Ando$+$03, Horiuchi$+$09 ($T_{\nu}=6$ MeV), Galais$+$10, and
Kresse$+$21 models Ando  _et al._ (2003); Horiuchi  _et al._ (2009); Galais
_et al._ (2010); Kresse  _et al._ (2021). One common feature of these
optimistic models is that they use the highest cosmic star formation history
allowed by observations. Additionally, the Kresse+21 predictions Kresse  _et
al._ (2021) are based on neutrino emission models derived from state-of-the-
art supernova simulations, and take a wide range of progenitors into account.
The ability of this analysis to probe such realistic models makes the
prospects for future experiments particularly promising. Currently, an excess
of about $1.5\sigma$ has been observed in the data and the resulting $90\%$
C.L. limits on the DSNB flux improve on the constraints obtained in Ref. Bays
_et al._ (2012) for SK-I,II,III. These results confirm the exclusion of the
Totani$+$95 model Totani and Sato (1995) and of the most optimistic
predictions of the Kaplinghat model Kaplinghat  _et al._ (2000).

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

Another striking feature observed in Fig. 29 is the stability of the observed
and expected limits on the DSNB flux over all models. This stability results
from the limited sensitivity of the current analysis to the DSNB spectral
shape. In order to translate flux limits into constraints on astrophysical
parameters, we therefore ignore subtle effects associated with e.g. black-hole
formation or multiple classes of progenitors, and use the simple Fermi-Dirac
neutrino emission model described in Section III, as in Ref. Bays  _et al._
(2012). The associated two-dimensional $90\%$ C.L. limit on the neutrino
emission temperature and on the DSNB rate for $E_{\nu}>17.3$ MeV is shown in
Fig. 30, for the combined SK-I,II,III,IV analysis. This figure also shows the
expected rates for different effective neutrino emission temperatures. Here,
the band accounts for the current uncertainty on the cosmic star formation
rate, using the low, fiducial, and high estimates from Ref. Horiuchi  _et al._
(2009). Since we used two-dimensional exclusion contours instead of computing
limits for individual models, the limits appear weaker than the ones shown in
Fig. 29. These limits are on par with the SK-I,II,III analysis Bays  _et al._
(2012). In addition to the rate limits, whose interpretation is complicated by
the high SK analysis threshold, we also show $90\%$ C.L. exclusion contours
for the total $\bar{\nu}_{e}$ energy and the effective neutrino temperature in
Fig. 31. This figure also shows the regions of parameter space allowed by the
Kamiokande II and IMB experiments for the supernova SN1987A Jegerlehner  _et
al._ (1996). While the neutrino temperature for this supernova is particularly
low, future experiments should be able to probe at least the associated IMB
region.

Figure 30: The $90\%$ C.L. excluded region for the DSNB rate for $E_{\nu}>17.3$ MeV and the effective neutrino temperature, for the blackbody models described in Section III (red surface). To guide the eye, we also show the one-dimensional exclusion limit for individual neutrino temperatures. The blue band shows the predictions for blackbody models, for the low, fiducial, and high star formation rates used in Ref. Horiuchi  _et al._ (2009). Figure 31: The $90\%$ C.L. excluded region for the total $\bar{\nu}_{e}$ energy emitted during a supernova and the effective neutrino temperature, for the blackbody models described in Section III (red surface). To guide the eye, we also show the one-dimensional exclusion limit for individual neutrino temperatures. We also show the $90\%$ C.L. contours corresponding to the range of parameters allowed for SN1987A. Table 8: Best-fit values and the $90\%$ C.L. upper limits on the DSNB fluxes (in cm${}^{-2}\cdot$sec-1) for the theoretical models for phases SK-I to IV as well as for the combined analysis. Here the upper limits are given for $E_{\nu}>17.3$ MeV. | Best fit | 90% CL limit | Pred.
---|---|---|---
Model | SK4 | All | SK1 | SK2 | SK3 | SK4 | All |
Totani+95 Constant |  2.5${}^{+1.4}_{-1.3}$ |  1.3${}^{+0.9}_{-0.9}$ | 2.3 | 6.3 | 7.0 | 4.5 | 2.6 | 4.67
Kaplinghat+00 HMA (max) |  2.6${}^{+1.5}_{-1.3}$ |  1.3${}^{+0.9}_{-0.9}$ | 2.3 | 6.7 | 7.1 | 4.7 | 2.6 | 3.00
Horiuchi+09 6 MeV, max |  2.6${}^{+1.4}_{-1.3}$ |  1.3${}^{+0.9}_{-0.9}$ | 2.4 | 6.0 | 7.0 | 4.6 | 2.6 | 1.94
Ando+03 (updated 05) |  2.7${}^{+1.5}_{-1.4}$ |  1.4${}^{+0.9}_{-0.9}$ | 2.3 | 6.6 | 7.2 | 4.7 | 2.7 | 1.74
Kresse+21 (High, NH) |  2.7${}^{+1.5}_{-1.3}$ |  1.4${}^{+0.9}_{-0.9}$ | 2.3 | 6.7 | 7.2 | 4.7 | 2.7 | 1.57
Galais+09 (NH) |  2.5${}^{+1.4}_{-1.3}$ |  1.3${}^{+0.9}_{-0.9}$ | 2.3 | 6.3 | 7.0 | 4.5 | 2.6 | 1.56
Galais+09 (IH) |  2.6${}^{+1.4}_{-1.3}$ |  1.3${}^{+0.9}_{-0.9}$ | 2.3 | 6.4 | 7.0 | 4.5 | 2.6 | 1.50
Horiuchi+18 $\xi_{2.5}=0.1$ |  2.6${}^{+1.4}_{-1.3}$ |  1.4${}^{+0.9}_{-0.9}$ | 2.4 | 6.1 | 7.1 | 4.6 | 2.7 | 1.23
Kresse+21 (High, IH) |  2.7${}^{+1.5}_{-1.3}$ |  1.4${}^{+0.9}_{-0.9}$ | 2.3 | 6.7 | 7.1 | 4.7 | 2.7 | 1.21
Kresse+21 (Fid, NH) |  2.7${}^{+1.5}_{-1.3}$ |  1.4${}^{+0.9}_{-0.9}$ | 2.3 | 6.8 | 7.2 | 4.7 | 2.7 | 1.20
Kresse+21 (Fid, IH) |  2.7${}^{+1.5}_{-1.3}$ |  1.4${}^{+0.9}_{-0.9}$ | 2.3 | 6.8 | 7.2 | 4.7 | 2.7 | 1.02
Kresse+21 (Low, NH) |  2.7${}^{+1.5}_{-1.4}$ |  1.4${}^{+0.9}_{-0.9}$ | 2.3 | 6.8 | 7.2 | 4.8 | 2.7 | 0.96
Tabrizi+21 (NH) |  2.7${}^{+1.5}_{-1.3}$ |  1.4${}^{+0.9}_{-0.9}$ | 2.4 | 6.6 | 7.1 | 4.7 | 2.7 | 0.92
Kresse+21 (Low, IH) |  2.7${}^{+1.5}_{-1.4}$ |  1.4${}^{+0.9}_{-0.9}$ | 2.3 | 6.8 | 7.2 | 4.8 | 2.7 | 0.84
Lunardini09 Failed SN |  2.8${}^{+1.5}_{-1.4}$ |  1.4${}^{+0.9}_{-0.9}$ | 2.4 | 6.8 | 7.3 | 4.8 | 2.8 | 0.73
Hartmann+97 CE |  2.6${}^{+1.4}_{-1.3}$ |  1.3${}^{+0.9}_{-0.9}$ | 2.3 | 6.5 | 7.1 | 4.6 | 2.6 | 0.63
Nakazato+15 (max, IH) |  2.7${}^{+1.5}_{-1.4}$ |  1.4${}^{+1.0}_{-0.9}$ | 2.4 | 6.5 | 7.2 | 4.8 | 2.7 | 0.53
Horiuchi+18 $\xi_{2.5}=0.5$ |  2.7${}^{+1.5}_{-1.4}$ |  1.3${}^{+0.9}_{-0.9}$ | 2.2 | 7.1 | 7.1 | 4.8 | 2.6 | 0.55
Horiuchi+21 |  2.1${}^{+1.3}_{-1.2}$ |  1.2${}^{+0.9}_{-0.9}$ | 3.4 | 4.3 | 5.9 | 3.9 | 2.5 | 0.28
Malaney97 CGI |  2.7${}^{+1.5}_{-1.3}$ |  1.3${}^{+0.9}_{-0.9}$ | 2.3 | 6.8 | 7.1 | 4.7 | 2.6 | 0.26
Nakazato+15 (min, NH) |  2.8${}^{+1.5}_{-1.4}$ |  1.4${}^{+1.0}_{-0.9}$ | 2.3 | 6.8 | 7.2 | 4.8 | 2.7 | 0.19
Table 9: Best-fit values and the $90\%$ C.L. upper limits on the DSNB rates (in events$\cdot$year-1) for the theoretical models for phases SK-I to IV as well as for the combined analysis. Here the upper limits are given for $E_{\nu}>17.3$ MeV. | Best fit | 90% CL limit | Pred.
---|---|---|---
Model | SK4 | All | SK1 | SK2 | SK3 | SK4 | All |
Totani+95 Constant |  4.5${}^{+2.5}_{-2.2}$ |  2.3${}^{+1.6}_{-1.5}$ | 4.2 | 11.2 | 12.2 | 7.9 | 4.6 | 8.35
Kaplinghat+00 HMA (max) |  4.4${}^{+2.4}_{-2.2}$ |  2.2${}^{+1.5}_{-1.4}$ | 3.9 | 11.2 | 11.7 | 7.7 | 4.4 | 5.09
Horiuchi+09 6 MeV, max |  4.6${}^{+2.5}_{-2.3}$ |  2.4${}^{+1.7}_{-1.6}$ | 4.4 | 11.7 | 12.5 | 8.2 | 4.8 | 3.54
Ando+03 (updated 05) |  4.7${}^{+2.5}_{-2.4}$ |  2.4${}^{+1.6}_{-1.6}$ | 4.2 | 11.8 | 12.4 | 8.2 | 4.7 | 3.09
Kresse+21 (High, NH) |  4.5${}^{+2.5}_{-2.3}$ |  2.4${}^{+1.6}_{-1.5}$ | 4.1 | 11.7 | 12.2 | 8.1 | 4.6 | 2.76
Galais+09 (NH) |  4.4${}^{+2.4}_{-2.2}$ |  2.2${}^{+1.6}_{-1.4}$ | 4.0 | 11.0 | 11.9 | 7.8 | 4.5 | 2.74
Galais+09 (IH) |  4.4${}^{+2.4}_{-2.2}$ |  2.2${}^{+1.6}_{-1.4}$ | 4.0 | 11.0 | 11.9 | 7.7 | 4.5 | 2.62
Horiuchi+18 $\xi_{2.5}=0.1$ |  4.8${}^{+2.6}_{-2.4}$ |  2.5${}^{+1.7}_{-1.6}$ | 4.5 | 12.1 | 12.8 | 8.4 | 4.9 | 2.29
Kresse+21 (High, IH) |  4.6${}^{+2.5}_{-2.3}$ |  2.4${}^{+1.6}_{-1.5}$ | 4.1 | 11.7 | 12.2 | 8.1 | 4.6 | 2.14
Kresse+21 (Fid, NH) |  4.5${}^{+2.5}_{-2.2}$ |  2.3${}^{+1.5}_{-1.5}$ | 4.0 | 11.7 | 11.9 | 7.9 | 4.5 | 2.06
Kresse+21 (Fid, IH) |  4.5${}^{+2.5}_{-2.2}$ |  2.3${}^{+1.5}_{-1.5}$ | 4.0 | 11.6 | 11.9 | 7.9 | 4.5 | 1.75
Kresse+21 (Low, NH) |  4.5${}^{+2.5}_{-2.2}$ |  2.2${}^{+1.6}_{-1.4}$ | 3.9 | 11.6 | 11.8 | 7.9 | 4.5 | 1.65
Tabrizi+21 (NH) |  4.6${}^{+2.5}_{-2.3}$ |  2.4${}^{+1.6}_{-1.5}$ | 4.2 | 11.8 | 12.3 | 8.2 | 4.7 | 1.64
Kresse+21 (Low, IH) |  4.5${}^{+2.5}_{-2.2}$ |  2.2${}^{+1.6}_{-1.4}$ | 4.0 | 11.6 | 11.9 | 7.9 | 4.5 | 1.43
Lunardini09 Failed SN |  5.0${}^{+2.7}_{-2.5}$ |  2.6${}^{+1.7}_{-1.7}$ | 4.5 | 12.7 | 13.1 | 8.8 | 5.1 | 1.36
Hartmann+97 CE |  4.4${}^{+2.4}_{-2.2}$ |  2.2${}^{+1.6}_{-1.4}$ | 4.0 | 11.2 | 11.9 | 7.7 | 4.5 | 1.09
Nakazato+15 (max, IH) |  5.0${}^{+2.7}_{-2.5}$ |  2.5${}^{+1.8}_{-1.6}$ | 4.5 | 12.1 | 13.0 | 8.7 | 5.0 | 0.99
Horiuchi+18 $\xi_{2.5}=0.5$ |  4.2${}^{+2.3}_{-2.1}$ |  2.1${}^{+1.5}_{-1.4}$ | 3.5 | 11.4 | 10.9 | 7.4 | 4.1 | 0.87
Horiuchi+21 |  3.7${}^{+2.2}_{-2.1}$ |  2.0${}^{+1.6}_{-1.6}$ | 6.0 | 7.5 | 10.2 | 6.8 | 4.2 | 0.49
Malaney97 CGI |  4.3${}^{+2.4}_{-2.1}$ |  2.2${}^{+1.5}_{-1.4}$ | 3.8 | 11.2 | 11.3 | 7.6 | 4.3 | 0.44
Nakazato+15 (min, NH) |  4.7${}^{+2.5}_{-2.4}$ |  2.4${}^{+1.6}_{-1.6}$ | 4.0 | 11.8 | 12.1 | 8.2 | 4.6 | 0.34

###  VIII.2 Prospects at future experiments

Recent doping of SK with gadolinium (SK-Gd) as well as the advent of Hyper-
Kamiokande (HK) in the more distant future open new promising perspectives for
the DSNB search. Using gadolinium will notably allow to considerably improve
neutron tagging and hence eliminate most spallation backgrounds and lower the
threshold of the analysis. As shown in Fig. 24, however, multiple backgrounds
involving neutrons can still dominate over the DSNB up to about $15.5$ MeV.
Reactor antineutrinos, in particular, provide an irreducible background up to
$9.5$ MeV. Up to about $11.5$ MeV, 9Li decays might then dominate, and
spallation reduction will remain a key component of future analyses. Locating
muon-induced showers using neutrons will become especially powerful at SK-Gd.
After eliminating spallation, atmospheric NCQE interactions will dominate over
a wide energy range and will call for new dedicated reduction techniques. At
SK-Gd, notably, the high visibility of the neutron capture signal will allow
to locate neutron capture without using the positron vertex, which could allow
to use the neutron travel distance as a discriminating observable. Finally,
while HK may not be doped with Gd, this considerably larger detector coupled
with the accelerator neutrino beams will allow to improve the characterization
of the NCQE rate, spectral shape, and neutron multiplicity, that are currently
associated with large systematic uncertainties Abe  _et al._ (2019). In the
future, further studies using the upcoming Intermediate Water Cherenkov
Detector (IWCD) will allow to considerably reduce neutron multiplicity
uncertainties by studying monochromatic beams, and will hence be a key piece
of the DSNB search program at HK. The DSNB study will also be extended by
searches at lower energies, achieved at future large-scale liquid
scintillation detectors, such as JUNO Sisti (2020), giving deeper insights
into spectrum shape over the wide energy range in combination with HK.

##  IX Conclusion

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

## Acknowledgment

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 (KNRC) funded by the Ministry of Science,
ICT, and Future Planning and the Ministry of Education (2018R1D1A3B07050696,
2018R1D1A1B07049158), the Japan Society for the Promotion of Science, the
National Natural Science Foundation of China under Grants No. 11620101004, the
Spanish Ministry of Science, Universities and Innovation (grant
PGC2018-099388-B-I00), the Natural Sciences and Engineering Research Council
(NSERC) of Canada, the Scinet and Westgrid consortia of Compute Canada, the
National Science Centre, Poland (2015/18/E/ST2/00758), the Science and
Technology Facilities Council (STFC) and GridPPP, UK, the European Union’s
Horizon 2020 Research and Innovation Programmeunder the Marie Sklodowska-Curie
grant agreement no. 754496, H2020-MSCA-RISE-2018 JENNIFER2 grant agreement
no.822070, and H2020-MSCA-RISE-2019 SK2HK grant agreement no. 872549.

##  Appendix A Spallation cut

###  A.1 Neutron cloud cut

For each cloud, we define a specific coordinate system whose origin is the
projection of the cloud barycenter on the muon track. For a given DSNB
candidate reconstructed vertex, $\ell_{l}$ is defined as the distance to this
origin along the muon track while $\ell_{t}$ is the transverse distance to the
track. For a given cloud multiplicity, a DSNB candidate observed less than
$dt$ after the associated muon, and within the ellipse defined by $\ell_{l}$
and $\ell_{t}$, is classified as spallation background. Conditions for the
neutron cloud cut are summarized in Table 10 for different multiplicities of
the cloud.

Table 10: Conditions for the neutron cloud cut as a function of the cloud multiplicity. A number followed by the ‘$+$’ symbol represents a lower bound on this multiplicity. Multiplicity | 2$+$ | 2$+$ | 2 | 3 | 4$-$5 | 6$-$9 | 10$+$
---|---|---|---|---|---|---|---
$dt$ [sec] | 0.1 | 1 | 30 | 60 | 60 | 60 | 60
$\ell_{l}$ [cm] | 1200 | 800 | 383 | 548 | 603 | 712 | 766
$\ell_{t}$ [cm] | 1200 | 800 | 219 | 268 | 379 | 490 | 548

###  A.2 Rectangular cuts

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

Table 11: Summary on the rectangular spallation cuts for events with exactly one tagged neutron. When a given cut has requirements on both $dt$ and $\ell_{t}$, at least one should be fulfilled. Muon type | Energy range |  $dt$ requirement |  $\ell_{t}$ requirement
---|---|---|---
All |  $E_{\rm rec}<23.5$ MeV |  $dt>0.001$ sec | -
$E_{\rm rec}<23.5$ MeV |  $dt>0.1$ sec |  $\ell_{t}>400$ cm
Misfit |  $15.5<E_{\rm rec}<19.5$ MeV |  $dt>1.5$ sec | -
Well-fitted single through-going |  $15.5<E_{\rm rec}<17.5$ MeV |  $dt>7$ sec |  $\ell_{t}>150$ cm
Stopping |  $15.5<E_{\rm rec}<19.5$ MeV |  $dt>0.05$ sec | -
Poorly-fitted stopping |  $15.5<E_{\rm rec}<17.5$ MeV |  $dt>6$ sec | -
All multiple |  $15.5<E_{\rm rec}<19.5$ MeV |  $dt>0.05$ sec | -
Table 12: Summary on the rectangular spallation cuts for events with tagged neutrons different from one. When a given cut has requirements on both $dt$ and $\ell_{t}$, at least one should be fulfilled. Here tighter rectangular cuts are applied to the muon categories generating most of the spallation —single through-going and multiple muons— to remove contributions from short-lived isotopes. Muon type | Energy range |  $dt$ requirement |  $\ell_{t}$ requirement
---|---|---|---
All |  $E_{\rm rec}<23.5$ MeV |  $dt>0.001$ sec | -
$E_{\rm rec}<23.5$ MeV |  $dt>0.1$ sec |  $\ell_{t}>400$ cm
Misfit |  $15.5<E_{\rm rec}<19.5$ MeV |  $dt>1.5$ sec | -
Well-fitted single through-going |  $15.5<E_{\rm rec}<17.5$ MeV |  $dt>7$ sec |  $\ell_{t}>150$ cm
$17.5<E_{\rm rec}<19.5$ MeV |  $dt>7$ sec |  $\ell_{t}>100$ cm
Stopping |  $15.5<E_{\rm rec}<19.5$ MeV |  $dt>0.05$ sec | -
All multiple |  $15.5<E_{\rm rec}<19.5$ MeV |  $dt>0.05$ sec | -
Well-fitted multiple |  $15.5<E_{\rm rec}<19.5$ MeV |  $dt>0.1$ sec | -
$15.5<E_{\rm rec}<19.5$ MeV |  $dt>1$ sec |  $\ell_{t}>400$ cm

###  A.3 Spallation PDFs and log-likelihood ratio

When producing the spallation and random PDFs, we separate the $dt$ and
$\ell_{t}$ bins to account for the possible correlations, 0$-$0.05,
0.05$-$0.5, and 0.5$-$30 sec for $dt$ (determined by the isotope half-lives in
Fig. 3), and 0$-$300, 300$-$1000, 1000$-$5000 cm for $\ell_{t}$. Fig. 32 shows
PDFs for each spallation variable from single through-going muons. Here the
results from the most spallation-rich region ($0<dt<0.05$ sec and
$0<\ell_{t}<300$ cm) are shown. We calculate the spallation log-likelihood
ratio with the procedure described in Section V.2. An example distribution is
shown for single through-going muons in Fig. 33. The cut criteria are tuned
depending on the reconstructed energy as well as the $dt$ and $\ell_{t}$ bins.
We use the common binning for $dt$, 0$-$0.05, 0.05$-$0.5, and 0.5$-$30 sec,
independent of the muon type. In contrast, for $\ell_{t}$, we separate bins
only for most impactful muon types, that is, single through-going and multiple
muons: 0$-$200, 200$-$300, 300$-$500, 500$-$1000, and $>$1000 cm for single
through-going, and 0$-$100, 100$-$200, 200$-$300, 300$-$500, 500$-$700,
700$-$1000, 1000$-$2000, 2000$-$3000, $>$3000 cm for multiple muons.

Figure 32: Spallation and random PDFs for $dt$ (top left), $\ell_{t}$ (top
right), $\ell_{l}$ (middle left), $Q_{\mu}$ (middle right), and $Q_{\rm res}$
(bottom), for the pre (blue) and post (red) samples. Results from the single
through-going muons for the $dt$ region of 0$-$0.05 sec and the $\ell_{t}$
region of 0$-$300 cm are shown here. Figure 33: Spallation log-likelihood
ratio of the single through-going muon for $7.5<E_{\rm rec}<9.5$ MeV. The
values smaller or larger than $-$30 or 30 are filled in the underflow or
overflow bins, respectively.

###  A.4 Spallation cut performance estimation

In this section we present our methodology to estimate the impact of
spallation cuts on non-spallation “signal” events (DSNB, solar, reactor, and
atmospheric neutrino interactions), as well as on spallation backgrounds. In
what follows, we treat 9Li decays separately as they provide the dominant
$\beta$$+$$n$ signature associated with spallation.

As explained in Section V.2, random event efficiencies for multiple spallation
and neutron cloud cuts are estimated to be $98$% and $99.8$%, respectively. We
estimate the signal efficiencies for the rectangular and likelihood cuts
described above by evaluating their impact on the post-sample, as described in
the previous section. The associated efficiency ($\epsilon_{\rm random}$) can
thus be expressed as:

$\displaystyle\epsilon_{\rm random}=N_{\rm post,after}/N_{\rm before},$ (22)

where $N_{\rm before}$ and $N_{\rm post,after}$ are the numbers of events in
the post-sample before and after spallation cuts. We validate this efficiency
by evaluating the impact of spallation cuts on events from solar neutrino
elastic interactions, using the opening angle $\theta_{\rm sun}$ defined in
Section IV as a discriminator.

To estimate the background rejection power of the spallation cuts, we use a
sample of DSNB candidate events verifying $\cos\theta_{\rm sun}<0$. In order
to suppress remaining radioactivity and atmospheric neutrino events, we also
apply cuts on the effective distance to the wall and the Cherenkov angle, that
will be defined in Section V.3, and require $E_{\rm rec}<15.5$ MeV. The
spallation remaining rate ($\epsilon_{\rm spall}$) is then estimated as the
fraction of events remaining after spallation cuts, after the expected number
of remaining atmospheric neutrino events from Monte-Carlo simulation has been
subtracted.

The procedure outlined above allows to accurately estimate the spallation
remaining rate for $E_{\rm rec}<$ 15.5 MeV —where spallation largely dominates
even after cuts, but not to assess the impact of spallation cuts on individual
isotopes. Using the spallation-dominated sample described above, we therefore
refine our approach as follows. To compute the remaining fraction of a
specific isotope, we attribute a random $dt$, drawn from the isotope decay
time distribution, to each pair in the sample, leaving all other observables
untouched. This is equivalent to renormalizing the spallation likelihood from
Equation 3 by the isotope’s exponential decay time distribution. Then, we
apply the likelihood cuts described above to the pairs with the modified $dt$,
and compute the fraction of remaining events, $\epsilon_{\mathrm{pre,iso}}$.
For each event, this fraction can be expressed as:

$\displaystyle\epsilon_{\mathrm{pre,iso}}=\epsilon_{\mathrm{iso,likeli}}^{\rm(pair)}\times\epsilon_{\mathrm{random,iso}}$
(23)

where $\epsilon_{\mathrm{iso,likeli}}^{\rm(pair)}$ is the probability for a
pair formed by a spallation event and its parent muon to pass the likelihood
cuts. Conversely, $\epsilon_{\mathrm{random,iso}}$ is the probability for all
other (uncorrelated) pairs between the event and preceding muons to pass these
cuts. This probability can be computed using a post-sample analogous to the
one described in Section V.2.2 and, since we reassigned $dt$ for all pairs in
each event, will be artificially lower than the signal efficiency
$\epsilon_{\mathrm{random}}$. The probability for a given isotope decay to
pass the likelihood cuts ($\epsilon_{\mathrm{iso,likeli}}$) is thus given by:

$\displaystyle\epsilon_{\mathrm{iso,likeli}}=\epsilon_{\mathrm{pre,iso}}/\epsilon_{\mathrm{random,iso}}\times\epsilon_{\mathrm{random}}.$
(24)

The fraction of remaining events after both preselection and likelihood cuts
($\epsilon_{\mathrm{iso,spall}}$) can then be estimated as:

$\displaystyle\epsilon_{\mathrm{iso,spall}}=\epsilon_{\mathrm{pre,iso}}/\epsilon_{\mathrm{random,iso}}\times\epsilon_{\mathrm{random}}\times\epsilon_{\mathrm{presel}},$
(25) $\displaystyle\ \ \ \ \ \ \ \ \ \
\mathrm{where}\quad\epsilon_{\mathrm{presel}}=\frac{\epsilon_{\mathrm{spall}}}{\epsilon_{\mathrm{likeli}}}.$
(26)

Here, $\epsilon_{\mathrm{spall}}$ is the spallation remaining rate for all
isotopes after spallation cuts and $\epsilon_{\mathrm{likeli}}$ is the rate
after likelihood cuts only. This ratio yields $\epsilon_{\mathrm{presel}}\sim
70$% while the preselection cut alone removes about half of the background,
highlighting the significant overlap between preselection and likelihood cuts.

This procedure is used throughout this study to estimate the remaining
fractions of various isotopes. For $E_{\rm rec}<$15.5 MeV, we will notably use
it to estimate the impact of spallation cuts on the 9Li background. Above 15.5
MeV, we will estimate the remaining fractions of 8Li and 8B, two isotopes that
are particularly difficult to remove due to their long half-lives. This
estimate will allow us to model the spallation background spectrum for the
spectral analysis.

In order to estimate the remaining amounts of individual isotopes without
relying on a simulation we made two major assumptions. First, that the time
difference between spallation events and their parent muon is not correlated
with other spallation observables. Second, that the impact of spallation cuts
above 15.5 MeV can be estimated using a sample of lower energy events. Both
assumptions would be verified if the $\ell_{t}$, $\ell_{l}$, $Q_{\mu}$, and
$Q_{res}$ spallation observables described in Section V.2.2 were similarly
distributed for all isotopes. This is however not the case as their values
will depend on the isotope production mode. To estimate the effect of this
isotope dependence, we split the sample described above into four 2-MeV bins
from 7.5 to 15.5 MeV, that will have different isotope compositions. For a
given spallation cut, we then compute the spallation remaining rates for the
four subsamples and take the largest deviation from the average value as our
systematic uncertainty. This uncertainty depends on the isotope considered but
typically lies around 50%.

##  Appendix B 241Am/Be calibration study for neutron tagging

As described in Section V.4.3, in order to validate the efficiency of the
neutron tagging procedure on real data, we applied neutron tagging to data
collected in the presence an 241Am/Be source, which, surrounded by a BGO
scintillator, produces a prompt signal and a delayed neutron capture signal as
part of its radioactive decay chain. However, the data collected with the
241Am/Be source is not fully modeled by the simulation software; notably, no
modeling of the BGO scintillator is available. The Monte-Carlo simulation that
we use to model 241Am/Be events is hence similar to the IBD simulation.
Because of this limitation, a direct comparison between the neutron tagging
efficiency in simulation and calibration data would be overly pessimistic, as
scintillation in the BGO can produce prompt events that are not accompanied by
the production of a neutron, thus artificially reducing the observed
preselection efficiency of the neutron candidates (see Section V.4.1). To
account for this, we treated uncertainty associated with the preselection
separately from that of the BDT selection, following the procedure described
in Ref. Abe  _et al._ (2021a).

The preselection efficiency is based on the number of hits in a given time
window and is hence particularly sensitive to the PMT photon detection
efficiency (quantum efficiency, QE). Hence, we determine the associated
systematic uncertainty by determining the value of QE for which the predicted
$N_{10}$ distribution fits the observations best. The $N_{10}$ predicted and
observed distributions are shown in Fig. 34 for the best-fit QE, that is found
to be 2.1% larger than the QE used in the IBD simulation from Section III.2.
The preselection efficiency associated to this best-fit QE is 3.0% larger than
the efficiency predicted by our simulation. Moreover, the 1$\sigma$
uncertainty on this best-fit value, found using a $\chi^{2}$ fit, leads to a
2.1% uncertainty on the value of the preselection efficiency. The combined
preselection uncertainty is hence of 3.7%.

After the uncertainty on the preselection efficiency $\epsilon_{\rm IS}$, we
evaluate the uncertainty associated with the BDT selection cuts by computing
the relative efficiencies $\epsilon_{\rm rel}=\epsilon_{\rm
BDT+IS}/\epsilon_{IS}$ for the data and the Monte-Carlo simulation for
different values of the BDT discriminant. We determine these efficiencies by
fitting the timing distributions of neutron candidates by an exponential plus
a constant. Results for different BDT cuts are shown in Table 13 for the
different 241Am/Be datasets considered here. The maximum discrepancy between
data and Monte-Carlo lies around 12% irrespective of the BDT cut.

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

Table 13: Observed (‘Data’) and predicted (‘MC’) relative efficiencies (in %) for different BDT cuts and for the 241Am/Be calibration runs considered in this study. Data was taken in 2009 and 2016, at three different locations in the tank: in the center of the detector, near the barrel wall (Y12) and the top wall (Z15). BDT cut | Position | Year | Data | MC | (Data $-$ MC)/Data
---|---|---|---|---|---
0.40 | Center | 2009 | 69.4 | 76.2 |  $-$9.8
2016 | 80.0 | 72.3 | 9.5
Y12 | 2009 | 74.6 | 77.4 |  $-$3.7
2016 | 83.8 | 77.6 | 7.2
Z15 | 2009 | 87.6 | 77.1 | 12.0
2016 | 80.2 | 74.1 |  $-$7.7
0.90 | Center | 2009 | 50.0 | 55.3 |  $-$10.6
2016 | 56.1 | 52.8 | 5.8
Y12 | 2009 | 56.4 | 59.6 |  $-$5.7
2016 | 60.8 | 59.2 | 2.6
Z15 | 2009 | 65.8 | 58.3 | 11.4
2016 | 57.0 | 55.4 | 2.8
0.99 | Center | 2009 | 33.6 | 35.5 |  $-$5.8
2016 | 34.9 | 33.2 | 4.8
Y12 | 2009 | 36.9 | 37.8 |  $-$2.4
2016 | 38.2 | 37.1 | 2.8
Z15 | 2009 | 38.8 | 34.5 | 10.5
2016 | 31.9 | 32.2 |  $-$0.1
Figure 34: $N_{10}$ observed and expected distributions for neutrons. The
expected distribution corresponds to the PMT photon detection efficiency for
which the $N_{10}$ fits the observations best. The observed distribution has
been obtained by subtracting the $N_{10}$ background distribution from the
total observed distribution. This background distribution has been obtained
from data, using random triggers. The lower predicted number of events for
$N_{10}=6$ is due to harsher continuous dark noise reduction cuts, as
mentioned in Section V.4.1.

##  Appendix C Neutrons produced from atmospheric neutrino interactions

The neutron multiplicities predicted by the Monte-Carlo simulation for
atmospheric neutrinos are shown in Fig. 35. Note that a significant number of
atmospheric neutrino interactions produce one or more neutrons. However, these
neutrons are more energetic than for the IBDs of DSNB neutrinos and sometimes
travel further away from their production vertex. Since the current neutron
tagging algorithm requires neutrons to be close to their associated positron
vertex, it can also be used to reduce neutron-producing atmospheric
backgrounds. A detailed study of neutrons associated with atmospheric neutrino
backgrounds has been performed in Ref. Abe  _et al._ (2021a).

Figure 35: Predicted neutron multiplicity for atmospheric neutrino
backgrounds. Here we use the atmospheric neutrino Monte-Carlo simulation in
the 7.5$-$79.5 MeV range, after the noise reduction cuts.

##  Appendix D Spectral fitting

###  D.1 Modeling of atmospheric neutrino backgrounds

The probability distribution functions of the atmospheric neutrino background
spectra used in the spectral analysis are displayed in Fig. 36. As mentioned
in Section VII, we divide atmospheric neutrino backgrounds into four
categories: $\mu/\pi$, $\nu_{e}$ CC, decay electrons, and NCQE interactions.
The decay electron spectrum is obtained from data while the other spectra are
obtained from the atmospheric neutrino Monte-Carlo simulation.

In the extended maximum likelihood fit described in Section VII, atmospheric
background fluxes are treated as nuisance parameters. In Table 14, we compare
the predicted numbers of events in the four atmospheric background categories
to the best-fit values obtained when considering the Horiuchi+09 model
Horiuchi  _et al._ (2009). Since, as discussed in Section VII, atmospheric
events from e.g. pion-producing NC interactions could be categorized as either
NCQE or $\mu/\pi$, we quote ranges rather than numbers for these two
categories. For decay electrons and $\nu_{e}$ CC interactions, the predicted
and observed numbers of events are within at most 23% of each other.

Table 14: Predicted and best-fit numbers of atmospheric events using the Horiuchi+09 model Horiuchi  _et al._ (2009). The predicted numbers of events in the NCQE and $\mu/\pi$ categories are quoted as ranges, since certain types of NC events could belong to either category after reconstruction. | Decay-e |  $\nu_{e}$ CC | NCQE | $\mu/\pi$
---|---|---|---|---
Predicted | 218.5 | 121.4 | 81.9$-$138.4 | 10.6$-$67.2
Observed | 169.5 | 109.9 | 88.2 | 27.9

| $20^{\circ}<\theta_{\rm C}<38^{\circ}$ | $38^{\circ}<\theta_{\rm C}<50^{\circ}$ | $78^{\circ}<\theta_{\rm C}<90^{\circ}$
---|---|---|---
0 or $>$1 ntags |  |  |
1 ntag |  |  |

| $20^{\circ}<\theta_{\rm C}<38^{\circ}$ | $38^{\circ}<\theta_{\rm C}<50^{\circ}$ | $78^{\circ}<\theta_{\rm C}<90^{\circ}$
---|---|---|---
0 or $>$1 ntags |  |  |
1 ntag |  |  |

| $20^{\circ}<\theta_{\rm C}<38^{\circ}$ | $38^{\circ}<\theta_{\rm C}<50^{\circ}$ | $78^{\circ}<\theta_{\rm C}<90^{\circ}$
---|---|---|---
0 or $>$1 ntags |  |  |
1 ntag |  |  |

| $20^{\circ}<\theta_{\rm C}<38^{\circ}$ | $38^{\circ}<\theta_{\rm C}<50^{\circ}$ | $78^{\circ}<\theta_{\rm C}<90^{\circ}$
---|---|---|---
0 or $>$1 ntags |  |  |
1 ntag |  |  |

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

###  D.2 Modeling of spallation backgrounds

To model the spallation background PDF we vary the isotope composition of the
background as described in Equation 3 and take the nominal spectrum as the
average between the two extreme slopes. We parameterize this average by the
following function:

$\displaystyle S_{\rm spall}(E_{\rm rec})$
$\displaystyle=N\exp\left(-\frac{(E_{\rm
rec}+0.511~{}\text{MeV})^{\alpha}}{\beta}\right),$ (27)

where $N$ is a normalization factor and $\alpha,\beta$ are determined by a fit
and depend on the SK phase considered. Then, as described in Section VII.3.2,
we both vary the isotope composition of the spallation spectrum and distort it
to account for energy scale and resolution uncertainties. We then fit the
resulting extreme spectra by the function from Equation 27 times a third order
polynomial:

$\displaystyle S^{\prime}_{\rm spall}(E_{\rm rec})=$ $\displaystyle
N^{\prime}\exp\left(-\frac{(E_{\rm
rec}+0.511~{}\text{MeV})^{\alpha}}{\beta}\right)$
$\displaystyle\times\left[1+\mathcal{P}_{3}(E_{\rm rec})\right],$ (28)

and thus parameterize the 1$\sigma$ systematic error distorsions. The nominal
spectrum and its associated analytical fit are shown in Fig. 37 for SK-IV with
its 1$\sigma$ uncertainty range. Note that, in order to account for the
steepness of the spallation spectrum while avoiding unphysical results, we
impose $S^{\prime}_{\rm spall}(E_{\rm rec})=0$ for energies larger than the
minimal $E_{\rm rec}$ for which the distorted PDF is zero.

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

## References

  * Fukugita  _et al._ (1994) M. Fukugita _et al._ , Springer ISBN 978-4-431-67029-2 (1994).
  * Kotake  _et al._ (2006) K. Kotake _et al._ , Reports on Progress in Physics 69, 971 (2006).
  * Nakazato  _et al._ (2013) K. Nakazato _et al._ , The Astrophysical Journal Supplement Series 205, 2 (2013).
  * Tammann  _et al._ (1994) G. Tammann _et al._ , The Astrophysical Journal Supplement 92, 487 (1994).
  * Diehl  _et al._ (2006) R. Diehl _et al._ , Nature 439, 45 (2006).
  * Nakazato  _et al._ (2015) K. Nakazato _et al._ , The Astrophysical Journal 804, 75 (2015).
  * Malaney (1997) R. A. Malaney, Astroparticle Physics 7, 125 (1997).
  * Totani and Sato (1995) T. Totani and K. Sato, Astroparticle Physics 3, 367 (1995).
  * Kaplinghat  _et al._ (2000) M. Kaplinghat _et al._ , Physical Review D 62, 043001 (2000).
  * Lunardini (2009) C. Lunardini, Physical Review Letters 102, 231101 (2009).
  * Horiuchi  _et al._ (2018) S. Horiuchi _et al._ , Monthly Notices of the Royal Astronomical Society 475 (2018).
  * Kresse  _et al._ (2021) D. Kresse _et al._ , The Astrophysical Journal 909, 2 (2021).
  * Horiuchi  _et al._ (2021) S. Horiuchi _et al._ , Physical Review D 103, 043003 (2021).
  * Tabrizi and Horiuchi (2021) Z. Tabrizi and S. Horiuchi, Journal of Cosmology and Astroparticle Physics 2021, 05 (2021).
  * Galais  _et al._ (2010) S. Galais _et al._ , Physical Review D 81, 053002 (2010).
  * Horiuchi  _et al._ (2009) S. Horiuchi _et al._ , Physical Review D 79, 083013 (2009).
  * Ando  _et al._ (2003) S. Ando _et al._ , Astroparticle Physics 18, 307 (2003).
  * Hartmann and Woosley (1997) D. H. Hartmann and S. E. Woosley, Astroparticle Physics 7, 137 (1997).
  * bib (2005) _Next Generation of Nucleon Decay and Neutrino Detectors: Proceedings, Aussois, Savoie, France, 7-9 April 2005_ (2005).
  * Cokinos and Melkonian (1977) D. Cokinos and E. Melkonian, Physical Review C 15, 1636 (1977).
  * Malek  _et al._ (2003) M. Malek _et al._ (Super-Kamiokande Collaboration), Physical Review Letters 90, 061101 (2003).
  * Bays  _et al._ (2012) K. Bays _et al._ (Super-Kamiokande Collaboration), Physical Review D 85, 052007 (2012).
  * Zhang  _et al._ (2015) H. Zhang _et al._ (Super-Kamiokande Collaboration), Astroparticle Physics 60, 41 (2015).
  * Gando  _et al._ (2012) A. Gando _et al._ (KamLAND Collaboration), The Astrophysical Journal 745 (2012).
  * Aharmim  _et al._ (2006) B. Aharmim _et al._ (SNO Collaboration), The Astrophysical Journal 653 (2006).
  * Aharmim  _et al._ (2020) B. Aharmim _et al._ (SNO Collaboration), Physical Review D 102, 062006 (2020).
  * Agostini  _et al._ (2021) M. Agostini _et al._ (Borexino Collaboration), Astroparticle Physics 125, 102509 (2021).
  * Fukuda  _et al._ (2003) S. Fukuda _et al._ (Super-Kamiokande Collaboration), Nuclear Instruments and Methods in Physics Research, Section A 501, 418 (2003).
  * Kume  _et al._ (1983) H. Kume _et al._ , Nuclear Instruments and Methods in Physics Research, Section A 205 (1983).
  * Abe  _et al._ (2014a) K. Abe _et al._ (Super-Kamiokande Collaboration), Nuclear Instruments and Methods in Physics Research, Section A 737, 253 (2014a).
  * Nakahata  _et al._ (1999) M. Nakahata _et al._ (Super-Kamiokande Collaboration), Nuclear Instruments and Methods in Physics Research, Section A 421, 113 (1999).
  * Blaufuss  _et al._ (2001) E. Blaufuss _et al._ (Super-Kamiokande Collaboration), Nuclear Instruments and Methods in Physics Research, Section A 458, 638 (2001).
  * Nishino  _et al._ (2009) H. Nishino _et al._ , Nuclear Instruments and Methods in Physics Research, Section A 610, 710 (2009).
  * Watanabe  _et al._ (2009) H. Watanabe _et al._ (Super-Kamiokande Collaboration), Astroparticle Physics 31, 320 (2009).
  * Hosaka  _et al._ (2006) J. Hosaka _et al._ (Super-Kamiokande Collaboration), Physical Review D 73, 112001 (2006).
  * Cravens  _et al._ (2008) J. Cravens _et al._ (Super-Kamiokande Collaboration), Physical Review D 78, 032002 (2008).
  * Abe  _et al._ (2011) K. Abe _et al._ (Super-Kamiokande Collaboration), Physical Review D 83, 052010 (2011).
  * Abe  _et al._ (2016) K. Abe _et al._ (Super-Kamiokande Collaboration), Physical Review D 94, 052010 (2016).
  * Smy (2008) M. Smy, Proceedings of the 30th International Cosmic Ray Conference 5, 1279 (2008).
  * Salpeter (1955) E. E. Salpeter, The Astrophysical Journal 121, 161 (1955).
  * Keil  _et al._ (2003) M. T. Keil _et al._ , The Astrophysical Journal 590, 971 (2003).
  * Strumia and Vissani (2003) A. Strumia and F. Vissani, Physics Letters B 564 (2003).
  * Brun  _et al._ (1994) R. Brun _et al._ , Report No. CERN-W5013 (1994).
  * Ankowski  _et al._ (2012) A. Ankowski _et al._ , Physical Review Letters 108, 052505 (2012).
  * Abe  _et al._ (2014b) K. Abe _et al._ (T2K Collaboration), Physical Review D 90, 072012 (2014b).
  * Wan  _et al._ (2019) L. Wan _et al._ (Super-Kamiokande Collaboration), Physical Review D 99, 032005 (2019).
  * Abe  _et al._ (2019) K. Abe _et al._ (T2K Collaboration), Physical Review D 100, 112009 (2019).
  * Honda  _et al._ (2007) M. Honda _et al._ , Physical Review D 75, 043006 (2007).
  * Honda  _et al._ (2011) M. Honda _et al._ , Physical Review D 83, 123001 (2011).
  * bib (Flux) http://www.icrr.u-tokyo.ac.jp/~mhonda/nflx2011/index.html (HKKM2011 Flux).
  * Hayato (2009) Y. Hayato, Acta Physica Polonica B 40, 2477 (2009).
  * Abe  _et al._ (2020) K. Abe _et al._ (Super-Kamiokande Collaboration), arXiv:2012.03807 [hep-ex] (2020).
  * Jiang  _et al._ (2019) M. Jiang _et al._ (Super-Kamiokande Collaboration), Progress of Theoretical and Experimental Physics 2019, 053F01 (2019).
  * Battistoni  _et al._ (2007) G. Battistoni _et al._ , AIP Conference Proceedings 896, 31 (2007).
  * Li and Beacom (2014) S. W. Li and J. Beacom, Physical Review C 89, 045801 (2014).
  * Li and Beacom (2015a) S. W. Li and J. Beacom, Physical Review D 91, 105005 (2015a).
  * Li and Beacom (2015b) S. W. Li and J. Beacom, Physical Review D 92, 105033 (2015b).
  * Zhang  _et al._ (2016) Y. Zhang _et al._ (Super-Kamiokande Collaboration), Physical Review D 93, 012004 (2016).
  * bib (eact) https://github.com/Goldie643/SKReact (SKReact).
  * Baldoncini  _et al._ (2015) M. Baldoncini _et al._ , Physical Review D 91, 065002 (2015).
  * int (2005) Technical Reports Series No. 428 (International Atomic Energy Agency, Vienna, 2005).
  * Locke  _et al._ (2021) S. M. Locke _et al._ (Super-Kamiokande Collaboration), to appear (2021).
  * Carminati (2011) G. Carminati, in  _32nd International Cosmic Ray Conference_ (2011).
  * Carminati (2015) G. Carminati (Super-Kamiokande collaboration), Physics Procedia 61, 666 (2015).
  * Elnimr (2017) M. Elnimr (Super-Kamiokande collaboration), Journal of Physics Conference Series. 888, 012189 (2017).
  * Conner (1997) Z. Conner, Ph.D Thesis, Boston University (1997).
  * Desai (2004) S. Desai, Ph.D Thesis, Boston University (2004).
  * Abe  _et al._ (2017) K. Abe _et al._ (T2K Collaboration), Physical Review D 96, 092006 (2017).
  * Zhang (2012) H. Zhang, Ph.D Thesis, Tsinghua University (2012).
  * Zhang (2015) Y. Zhang, Ph.D Thesis, Tsinghua University (2015).
  * Nakano  _et al._ (2020) Y. Nakano _et al._ , Nuclear Instruments and Methods in Physics Research, Section A 977, 164297 (2020).
  * Chen and Guestrin (2016) T. Chen and C. Guestrin, in  _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_, KDD ’16 (ACM, New York, NY, USA, 2016) pp. 785–794.
  * Akutsu (2019) R. Akutsu, Ph.D Thesis, The University of Tokyo (2019).
  * Abe  _et al._ (2021a) K. Abe _et al._ (Super-Kamiokande Collaboration), to appear (2021a).
  * Abe  _et al._ (2021b) K. Abe _et al._ (Super-Kamiokande Collaboration), to appear (2021b).
  * Rolke  _et al._ (2005) W. A. Rolke _et al._ , Nuclear Instruments and Methods in Physics Research, Section A 551, 493 (2005).
  * Benhar  _et al._ (1994) O. Benhar _et al._ , Nuclear Physics A 579, 493 (1994).
  * Benhar  _et al._ (2005) O. Benhar _et al._ , Physical Review D 72, 053005 (2005).
  * Horowitz and Serot (1981) C. J. Horowitz and B. D. Serot, Nuclear Physics A 368, 503 (1981).
  * Maieron  _et al._ (2003) C. Maieron _et al._ , Physical Review C 68, 048501 (2003).
  * Caballero  _et al._ (2005) J. A. Caballero _et al._ , Physical Review Letters 95, 252502 (2005).
  * Gonzalez-Jimenez  _et al._ (2013) R. Gonzalez-Jimenez _et al._ , Physical Review C 88, 025502 (2013).
  * Amaro  _et al._ (2005) J. E. Amaro _et al._ , Physical Review C 71, 015501 (2005).
  * Amaro  _et al._ (2006) J. E. Amaro _et al._ , Physical Review C 73, 035503 (2006).
  * Capuzzi  _et al._ (1991) F. Capuzzi _et al._ , Nuclear Physics A 524, 681 (1991).
  * Meucci  _et al._ (2004) A. Meucci _et al._ , Nuclear Physics A 739, 277 (2004).
  * Meucci and Giusti (2014) A. Meucci and C. Giusti, Physical Review D 89, 057302 (2014).
  * Note (1) The latest DSNB search result from KamLAND Abe  _et al._ (2021c), while still unpublished, serves tighter or comparable limits below $E_{\nu}=13.3$ MeV.
  * Ashie  _et al._ (2005) Y. Ashie _et al._ (Super-Kamiokande Collaboration), Physical Review D 71, 112005 (2005).
  * Regis  _et al._ (2012) C. Regis _et al._ (Super-Kamiokande), Physical Review D 86, 012006 (2012).
  * Winter  _et al._ (2006) W. T. Winter _et al._ , Physical Review C 73, 025503 (2006).
  * Jegerlehner  _et al._ (1996) B. Jegerlehner _et al._ , Physical Review D 54, 1194 (1996).
  * Sisti (2020) M. Sisti (JUNO Collaboration), Journal of Physics Conference Series 1468, 012150 (2020).
  * Abe  _et al._ (2021c) S. Abe _et al._ (KamLAND Collaboration), arXiv:2108.08527 [astro-ph.HE] (2021c).