Exploring the Fate of Stellar Core Collapse with Supernova Relic Neutrinos

Core collapse of massive stars leads to different fates for various physical factors, which gives different spectra of the emitted neutrinos. We focus on the supernova relic neutrinos (SRNs) as a probe to investigate the stellar collapse fate. We present the SRN fluxes and event rate spectra at a detector for three resultant states after stellar core collapse, the typical mass neutron star, the higher mass neutron star, or the failed supernova forming a black hole, based on different nuclear equations of state. Then possible SRN fluxes are formed as mixtures of the three components. We also show the expected sensitivities at the next-generation water-based Cherenkov detectors, SK-Gd and Hyper-Kamiokande, as constraining the mixture fractions. This study provides a practical example of extracting astrophysical constraints through SRN measurement.


INTRODUCTION
Massive stars are likely to cause huge explosions with a gravitational collapse of their cores, as referred to as core-collapse supernovae (SNe). They are considered to be the major birthplace of neutrinos and therefore detecting neutrinos would promote understanding of the mechanism of SNe. Specifically, neutrinos could access physics at the most inner core, which is not possible via optical observations. So far, the detection of neutrinos from supernova bursts was only achieved for SN1987A (Alekseev et al. 1987;Bionta et al. 1987;Hirata et al. 1987), because of a low supernova rate as well as detection difficulties of neutrinos. As a complementary strategy to the direct detection of neutrino bursts, measuring the accumulated flux from all distant SNe does not require us to wait for the rare explosion within a detectable range. This integrated flux is referred to as supernova relic neutrinos (SRNs), or the diffuse supernova neutrino background (DSNB) in some articles. The flux shape of SRNs depends on multiple factors, including not only the original neutrino spectrum from an explosion but also cosmic evolution, etc. Therefore, SRNs are believed to provide valuable information for further understanding of these physics. However, the observation has yet to be achieved, while a recent search at the Super-Kamiokande (SK) water Cherenkov detector has entered the range of theoretical predictions (Abe et al. 2021). Figure 1 shows the latest experimental upper limits from SK (Abe et al. 2021) and KamLAND (Abe et al. 2022), in comparison to predictions from our models in this paper. Various detectors such as water Cherenkov, liquid scintillators, xenon-based, and others in the next few decades are expected to achieve the SRN observation by improving sensitivity dramatically (Priya & Lunardini 2017;Møller et al. 2018;de Gouvêa et al. 2020;Sawatzki et al. 2021;Li et al. 2022;Suliga et al. 2022). Accordingly, there is a higher need for the preparation to extract physical constraints when SRNs are actually detected.
In this study, we focus on the fate of stellar core collapse and its impact on the emitted neutrino spectrum. Total energy of neutrinos emitted from a core-collapse SN corresponds to the binding energy of a neutron star (NS) formed after the explosion (Sato & Suzuki 1987;Lattimer & Yahil 1989), which is determined by its mass and the nuclear equation of state (EOS). While the masses of NSs observed so far are typically ∼ 1.35M , high-mass NSs with 1.6M were also discovered (Alsing et al. 2018). If a black hole (BH) is formed by core collapse, the spectrum of emitted neutrinos becomes different from the case of the SN with an NS and the total emission energy depends strongly on the EOS (Sumiyoshi et al. 2006). Therefore, the mass distribution of remnant NSs and the fraction of BH-forming collapses are responsible for the event rate of SRNs, and  (Abe et al. 2021) and KamLAND (Abe et al. 2022). The SRNνe fluxes estimated under the assumption that neutrinos originate from a single component of the SN with a canonical-mass neutron star (fHNS = 0 and fBH = 0), the SN with a high-mass neutron star (fHNS = 1 and fBH = 0), and the failed SN (fBH = 1) are shown with red, blue, and green lines, respectively. The solid and dashed lines correspond to the results for normal and inverted neutrino mass hierarchies (NH and IH), respectively. Here the fluxes with three different nuclear equations of state in Table 1 are shown in the same style. Detailed descriptions about modeling are given in Section 2. which factor is more largely responsible would depend on the nuclear EOS. Hence, measuring the SRN spectrum could provide implications about the core-collapse fate. 1 This paper is structured as follows. We introduce recent theoretical studies and describe our models in Section 2 and give experimental knowledge and assumptions for the sensitivity study in Section 3. Then we present the results and discussion in Section 4. Finally, we conclude in Section 5. The flux of SRNs is related to the evolution of galaxies that determines the cosmic rate of stellar core-collapse 1 Note that uncertainties about many other factors affecting the SRN flux should be carefully considered to access the NS mass and BH formation fraction.
events (Totani et al. 1996). In order to obtain the cosmic rate of stellar core-collapse events, we utilized the cosmological simulation results of the Illustris-1 run, which is the flagship of the Illustris simulation project (Vogelsberger et al. 2014;Nelson et al. 2015). The cosmic core-collapse rate at a redshift of z is calculated as where the cosmic star formation rateρ * (z) is adapted from the data of the Illustris-1 run which are publicly released for each snapshot at the selected redshift. Following the Illustris model (Vogelsberger et al. 2013), we choose the initial mass function (IMF) of progenitors ψ IMF (M ) proposed by Chabrier (2003). Consequently, the maximum and minimum masses of progenitors that cause core collapse are set to M max = 100M and M min = 6M , respectively. Note that the minimum mass for stars ending as a type II SN is assumed to be 6M in Vogelsberger et al. (2013) although larger values (8-10M ) are adopted for M min in the previous SRN studies (Nakazato et al. 2015;Horiuchi et al. 2018;Kresse et al. 2021). This leads to the cosmic core-collapse rate at z = 0 of R CC (0) = 3.9×10 −4 Mpc −3 yr −1 , which is significantly higher than the typical value of the core-collapse SN rate obtained from observations (e.g., Madau & Dickinson 2014). Nevertheless, the cosmic core-collapse rate of around z = 1 is comparable to that of Mathews et al. (2014), which is often adopted in previous SRN studies. A recent study using the Gaia DR2 database suggests that the upper limit on white dwarf formation is ∼ 6M (Richer et al. 2021). 2 Based on these settings, we obtain R CC (z) ρ * (z)/(57.5M ),ρ * (z)/(84.6M ), andρ * (z)/(115M ) for M min = 6M , 8M , and 10M , respectively. On the contrary, according to Horiuchi et al. (2021), mergers of binary stars whose individual masses are originally below the core-collapse threshold can produce SN progenitors and enhance the SRN flux. In addition, the redshift evolution of the IMF could be taken into account for the estimation of the SRN event rate (Aoyama et al. 2021). We classify the fate of stellar collapse into four types: the SN with a low-mass NS, the SN with a high-mass NS, the faint SN, and the failed SN. The low-and high-mass NSs are divided by a baryon mass of M NS,b = 1.6M . 3 Faint SNe are successful explosions accompanied by a substantial fallback that leads to a BH formation at later times. Since they are considered to be rare, we will not take into account their contribution in the following study of this paper. In failed SNe, a BH is formed by continuous mass accretion without explosion. The fate of stellar core collapse strongly affects the neutrino emission. In the previous SRN studies (Nakazato et al. 2015;Horiuchi et al. 2018;Kresse et al. 2021), the fate is assumed to be determined by the progenitor model, and the IMF-weighted neutrino spectrum is used to calculate the SRN flux. Nevertheless, this approach may be inconvenient because the evolutionary scenario for progenitors and explosion mechanism of core-collapse SNe are still uncertain. In other words, the stellar core-collapse fate depends not only on the progenitor mass and metallicity but also on the evolutionary model.
In the present study, we consider three cases for the fate of stellar core collapse: the SN with a canonicalmass NS, the SN with a high-mass NS, and the failed SN forming a BH. 4 Consequently, we set their fractions directly as parameters and define the fate-weighted average neutrino number spectrum as where f BH is the fraction of failed SNe to the total number of stellar core-collapse events and f HNS is the fraction of SNe with a high-mass NS to the total number of SNe with an NS. Furthermore, we employ generic models of the neutrino spectrum dN CNS (E ν )/dE ν , dN HNS (E ν )/dE ν , and dN BH (E ν )/dE ν for the SN with a canonical-mass NS, the SN with a high-mass NS, and the failed SNe with a BH, respectively. We then describe the SRN flux as 3 According to the approximation by Lattimer & Yahil (1989), neutron stars with M NS,b = 1.6M are expected to have a gravitational mass of M NS,g ≈ 1.43M . Note that the relation between the baryon mass and gravitational mass of neutron stars depends on the EOS while the accuracy of this approximation would be at most ±2%. 4 After the first submission of this paper, Ekanger et al. (2022) showed the prediction of SRN event rate focusing on the NS mass.
where c is the speed of light and the cosmological constants are Ω m = 0.2726, Ω Λ = 0.7274, and H 0 = 70.4 km s −1 Mpc −1 (Vogelsberger et al. 2014). Here, the neutrino energy at a detector E ν is related to that at a source E ν as E ν = (1 + z)E ν when the source is located at a redshift z.
The generic models of the neutrino spectrum for each component are adopted from the results of numerical simulations. As for the SN with an NS, we refer to Nakazato et al. (2022) and employ the models with NS baryonic masses of M NS,b = 1.47M and 1.86M for the representatives of a canonical-mass NS and a high-mass NS, respectively. Furthermore, in order to obtain the time-integrated neutrino spectra, we prepare the neutrino light curves from an onset of core collapse to cooling of the proto-NS by following the method adopted in the Supernova Neutrino Database . The shock revival time used in this procedure is taken from Table 1 of Nakazato et al. (2022) while it is denoted as t init . The progenitors of our SN models are taken from Woosley & Weaver (1995); the models with M NS,b = 1.47M and 1.86M have progenitors of the solar metallicity with zero-age main-sequence masses of 15M and 40M , respectively. On the other hand, for the failed SN, we utilize the models in Nakazato et al. (2021), where the core collapse and neutrino emission are studied for the same progenitor model (an initial mass of 30M and metallicity of Z = 0.004) involving a BH formation as in Nakazato et al. (2013).
The neutrino signal from stellar core collapse depends on the nuclear EOS. In the case of the SN with an NS, the total energy of emitted neutrinos is determined by the gravitational binding energy of the remnant NS and the EOS affects the average energy of neutrinos (Nakazato et al. 2018) and the cooling timescale of the proto-NS (Nakazato & Suzuki 2019. In particular, although the Supernova Neutrino Database ) provides the neutrino spectra up to 20 sec, the neutrino emission continues for a longer period (Suwa et al. 2019) and the EOS dependence gets stronger for the later phase (Nakazato et al. 2022). In the case of the failed SN, the time-integrated neutrino spectrum depends on the maximum mass of the proto-NS because neutrinos are emitted until a BH formation (Sumiyoshi et al. 2006). In the previous studies (Lunardini 2009;Nakazato et al. 2015;Mathews et al. 2020), the EOS dependence is taken into account only for the component of the failed SN. In the present study, we also consider the EOS dependence for the SN with an NS. For this purpose, we employ three EOS models: Togashi (Togashi et al. 2017), LS220 (Lattimer & Swesty 1991), and Shen (Shen et al. 2011). Then, the neutrino spectra obtained with use of different EOS models are prepared for each component.
We show the average and total energies of our timeintegrated neutrino signal in Table 1. First, we compare the models with different fates. The average neutrino energy of the SN with a high-mass NS is almost the same as that of the SN with a canonical-mass NS. In contrast, the models with a high-mass NS give a larger total energy of neutrinos than those with a canonical-mass NS because the remnant NS has a larger binding energy available for the neutrino emission. As for the failed SN, the gravitational potential energy released from core collapse is not fully converted to the neutrino radiation since a BH is formed. Thus the total neutrino energy of the failed SN is comparable to that of the SN with an NS. Nevertheless, the failed SN has a higher average energy of neutrinos than the SN with an NS because heating due to matter accretion continues until a BH formation. Next, we consider the EOS dependence. Since the total energy of emitted neutrinos corresponds to a binding energy of the remnant NS, EOS models with a smaller NS radius have a larger total energy for the SN with an NS. Among our models, the Togashi and Shen EOSs have the smallest and largest NS radii, respectively. In the case of the failed SN, the amount of emitted neutrinos is larger for EOS models with a higher maximum mass because the time taken to form a BH and accordingly duration of the neutrino emission are longer. Note that the maximum mass of proto-NSs is different from that of cold NSs and depends on the entropy as shown in Figure 13 of Schneider et al. (2020).

Distribution of compact remnants
In our study, the distribution of compact remnants formed by stellar collapse is encapsulated in two parameters, f BH and f HNS . Observationally, a bimodal mass distribution is inferred for the NS in the binary system (Valentim et al. 2011;Özel et al. 2012;Kiziltan et al. 2013;Antoniadis et al. 2016;Alsing et al. 2018). The low-mass component of the distribution contains a majority of the sample and has a narrow peak at ∼1.35M . Canonical-mass NS models in our calculation are chosen to be consistent with this component (Table 1). On the other hand, the high-mass component, which mainly consists of the NS with a white dwarf companion, has a broad peak around 1.5-1.8M . The bimodal mass distribution is interpreted as a result of difference in the birth mass and/or additional accretion history. Antoniadis et al. (2016) claimed that a certain fraction of the observed high-mass NSs should have been born with a mass larger than 1.6M because the accretion efficiency from the binary companion was estimated to be low. In contrast, a stellar evolution model by Woosley et al. (2020) produced much fewer high-mass NSs than claimed by Antoniadis et al. (2016). Incidentally, Woosley et al. (2020) shows that the birth mass function of NSs is insensitive to the mass-loss rate.
In our calculation, f HNS corresponds to the fraction of the high-mass NS at birth, which is still not well constrained.
Failed SNe associating a BH formation would be observed as a disappearance of luminous stars without SNe in optical wavelengths (Kochanek et al. 2008). Recently, Neustadt et al. (2021) compiled 11 yr of the monitoring data for luminous stars in nearby galaxies and found the failed SN fraction to be 4-39% at a 90% C.L. provided that one failed SN was detected. This range is consistent with Woosley et al. (2020), whose model predicts a fraction of 9-32%. In this prediction, the fraction of BH-forming collapses is sensitive to the mass-loss rate. Because the mass-loss rate depends on the stellar metallicity, the failed SN fraction would be a function of the redshift; f BH is considered to be larger for higher redshifts in general (Nakazato 2013;Nakazato et al. 2015;Yüksel & Kistler 2015). Nevertheless, we assume that the fraction of failed SNe does not depend on redshift for simplicity.

Neutrino oscillation and detection
In water Cherenkov detectors, the main detection channel of SRNs is inverse beta decay (IBD) of electron antineutrinos:ν e + p → e + + n.
(4) Therefore, we focus on the flux ofν e throughout the article. For the models of neutrino spectrum described above (Table 1), the neutrino flavor oscillations are not taken into account. The terrestrial flux ofν e is written as where dΦν e (E ν )/dE ν and dΦ νx (E ν )/dE ν are fluxes of ν e and ν x at a core collapse, respectively. Note that we denote ν µ ,ν µ , ν τ , andν τ as ν x collectively because the difference in their emission processes is subtle. The effect of flavor conversions is encapsulated in a survival probability ofν e ,P . Following Nakazato et al. (2015), we consider the flavor conversions due to the matterinduced effect (the so-called the MSW effect; Wolfenstein 1978;Mikheyev & Smirnov 1985) that takes place in the stellar envelope. The flavor mixing depends on the neutrino mass hierarchy (Dighe & Smirnov 2000)  Note-For the successful SN models, MNS,g and RNS are the gravitational mass and radius of NSs, respectively. For the failed SN models, tBH is the time to a BH formation measured from the core bounce. Eν i and Eν i ,tot are the average and total energies of the time-integrated neutrino signal for νi, where νx represents the average of νµ,νµ, ντ , andντ .
asP ≈ 0.7 for the normal mass hierarchy (NH) and P ≈ 0 for the inverted mass hierarchy (IH). Also the self-induced (collective) effect can exchangeν e and ν x (e.g., Duan et al. 2010). Unfortunately, the self-induced conversion is complicated and still unknown because it is hard to solve the flavor evolution involving the selfinteraction terms in general. Thereby, the oscillation probabilities depend on the time and energy. Nevertheless, according to Lunardini & Tamborra (2012), the self-induced effect on the SRN spectrum is minor. In the present study, we do not consider the collective oscillation while further investigations are important. The SRN event rate is calculated as where σ IBD (E ν ) is the IBD cross section taken from Strumia & Vissani (2003). For simplicity, we set the positron energy to be E e + = E ν − ∆c 2 with a neutronproton mass difference, ∆, while the recoil effects that are provided in Strumia & Vissani (2003) also reduce the positron energy. N p represents the number of free protons contained in the detector.

DETECTORS AND BACKGROUND SOURCES
In this study, we focus on water-based Cherenkov detectors. This type of detectors includes SK, its upgrade with a gadolinium loading, SK-Gd, and the future large volume detector, Hyper-Kamiokande (HK). We consider the SRNν e search range and background sources at SK-Gd and HK based on the latest SK analysis with 2970 days of data presented in Abe et al. (2021), as they share many in common. The detector volume of SK-Gd is set to the same as SK (N p = 1.5 × 10 33 ) while that of HK is set to 8.4 times larger (N p = 12.6 × 10 33 ) as assuming their fiducial volumes. As a detector operation time, we consider 10 yr for SK-Gd, and 3, 5, and 10 yr for HK. The backgrounds are scaled by the assumed livetime and detector volume from Abe et al. (2021).
The detection channel is IBD of electron antineutrinos as mentioned above. At SK, positrons are detected via Cherenkov light emitted from their travel in the detector. Neutrons cannot be detected directly as they are electromagnetically neutral, but SK made use of a 2.2 MeV γ ray emitted from their capture on hydrogens after thermalization. Since the signal size of electromagnetic cascades from this 2.2 MeV γ ray is small and then difficult to distinguish from a huge amount of low energy backgrounds, a machine-learning technique was employed in the analysis. The signal efficiency with this method was 20−30%, depending on energy region, to reduce low backgrounds to a manageable size. 5 At SK-Gd, because γ rays emitted from neutron captures on gadolinium have a higher total energy (∼8 MeV) enough to be detected unambiguously, the neutron tagging efficiency is expected to be dramatically improved, depending on the gadolinium loading rate. 6 In this paper, 5 We scale from the 30% neutron tagging at SK, as the efficiency around this was achieved for the chosen criteria for the energy region of current interest in Abe et al. (2021). 6 The capture rate is ∼50% with 0.01% gadolinium in water and ∼90% with 0.1% according to Li et al. (2022).
we assume a 70% neutron tagging efficiency for SK-Gd, while a higher efficiency may be achieved in the experiment. As for HK, since there are still many unclear factors, we assume the same tagging efficinecy as that for SK. 7 We assume the energy resolutions at both of SK-Gd and HK are the same as that at SK. The analysis energy window for the SRN search was set to 9.3 < E ν < 31.3 MeV in Abe et al. (2021). In this paper, we show the results from two different energy ranges: 13.3 < E ν < 31.3 MeV and 17.3 < E ν < 31.3 MeV. The background source has its characteristic energy region in the search window; reactor neutrinos, solar neutrinos, and cosmic ray muon spallation events dominate in the lower side (E ν 17.3 MeV), while the atmospheric neutrino background becomes an almost only one contributor at higher energies (E ν 17.3 MeV).
There are two types of spallation backgrounds. One is isotopes, such as 9 Li, that make almost identical signitures as IBD via β + n decay. The other is isotopes whose decay products do not include neutrons. Most of these can be easily rejected by the neutron tagging technique, but some make an accidental pair with a 2.2 MeV γ-ray-like signal and then mimic IBD events. These accidental backgrounds are expected to be reduced dramatically at SK-Gd, thanks to much better visuality of ∼8 MeV γ rays. We therefore assume a 10 times smaller amount of accidental backgrounds at SK-Gd. In contrast, cosmic ray muon spallation may increase at HK as overburden at the planned construction mountain is smaller than that at the SK site. Nevertheless, since there are still a lot of unclear factors, we assume the same amount of spallation at HK. In the later study, note that the background assumption for the results with 13.3 < E ν < 31.3 MeV would be underestimated for this reason. The results with 17.3 < E ν < 31.3 MeV are not affected as this region is almost spallation-free.
In the SRN search at SK, atmospheric neutrino backgrounds are classified into two categories as neutralcurrent quasielastic (NCQE) interactions and the others (non-NCQE). NCQE interactions involve nuclear γray emission via primary neutrino and following neutron reactions , which mimic the positron signal from IBD. Charged-current interactions of atmospheric electron antineutrinos give the exactly same final state as IBD. Charged-current muon neutrino interactions would become backgrounds when the final state muons have less energy than the Cherenkov threshold 7 The photosensors for HK have a better performance than those for SK, but the neutron tagging performance depends also on the detector size, the water purity, the coverage of sensors, etc.
but electrons from their decays are visible in the detector. Some of neutral-current interactions involve a low energy pion, which is also invisible while an electron from its decay is visible, and then look similar to signal. We assume the same fractions of these atmospheric backgrounds at SK-Gd and HK as SK.
We estimate the expected sensitivity to the SRNν e flux by using pseudo experiments following the method in Abe et al. (2021); we vary background events by their statistical and systematic uncertainties then obtain the upper bound on the number of events at a certain C.L., N lim . The expected upper bound on theν e flux is then calculated as where T is the livetime of the detector operation,σ IBD is the IBD cross section at a mean energy in the region, and sig is the signal efficiency. We use the IBD cross section at 24.3 (22.3) MeV for the analysis with an energy range of 17.3 (13.3) < E ν < 31.3 MeV. This treatment may cause a change in the result for a nonlinear shape of the differential cross section, but the impact on the final results is limited. Here we assume multiple assumptions on systematic uncertainties. In comparison to a 60% uncertainty at SK, we use a ×0.7 size of uncertainty on NCQE for SK-Gd 10 yr. On the other hand, we assume 30%, 20%, and 10% for HK 3 yr, HK 5 yr, and HK 10 yr, respectively. 8 These are based on the expected efforts in accelerator neutrino and nuclear experiments as described in Abe et al. (2019). In Abe et al. (2021), the non-NCQE background was estimated by using the sideband region, whose uncertainty is then mostly statistical. We assume a ×0.7 size of uncertainty compared to the SK one on this for SK-Gd 10 yr, while use 10% for HK 3 yr, 8% for HK 5 yr, and 5% for HK 10 yr. We take the same size of systematic uncertainties for the other backgrounds from Abe et al. (2021). The uncertainties for these other backgrounds may be improved in the future, but the impact is subdominant in the current study. Each condition for different detectors is summarized in Table 2.  Note-The SK-IV case is shown as a reference. The efficiencies and uncertainties are actually energy dependent in the SK-IV analysis, while we assume representative values and extrapolate to each case from these.
Here we discuss the SRN fluxes evaluated under the assumption that SRNs originate from a single component listed in Table 1. In relation to Equation (2), the models of the SN with a canonical-mass NS, the SN with a high-mass NS, and the failed SN with a BH correspond to (f HNS , f BH ) = (0, 0), (f HNS , f BH ) = (1, 0), and f BH = 1, respectively. Theν e fluxes of SRNs obtained in this way are shown as a function of neutrino energy in Figure 1, where the 90% C.L. observed upper limits on the extraterrestrialν e flux from the latest experimental searches at SK (Abe et al. 2021) and Kam-LAND (Abe et al. 2022) are also shown for comparison. Whereas some models of the failed SN are disfavored by the observed limits, a feasible flux should be obtained by mixing all three components as in Equation (2).
In Figure 2, the integrated SRNν e fluxes with E ν > 17.3 MeV for the single component models are shown together with the best-fit results and 90% C.L. upper limits in Abe et al. (2021). 9 For the SN models with an NS, the integrated flux varies by the EOS and remnant and its dependence is similar to that of the total energy of emitted neutrinos in Table 1. On the other hand, the failed SN models have a larger flux in E ν > 17.3 MeV for their total emission energy because their average neutrino energy is higher than that of the SN models with an NS. Note that the fluxes of the present theoretical models are larger than those of Nakazato et al. (2015) (see also Figure 29 of Abe et al. 2021) because the conversion factor fromρ * (z) to R CC (z) is larger.
The event rate spectra at an SK size water detector (22.5 kton) per year are shown in Figure 3 for the single component models. If the same EOS and neutrino mass 9 In Abe et al. (2021), the best-fit values with ±1σ and the corresponding upper limits were evaluated assuming normalized spectra taken from several theoretical models. While we show in Figure 2 the best-fit result associated with the models of Nakazato et al. (2015), 1.4 +1.0 −0.9 cm −2 sec −1 (same for both of the NH and IH cases), this value is insensitive to the normalized spectrum provided that realistic models are adopted.
hierarchy are assumed, the SN models with a high-mass NS show a higher event rate than the SN models with a canonical-mass NS, while their normalized spectra look similar. This is again consistent with the strong NS mass dependence of the total neutrino energy and the weak NS mass dependence of the average neutrino energy as shown in Table 1. However, the contribution from the failed SN is not so simple. In the case of the Shen EOS models, the failed SN leads to a much higher SRN event rate than the SN with an NS. The Shen EOS has a higher maximum mass of NSs and a larger NS radius. Therefore, the amount of neutrinos emitted from a BH formation is larger and that from an NS formation is smaller. In contrast, the LS220 EOS models show the opposite trend; the failed SN model has a lower event rate due to the low maximum mass of NSs (see Figure 13 of Schneider et al. 2020). Whether the NS or BH leads to a higher SRN event rate depends on combination of the NS radius and maximum mass. As for the failed SN models, the difference between the cases with NH and IH is larger. This is because neutrinos are mainly emitted by the accretion of matter in the case of the failed SN and, thereby,ν e and ν x have different spectra. In this process, the total energy is larger forν e and the average energy is higher for ν x (Table 1).

Experimental sensitivity to the SRN flux
Here we combine the contributions of all three components. Then, dependence on f BH and f HNS is investigated with the sensitivity study. In Figure 4, the 2σ C.L. expected sensitivities on theν e flux at SK-Gd over 10 yr, HK over 3 yr, and HK over 5 yr are shown on top of the integrated SRNν e flux for different energy regions, (1) 13.3 < E ν < 31.3 MeV and (2) 17.3 < E ν < 31.3 MeV, with NH. Figure 5 shows the corresponding distributions with IH. We can recognize that the sensitivity is improved by lowering the energy threshold to 13.3 MeV in some cases. This is thanks to the latest extensive study on the background in the  Table 1. The left and right panels show the results for NH and IH, respectively. The orange solid lines and shaded regions represent the best-fit values and their ±1σ uncertainties and the red dashed-dotted lines the 90% observed upper limits from the latest SK analysis (Abe et al. 2021) associated with the spectral models in Nakazato et al. (2015). Note that the SRN flux for each component may change for various assumptions such as Mmin and BH formation model as discussed in Section 4.3. low energy regime as described in Section 3. Nevertheless, for some parameter sets, the sensitivity is better for 17.3 < E ν < 31.3 MeV. This is likely for the cases with a large f BH because the failed SN models have relatively large fluxes in E ν > 17.3 MeV due to their hard spectra. Therefore, in the following, we plot the projected sensitivity by combining the sensitivities for 13.3 < E ν < 31.3 MeV and 17.3 < E ν < 31.3 MeV. In Figure 6, the detectable regions at SK-Gd over 10 yr are shown for different EOS models. We find that the region at a higher f HNS is more likely to be detected. This is because a larger amount of neutrinos are emitted from a high-mass NS than a canonical-mass NS for every EOS model while their normalized SRN spectra look similar (Figure 3). In contrast, dependence on f BH is different among the EOS models. Since the SRN flux of the failed SN with a BH is lower than that of the SN with a high-mass NS for the case of the LS220 EOS (Figure 2)  not always enhance the possibility of the SRN detection; it depends on the EOS. In any case, whether SK-Gd can reach the 2σ level over 10 yr depends on f BH and f HNS . Thus, within 10 yr, SK-Gd would provide some implications about the distribution of compact remnants formed by stellar core collapse.
We now turn to the SRN observation with HK. Figure 7 shows the detectable regions at HK over 10 yr for different C.L., EOS models, and neutrino mass hierarchies. The results show that the whole parameter space could be tested at a 2σ level for all cases and a large region of the parameter space is covered at a 3σ level. A higher significance could be achieved for only some cases. Note that our assumptions for the HK study still contains a lot more uncertain factors than those for SK-Gd.

Discussion
The gravitational energy released in a form of neutrinos during the failed process strongly depends on the progenitor structure (O'Connor & Ott 2011;Horiuchi et al. 2018;Schneider et al. 2020;Kresse et al. 2021), thus constraints on the fractional parameters shown in Section 4.2 may change for different models of the failed SN forming a BH. In Kresse et al. (2021), the failed SN case that forms a BH with a significant delay ( 1 sec) is studied in which about twice larger neutrino emission energy can be expected than the fast BH formation case. We therefore double the neutrino emission amount of our BH model and study contours with this assumption to investigate systematic uncertainty for the delayed BH formation. The results for HK over 10 yr are shown in Figure 8. In general, as expected, f BH can be constrained more strictly for the increased neutrino     amount, which is more obvious for the Shen EOS whose neutrino emission from BH is greater (see Figure 3). Our SRN flux calculation is based on the minimum mass of progenitors being M min = 6M as described in Section 2.1, while larger values are used in many other studies. In order to investigate the impact of M min on the current study, we calculate the SRN flux with M min = 8M and draw contours on the fractional parameters in Figure 9 for the case of HK over 10 yr. The significance is lower for all cases because of the smaller number of progenitors that contribute to neutrino emission, while the large parameter space can still be tested at a 2σ level.
We showed the SRN flux for the different EOS and neutrino mass hierarchy cases, while these could be constrained by external studies. The constraints on the EOS are provided by the mass measurements of heavy NSs (Antoniadis et al. 2013;Arzoumanian et al. 2018;Cromartie et al. 2020;Romani et al. 2021), the NS radius from the LIGO-Virgo detection of GW170817 (Abbott et al. 2018), and the NS mass and radius from the NICER observation (Riley et al. 2019;Miller et al. 2019Miller et al. , 2021. As shown in Figure 1 of Nakazato et al. (2022), the Togashi EOS is more favored than the other EOSs; the LS220 EOS has the maximum mass of cold NSs (2.04M ) that is lower than the mass of heavy NSs already discovered and the NS radius of the Shen EOS is larger than the estimation from the gravitational wave data. Further details on the EOS and its constraints are seen in Nakazato et al. (2022). NH is currently favored by the neutrino oscillation measurements (see Jiang et al. (2019) for example), while the statistical significance is not enough to make a conclusion. Planned oscillation experiments are expected to confirm this with a higher confidence (Abe et al. 2018;Abi et al. 2020;Aiello et al. 2022).
In the present study, we assume that the fraction of failed SNe (f BH ) is independent of the redshift. However, it would not be the case because the fate of massive stars is determined by metallicity. The metallicity distribution function of progenitors at the selected redshift can also be derived from the data of Illustris-1 result (Nelson et al. 2015), which is adopted for the cosmic star formation rate in our study, and the redshift dependence of f BH will be a subject of future investigations.
Our sensitivity studies are based on realistic assumptions extrapolated from the SK results, but further ef-forts could enhance the detector's power to test the parameter space. For example, reducing the atmospheric NCQE background would improve the sensitivity dramatically at lower energies. This could be achieved by a more precise reconstruction of the neutron vertex. In addition, external data about neutrino interaction cross sections and isotope productions in a muon spallation would help to predict backgrounds more precisely. This becomes more important especially at larger volume detectors such as HK.

CONCLUSION
In this paper, we focused on SRNs as a probe to investigate the fate of stellar core collapse. We prepared core-collapse SN models for three different fates: the SN with a canonical-mass NS, the SN with a high-mass NS, and the failed SN with a BH, for three different nuclear EOSs: Togashi, LS220, and Shen. These models show unique shapes in the resulting neutrino spectrum. We introduced two parameters, f BH (the fraction of failed SNe to the total number of stellar core-collapse events) and f HNS (the fraction of SNe with a high-mass NS to the total number of SNe with an NS), to mix these three fate components to form the SRN flux. We performed the sensitivity study at the next-generation water-based Cherenkov detector, SK-Gd and HK, as extrapolating from the latest SK analysis to achieve realistic assumptions. As a result, we found that the sensitivity depends mainly on f HNS if the delayed BH formation cases are neglected and the EOS has a low maximum mass of NSs as is the case for the LS220 EOS. In contrast, for the EOS with a higher maximum mass such as the Shen EOS, the impact of f BH gets stronger. We also showed that the analysis in different energy ranges would constrain different areas in a f BH -f HNS space. The sensitivity with HK over 10 yr shows that most of the area in this space could be tested at a 3σ level. The results of this paper help to provide implications about the fate of stellar core collapse in compound with other studies such as EOS and neutrino mass hierarchy.