Diffuse Supernova Neutrino Background Search at Super-Kamiokande

Core-collapse supernovae (CCSNe) are among the most cataclysmic phenomena in the Universe and are essential elements of dynamics of the cosmos. Their underlying mechanism is however still poorly understood, as characterizing it would require intricate knowledge of the core of the collapsing star. Information about this core could be accessed by detecting neutrinos emitted by supernova bursts, whose luminosity and energy spectra closely track the different steps of the CCSN mechanism in a neutrino-heating scenario (see Refs. Fukugita et al. (1994); Kotake et al. (2006); Nakazato et al. (2013), for example). Existing neutrino experiments, however, are mostly sensitive to supernova bursts occurring in our galaxy and its immediate surroundings, that are extremely rare (a few times per century in our galaxy Tammann et al. (1994); Diehl et al. (2006)). Learning about the aggregate properties of supernovae in the Universe hence requires observing the accumulation of neutrinos from all distant supernovae. The integrated flux of these neutrinos forms the diffuse supernova neutrino background, or DSNB (DSNB neutrinos are also referred to as supernova relic neutrinos (SRNs) in various publications).

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

While a significant fraction of DSNB neutrinos are expected to have energies lower than 10 MeV —as shown in Fig. 1—, the $\mathcal{O}(1)$ MeV region is hard to probe by current experiments due to overwhelming backgrounds from reactor antineutrinos, solar neutrinos, and radioactivity. Experimentally the DSNB signal has therefore been searched for at energies of $\mathcal{O}(10)$ MeV. At these energies, the dominant detection channel in most experiments is the inverse beta decay (IBD) of electron antineutrinos ($\bar{\nu}_{e}+p\rightarrow e^{+}+n$), and contributions from subleading channels such as electron neutrino elastic scattering or charged current interactions on oxygen can be neglected. This process produces an easily identifiable positron and a neutron (Fig. 2). Here the neutron does not emit Cherenkov or scintillation light directly but its capture on a proton leads to emission of a $2.2$ MeV photon. In pure water, the characteristic timescale for neutron moderation and capture is of $204.8\pm 0.4~{}\mu$s Cokinos and Melkonian (1977). Pairing the “prompt” positron signal with the “delayed” photon emission from neutron capture is hence a key component of many analyses.

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

In this study, we draw on the previous SK analyses to present two DSNB searches for antineutrino energies ranging from $9.3$ to $81.3$ MeV, with significantly improved background modeling and reduction techniques. In the $9.3$ to $31.3$ MeV range, we derive differential upper limits on the $\bar{\nu}_{e}$ flux independently from the DSNB model, following the strategy outlined in Ref. Zhang et al. (2015) using a $22.5\times 2970$-kton$\cdot$day data set. In the $17.3$ to $81.3$ MeV range we constrain a wide variety of DSNB models using spectral fits analogous to the ones described in Ref. Bays et al. (2012). We then combine the results of this analysis with the ones obtained in Ref. Bays et al. (2012) for the former SK phases, thus analyzing $22.5\times 5823$-kton$\cdot$days of data. This unprecedented exposure will allow to probe the DSNB with an unmatched sensitivity.

The rest of this article proceeds as follows. First, the SK detector and the specific features of its data acquisition system are described in Section II. Then, details about modeling of the DSNB signal and the different backgrounds are given in Sections III and IV, respectively. We then describe the data reduction process in Section V. The procedures associated with the DSNB model-independent and spectral analyses are described in Sections VI and VII, respectively. Finally we discuss the current constraints on the DSNB flux and future opportunities at SK with gadolinium and Hyper-Kamiokande in Section VIII before concluding in Section IX.

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

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

SK-IV has been the longest data taking phase of the experiment, lasting about $10$ years. During this term, the PMT gain has steadily increased by about 15% over the whole SK-IV period. Additionally, since the energy scale change due to the gain change for this analysis was about 3 %, that effect was corrected in order to maintain the stability of the detector. A key ingredient of our DSNB study is therefore a simulation of antineutrino IBD processes that accounts for the evolution of the detector’s properties throughout the entire SK-IV period.

In this study we simulate IBD processes uniformly distributed in the entire inner detector, in order to account for events close to the ID wall being misreconstructed inside the fiducial volume. For each vertex we generate a set of three momentum vectors corresponding to an antineutrino, a positron, and a neutron. In a view of the wide diversity of DSNB models, we generate uniformly distributed positron energies in 1$-$90 MeV and renormalize events later to model specific spectra. This procedure will also allow to use this simulation to model backgrounds associated to other IBD processes or $\beta$ decays, as will be mentioned later. We then simulate detection signals using a dedicated simulation, based on GEANT3 Brun et al. (1994), that reproduces the properties of the water in SK, as well as models the PMT properties and electronic response. Detector properties such as water transparency and PMT noise are being measured daily, allowing to accurately model the time-dependent noise contamination of the positron signal. Neutron capture, however, produces a $2.2$ MeV $\gamma$ ray whose light yield is below the SK trigger threshold and particularly difficult to distinguish from PMT dark noise and low energy radioactive decays. In order to develop robust neutron identification techniques we therefore collect noise samples from data using a random wide trigger, and inject them into simulation results, after the positron activity peak.

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

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

In spite of its 2700-m water-equivalent overburden, SK is still exposed to cosmic ray muons with a rate of $\sim$2 Hz. Interactions of these muons in water produce electromagnetic and hadronic showers. Spallation of oxygen nuclei induced by the muons or by secondary particles can then lead to the production of radioactive isotopes, whose decays can be misidentified as IBD events in the DSNB search window. Below about 20 MeV, the associated background is $10^{6}$ times higher than the DSNB flux predictions, making spallation reduction an essential aspect of this analysis. Since most spallation isotope decays only produce $\beta$ and $\gamma$ rays, spallation backgrounds can be significantly reduced using neutron tagging. Nonetheless, accidental pairing between prompt $\beta$ or $\gamma$ events and PMT hits due to dark noise or intrinsic radioactivity will allow a sizable number of spallation events to pass as IBDs. Furthermore, a few isotopes undergo a $\beta$+$n$ decay that mimics the IBD signal, and thus cannot be removed using neutron tagging. Efficient spallation reduction therefore requires both dedicated spallation reduction techniques tailored to the properties of the different isotopes, and neutron tagging cuts.

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

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

Antineutrinos emitted from nearby nuclear reactors can undergo IBDs in SK and therefore be mistaken as the DSNB signal at low energies. We estimate the reactor antineutrino flux using the SKReact code bib (eact) developed for the SK reactor neutrino analysis. This code uses the reactor neutrino model from Ref. Baldoncini et al. (2015) based on IAEA records int (2005) to evaluate the electron antineutrino fluxes for any reactor in the world, and simulates oscillation effects. The resulting flux prediction also accounts for the time dependence of the reactor antineutrino flux due to changes in reactor activity. The expected total reactor neutrino flux at SK can thus be predicted up to $E_{\nu}^{\text{max}}=9.2$ MeV. The expected reactor neutrino spectrum is obtained from this predicted flux using the IBD simulation, and is shown in Fig. 5. Above 6 MeV, as can be seen in this figure, the reconstructed spectrum is primarily shaped by resolution effects. Since the reactor antineutrino spectrum has been constrained up to 12 MeV by Daya Bay An et al. (2017), possible antineutrino contributions above $E_{\nu}^{\text{max}}$ can be safely neglected in this study. In the 7.5$-$9.5 MeV reconstructed kinetic electron energy range, the reactor antineutrino background is about 6 times higher than the DSNB signal predicted by the Horiuchi$+$09 6 MeV model Horiuchi et al. (2009). Reactor antineutrino backgrounds thus effectively set a lower energy threshold for the DSNB search.

Electron neutrinos from the Sun cause elastic interactions with electrons in SK. In the DSNB search window, the ⁸B and hep solar neutrino fluxes represent an important background up to $\sim$20 MeV. This background can be drastically reduced by neutron tagging, although a small fraction of electrons from solar neutrinos will be accidentally paired with background fluctuations or radioactive decays, like for the spallation backgrounds. In the absence of neutron tagging, a solar neutrino interaction can still be identified by comparing the reconstructed directions of the scattered electron to the direction of the Sun at the detection time. The angle between these two directions, called the opening angle, is expected to be close to zero for solar electron neutrinos scattering elastically. It is however smeared by reconstruction effects that need to be modeled accurately. In this study, we use the simulation designed for the latest SK-IV solar neutrino analysis Abe et al. (2016) to model solar neutrino spectra and evaluate the impact of opening angle cuts.

Finally, note that in this analysis we do not study the impact of solar antineutrino production, since the rate predicted by the Standard Model (SM) is negligible. Antineutrino production via beyond the SM processes would yield a signal that can only be distinguished from the DSNB signature by spectral shape studies. So far, only upper limits on the solar antineutrino flux have been derived at SK, notably in Ref. Abe et al. (2020).

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

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

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

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

The introduction of a new trigger scheme at SK-IV has made characterizing IBDs via neutron tagging possible, by requiring the prompt and delayed signals be detected within 500 $\mu$s of each other. In what follows, we will devise a neutron tagging algorithm tailored to the DSNB analysis, inspired by previous SK studies from Refs. Zhang (2012, 2015).

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

In this analysis we use the method described in Ref. Bays et al. (2012) and apply $\cos\theta_{\rm sun}$ cuts in the four MSG bins shown in Fig. 17. We estimate the impact of these cuts on DSNB signal events by assuming that these events have a uniformly distributed $\theta_{\rm sun}$, and we obtain the associated MSG distributions using the IBD Monte-Carlo simulation. Conversely, to estimate the number of solar neutrino events above 15.5 MeV, we evaluate the number of solar events in 13.5-15.5 MeV from data and extrapolate the result to higher energies using the solar neutrino MC simulation from Ref. Abe et al. (2021c), a procedure already used in the DSNB search at SK-I,II,III Bays et al. (2012). Here, we can safely assume the $\cos\theta_{\text{sun}}$ distribution for non-solar event to be uniform Abe et al. (2016), and therefore approximate the number of solar events in 13.5-15.5 MeV by:

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

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

Here we discuss the behavior of the observed data after the cuts. In particular, the impact of spallation, $d_{\mathrm{eff}}$, and the other positron candidate selection cuts on data is shown in Fig. 19. As shown in this figure, for reconstructed energies below 20 MeV contributions from spallation and radioactivity near the wall are significant. To validate our background modeling and neutron tagging procedure we also compare the observed data to the background expectations after the noise reduction, multiple spallation cuts, DSNB positron candidate selection, and neutron tagging cuts. To simplify the interpretation of the data, we impose the same neutron tagging cut, with a 20% signal efficiency, for all energies. While most spallation cuts are removed in order to focus on neutron tagging, multiple spallation cuts are conserved in order to avoid contaminating the neutron detection window with isotope decays —an effect not accounted in the mistag rate predictions. The predicted and observed event spectra are shown in Fig. 20 and agree within uncertainties.

After the cuts derived in Section V.6 are applied, we find 102 events with exactly one neutron in 7.5$-$79.5 MeV. For these events, the distribution of the time difference between the neutron and the prompt event is shown in Fig. 21. Fitting this distribution by an exponential plus a constant, as is done in the ²⁴¹Am/Be calibration, yields a time constant of $\tau=154\pm 104~{}\mu$s, which is compatible with the expected neutron capture time in water. In the 15.5$-$79.5 MeV region, that will be used for the spectral analysis, we also find 248 events with zero neutrons and 9 events with more than one neutron. The observation times and reconstructed vertices for these events are shown in Figs. 22 and 23. No event time cluster has been observed and the time-dependence of the event rate is consistent with the livetime variations over the SK-IV period. Events are also uniformly distributed all over the FV, irrespective of the number of neutrons. The space and time distributions of the selected events are therefore consistent with what is expected for DSNB candidates.

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

In order to evaluate the number of atmospheric neutrino background events, we define three regions of parameter space based on the reconstructed Cherenkov angle: one signal region with $\theta_{\rm C}\in[38,50]^{\circ}$, that will contain most of the signal and the irreducible backgrounds, and two sidebands with $\theta_{\rm C}\in[20,38]^{\circ}$ and $\theta_{\rm C}\in[78,90]^{\circ}$. The low Cherenkov angle region will be populated with mostly atmospheric backgrounds involving visible muons and pions while the high angle region will be mostly populated by NCQE atmospheric neutrino events with multiple $\gamma$ rays. Finally, we separate events with exactly one identified neutron from the others, thus defining an “IBD-like” and a “non IBD-like” region. Note that due to the low efficiencies of the neutron tagging cuts the non IBD-like region is expected to contain a sizable amount of signal. Our analysis will hence involve six regions of parameter space: two signal regions with intermediate values of the Cherenkov angle, and four sidebands with low and high Cherenkov angle values, as summarized in Table 6.

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

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

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.

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

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

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.