Relic neutrino background from cosmological supernovae

Present and future observations of supernova relic neutrinos (SRNs), i.e., a cosmological neutrino background from past core-collapse supernova explosions, potentially give us useful information concerning various fields of astrophysics, cosmology and particle physics. We review recent progress of theoretical and observational studies of SRNs, particularly focusing on the detectability, and also on implications for cosmic star formation history and neutrino physics.


Introduction
A core-collapse supernova explosion is one of the most spectacular events in astrophysics, and it attracts a great deal of attention from many physicists and astronomers. It also produces a number of neutrinos and 99% of its gravitational binding energy is transformed into neutrinos; detection of the galactic supernova neutrino burst by ground-based large waterČerenkov detectors, such as Super-Kamiokande (SK) and Sudbury Neutrino Observatory (SNO), would provide valuable information on the nature of neutrinos as well as supernova physics. In addition, because supernova explosions have occurred very commonly in both the past and present universe, tracing the cosmic star-formation rate (SFR), they should have emitted a great number of neutrinos, which now make a diffuse background, i.e., supernova relic neutrinos (SRNs). The physics involved in SRN ranges quite widely-from cosmic SFR and supernova physics to neutrino properties as elementary particles. Therefore, detecting SRNs or even setting limits on their flux can give us quite useful and unique implications for various fields of astrophysics, cosmology and particle physics. Detectability of SRNs in various detectors and its implications have been discussed in many theoretical papers [1]- [19] from various points of view, which we discuss in detail in the following part of the present paper.
Neutrino properties can also be probed by the SRN observation in principle. It is pointed out that the resulting SRN signal depends on neutrino oscillation models [14]; particularly, if the mass hierarchy is inverted, i.e., the first mass eigenstate that most strongly couples to electron flavour is heavier than the third state, m 1 > m 3 , and the value of θ 13 , which has not yet been well constrained, is sufficiently large, then the flavour conversion during propagation inside the supernova envelope could considerably change the neutrino signal. Another interesting possibility is that a stringent constraint can be set on neutrino decay models from the SRN observation as first pointed out by Ando [15] and subsequently studied by Fogli et al [18]. Non-radiative neutrino decay can be induced by the interaction between neutrinos and massless or very light particles such as Majoron. The strongest lower limit to the neutrino lifetime-to-mass ratio is obtained from the solar neutrino observation [39]- [41] and meson decay [42]- [44] to be τ/m 10 −4 s eV −1 . SRN observation, on the other hand, is sensitive to ∼10 10 s eV −1 , which is many orders of magnitude larger than the current limit. This is because of the much longer baseline between the Earth and cosmological supernovae, from which SRNs are emitted, compared to 1 AU in the case of solar neutrino observations. Neutrino decay could also change the detected signal from high-energy astrophysical objects [45] or the galactic supernova explosions [46]- [48], and could alter usual discussions on the early universe and structure formation [49] as well as on supernova coolings [50]- [52]; but we stress that, in principle, we can obtain the most stringent constraints on the neutrino lifetime from the SRN observations. This paper is organized as follows. We first introduce the formulation for calculating the SRN flux in section 2. Models of cosmic SFR, supernova neutrino spectrum calculated by various groups are given in the same section. Neutrino oscillation inside the supernova envelope is briefly described there. SRN flux and event rate at ground-based detectors calculated with our reference models are given in section 3. Section 4 is devoted to a detailed discussion on background events against the SRN detection and detectability at various detectors. Current observational limit by SK, which is obtained by a statistical argument including background events, is also briefly summarized in the same section. We discuss implications of the SRN observation for cosmic SFR and neutrino properties such as oscillation and decay in sections 5 and 6, respectively. Finally, we conclude this paper by giving a brief summary in section 7. Throughout this paper, we only consider electron anti-neutrinos (ν e ) at detectors because this kind of flavour is most easily detected by the waterČerenkov detectors, on which we mainly focus.

Formulation
The present number density of SRNs (ν e ), whose energy is in the interval E ν ∼ E ν + dE ν , emitted in the redshift interval z ∼ z + dz, is given by where E ν = (1 + z)E ν is the energy of neutrinos at redshift z, which is now observed as E ν ; R SN (z) represents the supernova rate per comoving volume at z, and hence the factor (1 + z) 3 should be multiplied to obtain the rate per physical volume at that time; dN ν /dE ν is the number spectrum of neutrinos emitted by one supernova explosion; and the factor (1 + z) −3 comes from the expansion of the universe. The Friedmann equation gives the relation between t and z as and we adopt the standard CDM cosmology ( m = 0.3, = 0.7, and H 0 = 70h 70 km s −1 Mpc −1 ). 3 We now obtain the differential number flux of SRNs, dF ν /dE ν , using the relation dF ν /dE ν = c dn ν /dE ν : where we assume that gravitational collapses begin at the redshift z max = 5.

Models for cosmic SFR
As our reference model for the SFR, we adopt a model that is based on recent progressive results of rest-frame UV, NIR Hα, and FIR/sub-millimetre observations; a simple functional form for the SFR per unit comoving volume is given as [53] ψ * (z) = 0.32f * h 70 exp(3.4z) exp(3.8z) + 45 where f * is a factor of order unity, which we illustrate below. Figure 1 shows the SFR ψ * (z) with the various data points from rest-frame UV [22,23,30], Hα line [33]- [35], and FIR/submillimetre [31,32] observations; these data points are not corrected for dust extinction. In the local universe, all studies show that the comoving SFR monotonically increases with z out to a redshift of at least 1. Although there is such a general observational tendency, even the local (z = 0) SFR density is far from precise determination; it ranges fairly widely as ψ * (0) = (0.5-2.9)×10 −2 h 70 M year −1 Mpc −3 [54]. For this reason, we introduce the correction factor f * . Our reference model (4) with f * = 1 is consistent with mildly dust-corrected UV data at low redshift; on the other hand, it may underestimate the results of the other wave band observations. In fact, it predicts the local SFR value of 0.7 × 10 −2 h 70 M year −1 Mpc −3 , which is close to the lower limit of the estimation by Baldry and Glazebrook [54], but we stress that the SRN flux given below using this model can be simply applied to the other cases if we adjust the correction factor f * , and would be quite general. Calculations by Strigari et al [16] are based on the SFR model giving a local value of 1.6 × 10 −2 M year −1 Mpc −3 , which corresponds to f * = 2.3 in our notation. Although the SFR-z relation generally tends to increase from z = 0 to ∼1, behaviours at the higher redshift region z > 1 are not clear at all. Ando et al [12] also investigated the dependence on the several adopted SFR models, which were only different at high-redshift regions (z 1.5); our reference model (4) was referred to as the 'SF1' model there. Cosmic SFR as a function of redshift. Data points are given by restframe UV ( ) [22,23,30], NIR Hα (×) [33]- [35], and FIR/sub-millimetre ( ) [31,32] observations. The solid line represents our reference model given by equation (4) with f * = 1. The standard CDM cosmology is adopted They showed that the SRN flux at E ν > 10 MeV is highly insensitive to the difference among the SFR models (owing to the energy redshift, as discussed in section 3.1). We obtain the supernova rate (R SN (z)) from the SFR by assuming the Salpeter IMF (φ(m) ∝ m −2.35 ) with a lower cutoff around 0.5M , and that all stars with M > 8M explode as core-collapse supernovae, i.e., Here we assume that the IMF does not change with time, which may be a good approximation provided there are no significant correlations between the IMF and the environment in which stars are born; extant evidence seems to argue against such correlations over the redshift range of interest (z 2) [55]. The resulting local supernova rate evaluated with f * = 1 agrees within errors with the observed value of R SN (0) = (1.2 ± 0.4) × 10 −4 h 3 70 year −1 Mpc −3 (see e.g., [56] and references therein). In fact, the totally time-integrated neutrino spectrum from massive stars ( 30M ) could be very different from the models that we use (and give in the next subsection), possibly because of, e.g., black hole formation. However, the conversion factor appearing in equation (5) is highly insensitive to the upper limit of the integral in the numerator; for instance, if we change the upper limit of the numerator to 25M , the factor becomes 0.010M −1 , which is only slightly different from the value in equation (5).

Neutrino spectrum from supernova explosions
For the neutrino spectrum from each supernova, we adopt three reference models by different groups, i.e., simulations by the Lawrence Livermore (LL) group [57] and Thompson, Burrows and Pinto [58] (hereafter TBP), and the MC study of spectral formation by Keil, Raffelt and Janka [59] (hereafter KRJ). In this field, however, the most serious problem is that the recent sophisticated hydrodynamic simulations have not obtained the supernova explosion itself; the shock wave cannot penetrate the entire core. Therefore, many points still remain controversial, e.g., the average energy ratio among neutrinos of different flavours, or how the gravitational binding energy is distributed to each flavour. All these problems are quite serious for our estimation, since the binding energy released asν e changes the normalization of the SRN flux, and the average energy affects the SRN spectral shape. Traditionally, neutrino spectrum is assumed to have a Fermi-Dirac spectral shape with Tν e 5 MeV and T ν x 8 MeV, where ν x represents non-electron neutrinos and anti-neutrinos, as adopted in many studies including recent ones [7]- [11], [13,16,19]. This approximation is roughly consistent with the LL model, although slight difference exists at both high-and low-energy regions [12,57]. The other two models (TBP and KRJ) that are more sophisticated, however, are not consistent with such simple treatment at all as we describe below or as already discussed in [17]. For all these reasons stated above, we believe that these three models from different groups will be complementary. The numerical simulation by the LL group [57] is considered to be the most appropriate for our estimation, because it is the only model that succeeded in obtaining a robust explosion and in calculating the neutrino spectrum during the entire burst (∼15 s). According to their calculation, the average energy difference betweenν e and ν x was rather large and the complete equipartition of the binding energy was realized L ν e = Lν e = L ν x , where L ν α represents the released gravitational energy as α-flavour neutrinos. The neutrino spectrum obtained by their simulation is well fitted by a simple formula, which was originally given by KRJ as whereĒ ν is the average energy; the values of the fitting parameters for theν e and ν x spectrum are summarized in table 1.
Although the LL group succeeded in obtaining a robust explosion, their result has recently been criticized because it lacked many relevant neutrino processes that are now recognized as important. Thus, we adopt the recent result of another hydrodynamic simulation, the TBP one, which included all the relevant neutrino processes, such as neutrino bremsstrahlung and neutrino-nucleon scattering with nucleon recoil. Their calculation obtained no explosion, and the neutrino spectrum ends at 0.25 s after core bounce. In the strict sense, we cannot use their result as our reference model because the fully time-integrated neutrino spectrum is definitely necessary in our estimate. However, we adopt their result in order to confirm the effects of recent sophisticated treatments of neutrino processes in the supernova core on the SRN spectrum. The TBP calculations include three progenitor mass models, i.e., 11, 15 and 20M ; all of these models are well-fitted by equation (6), and the fitting parameters are summarized in table 1. The average energy for bothν e and ν x is much smaller than that by the LL calculation. Although we do not show this in table 1, it was also found that at least for the early phase of the core-collapse, the complete equipartition of the gravitational binding energy for each flavour was not realized. However, it is quite unknown whether these trends hold during the entire burst. In this study, we adopt the average energy given in table 1 as our reference model, while we assume perfect equipartition between flavours, i.e., Lν e = L ν x = 5.0 × 10 52 erg.
In addition, we also use the model by KRJ. Their calculation did not couple with the hydrodynamics, but it focused on the spectral formation of neutrinos of each flavour using an MC simulation. Therefore, the static model was assumed as a background of neutrino radiation, and we use their 'accretion phase model II', in which the neutrino transfer was solved in the background of a 150 ms postbounce model by way of a general relativistic simulation. The fitting parameters for their MC simulation are also summarized in table 1. Unlike the previous two calculations, their result clearly shows that the average energy of ν x is very close to that ofν e . It also indicates that the equipartition among each flavour was not realized, but rather L ν e Lν e 2L ν x . Since the totally time-integrated neutrino flux is unknown from such temporary information, in this case also, we assume perfect equipartition, Lν e = L ν x = 5.0 × 10 52 erg, as well as that the average energies are the same as those in table 1.

Neutrino spectrum after neutrino oscillation
The originalν e spectrum is different from what we observe asν e at Earth, owing to the effect of neutrino oscillation. Since the specific flavour neutrinos are not mass eigenstates, they mix with other flavour neutrinos during their propagation. The behaviour of flavour conversion inside the supernova envelope is well understood, because the relevant mixing angles and mass square differences are fairly well determined by recent solar [60,61], atmospheric [62], and reactor neutrino experiments [63]. The remaining ambiguities concerning the neutrino oscillation parameters are the value of θ 13 , which is only weakly constrained (sin 2 θ 13 0.1 [64]), and the type of mass hierarchy, i.e., normal (m 1 m 3 ) or inverted (m 1 m 3 ). We first discuss the case of normal mass hierarchy as our standard model; in this case, the value of θ 13 is irrelevant. The case of inverted mass hierarchy is addressed in section 6.1. In addition, other exotic mechanisms, such as resonant spin-flavour conversion (see [65] and references therein) and neutrino decay [15,18], which possibly change the SRN flux and spectrum, might work in reality, and these topics are also discussed later in sections 6.2 and 6.3, respectively.
The producedν e at the supernova core are coincident with the lightest mass eigenstatē ν 1 owing to the large matter potentials. Since this stateν 1 is the lightest in vacuum also, there are no resonance regions in which one mass eigenstate can change into another state, and thereforeν e at production arrives at the stellar surface asν 1 . Thus, theν e spectrum observed by The same as (a), but indicating contribution from various redshift ranges. LL is adopted as the supernova model. These figures are taken from [17].
the distant detector is where the quantities with superscript 0 represent those at production, U αi is the mixing matrix element between the α-flavour state and ith mass eigenstate, and observationally |U e1 | 2 = 0.7. In other words, 70% of the originalν e survives; on the other hand, the remaining 30% comes from the other component ν x . Therefore, both the originalν e and ν x spectra are necessary for the estimation of the SRN flux and spectrum; since the original ν x spectrum is generally harder than that of the originalν e , as shown in table 1, the flavour mixing is expected to harden the detected SRN spectrum.

Flux of SRNs
The SRN flux can be calculated by equation (3) with our reference models given in section 2. Figure 2(a) shows the SRN flux as a function of neutrino energy for the three supernova models, LL, TBP and KRJ. The flux of atmosphericν e , which becomes a background event for SRN detection, is shown in the same figure [66,67]. The SRN flux peaks at 5 MeV, and around this peak, the TBP model gives the largest SRN flux because the average energy of the originalν e is considerably smaller than in the other two models but the total released energy is assumed to be the same. On the other hand, the model gives a smaller contribution at high-energy regions, E ν > 10 MeV. In contrast, the high-energy tail of the SRN flux with the LL model extends  farther than with the other models, and it gives flux more than an order of magnitude larger at E ν = 60 MeV. This is because the high-energy tail was mainly contributed by the harder component of the original neutrino spectrum; in the case of the LL calculation, the average energy of the harder component ν x is significantly larger than that of the other two calculations, as shown in table 1. We show the values of the SRN flux integrated over the various energy ranges in table 2.
The total flux is expected to be (11-16)f * cm −2 s −1 for our reference models, although this value is quite sensitive to the shape of the assumed SFR, especially at high-z. The energy range in which we are more interested is high-energy regions such as E ν > 19.3 MeV and 11.3 MeV, because as discussed below, the background events are less critical and the reaction cross section increases as ∝E 2 ν . In such a range, the SRN flux is found to be (1. Thus, the uncertainty about the supernova neutrino spectrum and its luminosity gives at least a factor of 2-4 ambiguity compared to the expected SRN flux in the energy region of our interest. Figure 2(b) shows the contribution by supernova neutrinos emitted from various redshift ranges. At the high-energy region E ν > 10 MeV, the dominant flux comes from the local supernovae (0 < z < 1), while the low-energy side is mainly contributed by the high-redshift events (z > 1). This is because the energy of neutrinos that were emitted from a supernova at redshift z is reduced by a factor of (1 + z) −1 reflecting the expansion of the universe, and therefore high-redshift supernovae only contribute to low-energy flux. We also show the energy-integrated flux from each redshift range in table 2 in the case of the LL supernova model. From the table, it is found that in the energy range of our interest, more than 70% of the flux comes from local supernova explosions at z < 1, while the high-redshift (z > 2) supernova contribution is very small.

Event rate at waterČerenkov detectors
The waterČerenkov neutrino detectors have greatly succeeded in probing the properties of neutrinos as elementary particles, such as neutrino oscillation. The SK detector is one of these detectors, and its large fiducial volume (22.5 kton) might enable us to detect the diffuse background of SRNs. Furthermore, much larger waterČerenkov detectors such as HK and UNO   are being planned. SRN detection is most likely with the inverse β-decay reaction with protons in water,ν e p → e + n, and its cross section is precisely understood [68,69]. In our calculation, we use the trigger threshold of SK-I (before the accident). (b) Event rate of SRNs and invisible muon decay products. In both panels, LL is adopted as the original neutrino spectrum. These figures are taken from [12].
MeV. This clearly indicates that if the background events that hinder the detection are negligible, the SK has already reached the required sensitivity for detecting SRNs; with the future HK and UNO, a statistically significant discussion would be possible. This also shows that the current shortage of our knowledge concerning the original supernova neutrino spectrum and luminosity gives at least a factor of 2 (E ν > 10 MeV) to 4 (E ν > 18 MeV) uncertainty to the event rate at the high-energy range (actual detection range). We also summarize the contribution from each redshift range in the same table, especially for the calculation with the LL model. The bulk of the detected events will come from the local universe (z < 1), but the considerable flux is potentially attributed to the range 1 < z < 2.

Background events against detection
In the previous section, we calculated the expected SRN spectrum at the waterČerenkov detectors on the Earth, but the actual detection is quite restricted because of the presence of other background events. In this paper, we follow a detailed consideration of these backgrounds by Ando et al [12]. There are atmospheric and solar neutrinos, anti-neutrinos from nuclear reactors, spallation products induced by cosmic-ray muons and decay products of invisible muons. We show in figures 4(a) and (b) the flux and event rate of SRNs and these background events.
The flux of the atmospheric neutrinos is usually calculated using MC method including various relevant effects (flux of primary cosmic rays, solar modulation, geomagnetic field, interaction of cosmic rays in the air, and so on), and in that simulation one-dimensional approximation is used, i.e., after the interaction of primary cosmic ray particles with air nuclei, all the particles are assumed to be moving along the line of the momentum vector of the primary cosmic ray particles [70]. There are many authors who calculated the atmospheric neutrino flux (see [71]- [73] for recent one-dimensional results). We use in this paper the flux calculated by Gaisser et al [66,67]. More recently, results of three-dimensional flux calculations have been reported by several groups (see [74,75] and references therein). These calculations show slight increase of the flux at low-energy regions, although the flux below 100 MeV is not given. Therefore, it should be noted that there is possibility that the flux of atmospheric neutrinos in the relevant energy regime is higher than that adopted in this paper by about 20-30%. Solar neutrino flux is dominant in the energy range below 19 MeV. We use the flux predicted by the standard solar model (SSM) in figure 4 [76]. Since the solar neutrinos are notν e but ν e , the cross section for them is about two orders of magnitude smaller than that forν e . Furthermore, recoil electrons scattered by solar neutrinos strongly concentrate to the opposite direction of the Sun, in contrast to the isotropic distribution ofν e events. Therefore, the solar neutrinos are an avoidable background, not as critical as other events. At the same energy range corresponding to E ν 19 MeV, there is another serious background that becomes an obstacle for solar neutrino detection, i.e., spallation products induced by cosmic ray muons. The event rate of the spallation background is several hundred per day per 22.5 kton, and it is extremely difficult to reject all of these events. Therefore, the solar neutrino range (E ν 19 MeV), cannot be used as an energy window.
The third background which we must consider is anti-neutrinos from nuclear reactors. In each nuclear reactor, almost all the power comes from the fissions of the four isotopes, 235 U (∼75%), 238 U (∼7%), 239 Pu (∼15%) and 241 Pu (∼3%) [77]. Each isotope produces a unique electron anti-neutrino spectrum through the decay of its fission fragments and their daughters. Theν e spectrum from 235 U, 239 Pu and 241 Pu can be derived using the semi-empirical formula with which we fit the data of detected β-spectrum from fission by thermal neutrons [78,79]. ( 238 U undergoes only fast neutron fission and hence electron spectrum from 238 U cannot be measured by this kind of experiment.) Above 7 MeV, the number of β counts drops dramatically and fitting error becomes large. In addition, with this method, as we determine the maximum β energy and derive the energy distribution below that energy, it is difficult to estimate the errors at the high-energy range. While theν e spectra in [78,79] are given as tables, we use for simplicity a somewhat less accurate analytical approximation given in [80]. As a normalization factor we use energy-integratedν e flux at Kamioka, 1.34 × 10 6 cm −2 s −1 , which is the total flux from various nuclear reactors in Japan and Korea [77].
With the backgrounds that we discussed above and from figure 4(a), we expect the energy window of SRN events to be in the range 19-30 MeV. However, electrons or positrons from invisible muons are the largest background in the energy window from 19 to 60 MeV. This invisible muon event is illustrated as follows. Atmospheric neutrinos produce muons by interaction with the nucleons (both free and bound) in the fiducial volume. If these muons are produced with energies belowČerenkov radiation threshold (kinetic energy less than 53 MeV), then they will not be detected ('invisible muons'), but successively produced electrons and positrons from the muon decay will be visible. Since this muon decay signal will mimic thē ν e p → e + n process in SK, it is difficult to distinguish SRN from these events. The spectrum of this invisible muon events is shown in figure 4(b), compared with the SRN spectrum. Therefore, even at the remaining part of the energy window 19-30 MeV, there is a huge background due to the invisible muons, resulting in a 'no energy window'for the SRN detection at present.

Detectability at pure-waterČerenkov detectors
We expect that SRNs could most likely be detected at waterČerenkov detectors such as SK, because the largest class of volume can be realized, enabling the most statistically significant discussions. The most serious problem is that there is no energy window as we described in the previous subsection. In the energy range 19-30 MeV, the SRN event number is estimated to be N SRN = 0.73f * (V eff 22.5 kton −1 year −1 ) for the LL model, where we define the effective volume of detectors V eff as (fiducial volume)×(time)×(efficiency), which is the relevant quantity representing the detector performance. On the other hand, the background level would be as large as N bg ∼ 3.4 (V eff 22.5 kton −1 year −1 ), contributed mainly by the invisible muon decay products. Therefore, signal-to-noise ratio S/N can be written as plotted in figure 5 as a function of f * for several values of the effective volume. For detectors with the size of SK, one-year observation is found to be insufficient to obtain significant SRN detection. Further data-taking for about 10 years may reach S/N a few, but this number is still insufficient in statistical significance. Future mega-ton class detectors such as HK and UNO, now under planning, potentially give us a considerable number of SRN detection, by data-taking for a couple of years, as can be seen from the upper curve in figure 5.

Detectability at other detectors
In this subsection, we discuss the SRN detectability at SNO (the heavy waterČerenkov detector) and KamLAND (the scintillation detector), which is mainly discussed in [11,12] and [12,16], respectively. An advantage of these detectors is that we are able to identifyν e events using delayed coincidence signals; neutrons produced viaν e d → e + nn (SNO) andν e p → e + n (KamLAND) are tagged by deuterons, Cl or He (SNO) and by protons (KamLAND), resulting in γ-ray cascades, identified with the preceded positron signal. Using this technique, we can remove other backgrounds of non-ν e origin (solar neutrinos, invisible muons and spallation products), opening up an energy window in the range of 10-30 MeV. Unfortunately, because the detector volume is small (both detectors are about 1 kton), the expected SRN rate is quite small, 0.03f * year −1 for SNO and 0.1f * year −1 for KamLAND [12]. More recently, Strigari et al [16] estimated the event rate at KamLAND to be ∼0.4 year −1 , but they used the energy range above 6 MeV, simply because there have been no events seen above this energy at KamLAND [63]. In addition, the detectability at liquid argon detectors has been discussed in [19]; although the detection is still challenging, potential advantage is that these detectors can mainly capture ν e , which is difficult with other detectors.

Current observational limit and future prospects
The most stringent upper limit on the SRN flux is obtained by the observation for 1496 days (4.1 years) at the SK detectors [20]. This limit is obtained by the statistical analysis including the background events from atmospheric neutrinos and invisible muons, and is <1.2 cm −2 s −1 over the energy region of E ν > 19.3 MeV (90% CL). Comparing with the prediction for the same energy range given in table 2, we find that the current SK limit is only about a factor 2.5-8.5 larger than our prediction using the reference model for the cosmic SFR with f * = 1, depending on the adopted original neutrino spectrum. This strong constraint motivated many theoretical studies [13]- [16] and has been translated into constraints on various quantities as we further discuss in the following sections. Although the SK limit was derived by a careful analysis of the spectral data over the energy range E e > 18 MeV, we here show that a very rough statistical argument using equation (8) can fairly well reproduce the same limit. In equation (8), we consider setting a limit at 90% CL, which corresponds to 1.64σ level (or S/N < 1.64), as a result of no SRN detection for V eff = 22.5 kton · 4.1 year observation. Substituting these values into equation (8), it is solved for f * , resulting in the upper limit f * < 2.6 (90% CL). This value is then translated into the flux limit using the relation between the SRN flux and f * (for the LL model) given in table 2, i.e., F ν < 1.2 cm −2 s −1 for E ν > 19.3 MeV (90% CL), which is in remarkable agreement with the actual result reported by the SK group [20].
In the near future, sensitivity of waterČerenkov detectors for the SRN detection would be significantly improved by the promising technique proposed recently [21]. The basic idea is the same as the delayed coincidence technique actually adopted by SNO or KamLAND (see discussions in section 4.2), but GdCl 3 is dissolved into the pure-water of SK (or other future detectors), which enables us to actively identifyν e by capturing neutrons produced by theν e p → e + n reaction. Owing to this proposal, the range 10-30 MeV would be an energy window because we can positively distinguish theν e signal from other backgrounds such as solar neutrinos (ν e ), invisible muon events and spallation products. The neutron capture efficiency by Gd is estimated to be 90% with the proposed 0.2% admixture by mass of GdCl 3 in water, Table 4. Constraints on supernova rate and SFR models. Upper part shows the SK limit on the SRN observation [20], whereas the lower part shows the results of local supernova surveys [56] and observational inference of the local SFR [54].
[56] 0.9-1.9 (0.8-1.6) × 10 −4 - [54] 0.7-4.2 -(0.5-2.9) × 10 −2 and subsequently 8 MeV γ-cascade occurs from the excited Gd. The single-electron energy equivalent to this cascade was found to be 3-8 MeV by a careful simulation [81], and with the trigger threshold adopted in SK-I, only about 50% of such cascades can be detected actually. However, it is expected that SK-III, which will begin operation in mid-2006, will trigger at 100% efficiency above 3 MeV, with good trigger efficiency down to 2.5 MeV [21]. In that case most of the γ-cascades from Gd will be detected with their preceding signal of positrons. From this point on, we assume 100% efficiency; even if we abandon this assumption, it does not affect our physical conclusion, since the relevant quantity representing the detector performance is the effective volume V eff that already includes efficiency.

Constraints from the current observational limit
The current SK upper limit on the SRN flux is already stringent to give some physical or astronomical consequences. In this subsection, we discuss a constraint on the cosmic SFR from the SRN limit; similar arguments have been given in [13]. Comparing with our flux predictions using several supernova models (table 2), the SK limit can be directly translated into a bound on the correction parameter f * to the cosmic SFR introduced in (4), since the redshift dependence is roughly consistent among various observations. Further, the constraint on f * can then be used as that on the local supernova rate and SFR. The results are summarized in  [56], and the observationally inferred SFR ψ * (0) = (0.5-2.9) × 10 −2 M year −1 Mpc −3 [54]. In particular, for the local SFR, the SK limit on the SRN flux may rule out some fraction of the observationally inferred value, if we choose the LL model for the original neutrino spectrum; for the other two supernova models, the SK limit is very close to the current upper bound on the cosmic SFR by observations with the light. Therefore, we stress that the neutrino observation has already reached sensitivity to the cosmic SFR comparable with the usual and traditional approaches using light.

Performance of Gd-loaded detectors
Performance of the proposed Gd-loaded detectors as an SFR probe is of interest to us, and has recently been investigated in detail by Ando [17] using the MC simulation; in this subsection, we briefly introduce his discussion. Although we focus here on how far the SFR can be probed by SRN observation, the uncertainty from the supernova neutrino spectrum would give a fair amount of error. However, this problem can be solved if a supernova explosion occurs in our galaxy; the expected event number is about 5000-10 000 at SK, when a supernova neutrino burst occurs at 10 kpc, and it will enable a statistically significant discussion concerning the neutrino spectrum from supernova explosions. Even if there are no galactic supernovae in the near future, remarkable development of the supernova simulation can be expected with the growth of computational resources and numerical technique. With such developments, the supernova neutrino spectrum and luminosity may be uncovered, and the ambiguity is expected to be reduced significantly. Thus, in this paper we assume that the supernova neutrino spectrum is well understood and that our reference models are fairly good representatives of nature; we analyse the SFR alone with several free parameters. The basic procedure adopted in [17] is as follows. (1) The expected signal (spectrum) at a Gd-loaded detector in the range 10-30 MeV is simulated, assuming that there are no background events. In that process, the SFR given by equation (4) with f * = 1 and the LL model as neutrino spectrum are used for the generation of the SRN signal. (2) Then the SRN spectrum is analysed using the maximum likelihood method with two free parameters of the SFR and a set of the best-fit values for those parameters is obtained; they are concerned with the supernova rate as where R 0 SN represents the local supernova rate and α determines the slope of supernova rate evolution. Although a constant SFR at z > 1 is assumed, it is found that the result would be the same even if this assumption is changed. This is because the bulk of the detected event comes from local supernova as already shown. (3) 10 3 such MC simulations are performed and 10 3 independent sets of best-fit parameters are obtained. Then we discuss the standard deviation of the distributions of such best-fit parameter sets and the implications for the cosmic SFR.
Using this procedure, we discuss the performance of Gd-SK for 5 years, or an effective volume of 22.5 kton · 5 year. As the expected event number is only ∼10, the parameters R 0 SN and α both cannot be well determined at once. Therefore, one of those parameters should be fixed in advance at some value inferred from other observations. Firstly, the value of R 0 SN was fixed to be 1.2 × 10 −4 year −1 Mpc −1 , which was inferred from the local supernova survey [56], and the distribution of the best-fit values of the parameter α was obtained. The result of 10 3 MC simulations are shown in figure 6(a) as a histogram of the distribution of best-fit parameters α (solid histogram). The average value of these 10 3 values for α is found to be 2.7, and the standard deviation is 0.8, i.e., α = 2.7 ± 0.8. Secondly, in turn, the value of α was fixed to be 2.9 in order to obtain the distribution of best-fit values for the local supernova rate R 0 SN from the SRN observation. The result of 10 3 MC generations and analyses in this case is shown in figure 6(b). The average value for R 0 SN is 1.  Figure 7 shows the comparison between the supernova rate model, in which the parameter is inferred from the MC simulations, and the 'true' reference model; the cases of fixed R 0 SN and α are shown in figures 7(a) and (b), respectively. The allowed region at the 1σ level is located between the two dotted curves, while the solid curve represents the reference model, from which MC data were generated. Thus, with the Gd-SK detector we can roughly reproduce the supernova rate profile at z < 1 for 5 years' operation.  [17]. Now, we turn our attention to future mega-ton class detectors such as Gd-HK or Gd-UNO. With these detectors, the effective volume considered, 440 kton · 5 year, is expected to be realized in several years from the start of their operation. Firstly, the same analysis as illustrated above was performed; i.e., one of relevant parameters, α or R 0 SN , was fixed and the dependence on the remaining parameter was investigated. The result of these cases are also shown in figures 6(a) and (b) as dashed histograms, which give α = 2.5 ± 0.2 and R 0 SN = (1.0 ± 0.1) × 10 −4 year −1 Mpc −3 , respectively. The statistical errors are considerably reduced compared with the case of 22.5 kton · 5 year, because of the ∼20 times larger effective volume. Thus, future maga-ton detectors will probably pin down, within 10% statistical error, either the index of supernova rate evolution α or the local supernova rate R 0 SN if the other is known in advance. The dashed curves in figure 7 set the allowed region of the supernova rate at the 1σ level by the considered detectors, well reproducing the assumed model. In principle, we can determine both parameters by SRN observation, because R 0 SN is concerned with the absolute value of the flux alone but α is concerned with both the absolute value and the spectral shape; i.e., these two parameters are not degenerate with respect to each other. Thus, the same procedure was repeated but without fixing the values of α or R 0 SN . The distribution of 10 3 best-fit parameter sets of (α, R 0 SN ) was calculated assuming a detector with an effective volume of 440 kton · 5 year; the mean values and the standard deviations were found to be α = 3.5 ± 1.3 and R 0 SN = (8.8 ± 4.8) × 10 −5 year −1 Mpc −3 . Even though the effective volume is as large as 440 kton · 5 year, it is still insufficient for determining both parameters at once. Figure 7. Supernova rate as a function of redshift. In both panels, solid curves represent the reference model, from which MC data were generated. (a) The allowed region at the 1σ level, concerning the fitting parameter α with fixed R 0 SN , is shown as the area between the two dotted curves for an effective volume of 22.5 kton · 5 year and as the area between the two dashed curves for an effective volume of 440 kton · 5 year. (b) The same as (a) but for fitting parameter R 0 SN with fixed α. This figure is taken from [17].
Throughout the above arguments or those in [17], the reference SFR model with f * = 1 was assumed and they seem to be dependent on the adopted values of f * . However, the results can easily be applied to the case with other values of f * . This is because the statistical errors are inversely proportinal to the square root of the event number, which is simply proportional to f * V eff in our arguments; i.e., δα/ α ,δR 0 SN / R 0 SN ∝ (f * V eff ) −1/2 . In consequence, we find that the discussions given above and in [17] are quite general and can be easily modified in the other cases of f * using the following relations: As a next step, it would be useful to investigate the possibility of how far (to which z) we can probe the cosmic SFR by the SRN observations. As we have already shown, the main contribution to the SRN event rate at 10-30 MeV comes from the low-redshift region 0 < z < 1. Signals from further high-redshift universe would become enhanced if we could reduce the lower energy threshold E th of the SRN detection range. The value of E th is restricted to 10 MeV because at energy regions lower than this, there is a large background of reactor neutrinos; its removal is impossible with the current detection methods. Since the SK and HK detectors are (will be) located at Kamioka in Japan, they are seriously affected by background neutrinos from many nuclear reactors. If some large-volume detectors were built at a location free from such background, the lower threshold energy could be reduced, enabling us to probe the high-redshift supernova rate. Ando [17] also discussed the performance of future mega-ton class of detectors as a function of the value of E th and found that the behaviour of SFR at 1 < z < 2 affects more and more the detected event number at E th < E e < 30 MeV, as we reduce the threshold energy.

Inverted mass hierarchy
Throughout the above discussions, we have assumed normal hierarchy of neutrino masses (m 1 m 3 ). However, the case of inverted mass hierarchy has not yet been experimentally excluded, and we explore this possibility in this subsection following the discussion given by Ando and Sato [14]. In this case, flavour conversions inside the supernova envelope change dramatically, compared with the normal mass hierarchy already discussed in section 2.4. Sincē ν 3 is the lightest,ν e are created asν 3 , owing to a large matter potential. In that case, it is well known that at a so-called resonance point, there occurs a level crossing betweenν 1 andν 3 (for a more detailed discussion, see, e.g., [82]). At this resonance point, completeν 1 ↔ν 3 conversion occurs when the so-called adiabaticity parameter is sufficiently small compared to unity (it is said that resonance is 'non-adiabatic'), while conversion never occurs when it is large (adiabatic resonance). The adiabaticity parameter γ is quite sensitive to the value of θ 13 , i.e., γ ∝ sin 2 2θ 13 ; when sin 2 2θ 13 10 −3 (sin 2 2θ 13 10 −5 ), the resonance is known to be completely adiabatic (non-adiabatic) [82]. When the resonance is completely non-adiabatic (because of small θ 13 ), the situation is the same as in the case of normal mass hierarchy already discussed in section 2.4 (becauseν e at production becomesν 1 at the stellar surface), and theν e spectrum after oscillation is represented by equation (7). On the other hand, adiabatic resonance (due to large θ 13 ) forcesν e at production becomeν 3 when they escape from the stellar surface, and therefore the observed ν e spectrum is given by The second equality follows from the fact that the value of |U e3 | 2 is constrained to be much smaller than unity from reactor experiments [64]. Thus, equation (11) indicates that complete conversion takes place betweenν e and ν x . When the value of θ 13 is large enough to induce adiabatic resonance (sin 2 2θ 13 10 −3 ), the obtained SRN flux and spectrum should be very different from those obtained in sections 3.1 and 3.2. The SRN flux and event rate in this case were calculated with equations (3) and (11), and the results are summarized in tables 5 and 6. The values (with the LL model) shown in this table are consistent with the calculation by [14], in which numerically calculated conversion probabilities were adopted with some specific oscillation parameter sets (which include a model with inverted mass hierarchy and sin 2 2θ 13 = 0.04), as well as realistic stellar density profiles; they also included the shock wave propagation and the Earth matter effect in their calculations, but both were found to affect only by a few percent and thus irrelevant. The total flux becomes (9.4-14)f * cm −2 s −1 , somewhat smaller than the values given in table 2, because the total flux is dominated by the low-energy region. The fluxes at E ν > 19.3 MeV are enhanced to be (0. .94)f * cm −2 s −1 , but this is still below the current   (10), because they are inversely proportional to the square root of the event number.

Resonant spin-flavour conversion
If neutrinos have non-zero magnetic moment as large as 10 −12 µ B , where µ B is the Bohr magneton, it potentially changes a supernova neutrino signal owing to an additional effect of the flavour conversions. In particular, if neutrinos are Majorana particles, the interaction between the Majorana magnetic moment and supernova magnetic fields induces a spin-flavour conversion (e.g.,ν e ↔ν µ,τ , ν e ↔ν µ,τ ) resonantly (see [65] and references therein). This mechanism could potentially give quite a characteristic supernova neutrino signal at detectors, and is also expected to affect the SRN spectrum significantly. However, it is still premature to estimate the SRN flux including the spin-flavour conversions for several reasons. Firstly, the shock wave propagation can change the magnetic field structure as well as the density profile of supernova progenitors, both of which are essential for calculating flavour-conversion probabilities [83]. Because there is no reliable supernova simulation that succeeded in pushing the shock wave outside the core, we are not even at the stage to start the calculation. Furthermore, calculating how the propagating shock wave changes the magnetic field structure would be quite a difficult task. The second reason is that the spin-flavour conversion probabilities strongly depend on the metallicity of progenitor stars [84]. Because SRN is the accumulation of neutrinos from past supernova explosions, it should be quite natural that the poor metal stars give some contribution to the SRN flux. In consequence, it would be difficult to obtain some implication for the neutrino magnetic moment or supernova magnetic fields from the SRN detection at present; instead, a future galactic supernova neutrino burst might give some clues by its time profile or spectrum [65,85].

Decaying neutrinos
Non-radiative neutrino decay, which is not satisfactorily constrained, potentially and significantly changes the SRN flux and spectrum. The most stringent lower limit on the neutrino lifetimeto-mass ratio comes from the solar neutrino observations [39]- [41] as well as the meson decay experiments [42]- [44], which is τ/m 10 −4 s eV −1 . Since this limit is still very weak, neutrino decay could also change the detected signal from high-energy astrophysical objects [45] or the galactic supernova explosions [46]- [48], and could alter usual discussions on the early universe and structure formation [49] as well as on supernova coolings [50]- [52]. In this subsection, we discuss the potential consequence of the neutrino decay on the SRN flux, with which we can obtain the most stringent limit on the neutrino lifetime in principle, as pointed out by Ando [15].
Since the flux estimation of SRNs including the neutrino decay is somewhat complicated, we rather follow a simple and intuitive argument given in [15]; exact formulation was recently derived by Fogli et al [18]. Ando [15] assumed the following conditions: (i) daughter neutrinos are active species; (ii) neutrino mass spectrum is quasi-degenerate (m 1 ≈ m 2 ≈ m 3 ); and neutrino mass hierarchy is normal. This is because the SRN flux would be most strongly enhanced owing to the decay and it might give a larger flux than the current observational limit by SK [20]; in that case the flux limit can be translated into the limit on the neutrino lifetime. He also restricted his discussion only to the helicity-conserving mode ofν 3 →ν 1 andν 2 →ν 1 for simplicity, characterized by lifetimes τ 3 and τ 2 of each mode; in models discussed in [86]- [88], the above condition (ii) strongly suppresses the helicity-flip mode. In the actual calculation, instead of these lifetimes, he introduced the 'decay redshift' z d i (i = 2, 3) of the mass eigenstateν i as two free parameters. If the source redshift z is larger than the decay redshift z d i , all the neutrinosν i were assumed to decay; on the other hand, if z < z d i , completely survive. The relation between z d i and τ i /m can be written as where we assume E ν = 10 MeV as a typical neutrino energy when evaluating the lifetime τ i from z d i .   figure 8 shows the SRN flux for z d 2 = z d 3 = z max (=5), which represents the same flux as that without the neutrino decay. Then, the decay of the heaviest mass eigenstateν 3 was included by reducing the value of z d 3 , but keepingν 2 stable. When z d 3 = 1.0 (dotted curve), the SRN flux at the low-energy region (E ν 35 MeV) deviates from the nodecay model. This is because the neutrinos from supernovae at redshift larger than z d 3 = 1.0 are affected by theν 3 →ν 1 decay and it results in the increase ofν e . Since the neutrino energies are redshifted by a factor of (1 + z) −1 owing to an expansion of the universe, the decay effect can be seen at low energy alone. When the value of z d 3 is reduced to 10 −2 , the neutrinos even from the nearby sources are influenced by theν 3 →ν 1 decay, resulting in a deviation over the entire energy range as shown by the long-dashed curve in figure 8. If theν 2 →ν 1 decay is added, it further enhances the SRN flux. In table 7, the SRN flux integrated over the energy range of E ν > 19.3 MeV is summarized, for the each decay models. In the second column, the lifetimeto-mass ratio is indicated, which corresponds to each decay redshift, which is obtained using equation (12). The ratio between the predicted flux and the 90% CL upper limit given by the SK observation [20] is also shown in the fourth column. This result shows that several decaying models may have already been excluded or severely constrained by the current SRN limit by SK. It also proves that the SRN observation can potentially give us a lower limit on the neutrino lifetime as strong as ∼10 10 s eV −1 , which is many orders of magnitude stronger than that by the solar neutrino observations. In addition, future detectors such as the HK and UNO (may be loaded with Gd) is expected to greatly improve our knowledge of the SRN spectral shape as well as its flux; if a number of data were actually acquired by such detectors, the most general and model-independent discussions concerning the neutrino decay would become accessible.

Conclusions
As stressed several times in this paper, we are in an exciting era of the SRN observation. This is because the current SK upper limit [20] is just above several theoretical predictions using realistic models for the cosmic SFR and supernova neutrino spectrum, and also because of the promising technique proposed in [21] could greatly improve performance of the presently working or planned future detectors, especially for detecting SRNs. Having arrived at such a situation, we believe that it is time to review recent progress of theoretical and observational researches concerning SRNs from various points of view. Many physical and astronomical consequences can be derived from the constraints on the SRNs as we have reviewed throughout this paper. The basics of SRN calculation were given in sections 2 and 3. Models involved in the SRN calculations are those of the cosmic SFR and original supernova neutrino spectrum. Although, the SFR-z relation has been intensively studied in recent years, there remains a fair amount of uncertainty; in fact, the best estimate of the local SFR density ranges fairly widely as (0.5-2.9) × 10 −2 h 70 year −1 Mpc −3 [54]. In order to take this uncertainty into account, we have introduced the correction factor f * in equation (4), for which a wide range (0.7-4.2) is still allowed from the SFR observation. As for the original supernova spectrum, contrary to the traditional approach using the Fermi-Dirac distribution, we followed the approach taken by our group [12,17], i.e., adopting the results of numerical simulations. In particular, following [17], we used three neutrino spectra by different groups, i.e., LL [57], TBP [58] and KRJ [59], and performed comparison among models. Recent progress in neutrino experiments has proven that neutrinos are massive and mix among different flavours. This effect was also taken into account appropriately. As a result, it was found that the uncertainty concerning the original neutrino spectrum gives difference in the resulting SRN flux by at least a factor of 2-4. The expected event rate at a waterČerenkov detector with a fiducial volume of 22.5 kton was found to be 0.25-1f * year −1 (for E e > 19.3 MeV) and 1-2f * year −1 (for E e > 10 MeV), depending on these supernova models (see table 3).
Besides the flux estimation, a careful discussion about background events is definitely necessary in order to investigate the detectability. This has been thoroughly argued in [12] and we followed them in section 4. The most serious background comes from decay products of invisible muons, which are produced almost at rest inside the detector and that is why 'invisible.' Because of the background, the detection is severely restricted; for a pure-waterČerenkov detector with the size of SK, it would take about (may be more than) 10 years in order to reach a signal-tonoise ratio with a small value (see figure 5). We also discussed the performance of other detectors such as SNO and KamLAND in the same section, by mainly following [12,16]. The current observational upper limit by SK [20] and the proposed technique that potentially improves the performance of waterČerenkov detectors [21] were briefly introduced there. In addition, we showed that the SK limit on the SRN flux can be fairly well reproduced by a simple statistical argument using equation (8).
Section 5 was devoted to the discussion on the current status and future prospects for obtaining limits on the cosmic SFR from the SRN observations. We first showed that the current SRN limit from SK [20] already sets a stringent limit on the local SFR value, as discussed in [13]. Surprisingly, it might rule out a part of local SFR value allowed by the current observations using the light as shown in table 4. Performance of future detectors loaded with Gd as proposed in [21] was investigated in detail using the MC simulation, following the discussion in [17]. The distribution of the best-fit values of two parameters concerning SFR (9) was obtained by the generated 10 3 MC data and by accompanying likelihood analyses. If one of the two free parameters was fixed somehow in advance using other observations, then the SRN observation can constrain the remaining parameter by accuracy of ∼0.3 (f * V eff /22.5 kton · 5 year) −1/2 ; this number would soon become quite significant if the future mega-ton class detectors such as HK or UNO loaded with Gd started data-taking.
Finally, in section 6, SRN constraints on particle physics models of neutrinos were discussed. Neutrino oscillations in the case of inverted mass hierarchy [14], resonant spinflavour conversions induced by the neutrino magnetic moment (see [65] and references therein) and neutrino decay [15,18] are expected to considerably change the SRN signal at the detectors. In particular, for the neutrino decay model, the SRN observation potentially gives the most stringent constraint on the neutrino lifetime, compared with a galactic supernova [46]- [48], high-energy astrophysical objects [45] as well as the most reliable (but weak) limit at present by the solar neutrino observations [39]- [41].