Towards a complete reconstruction of supernova neutrino spectra in future large liquid-scintillator detectors

In this paper, we show how to carry out a relatively more realistic and complete reconstruction of supernova neutrino spectra in the future large liquid-scintillator detectors, by implementing the method of singular value decomposition with a proper regularization. For a core-collapse supernova at a distance of $10~{\rm kpc}$ in the Milky Way, its $\overline{\nu}^{}_e$ spectrum can be precisely determined from the inverse beta-decay process $\overline{\nu}^{}_e + p \to e^+ + n$, for which a $20~{\rm kiloton}$ liquid-scintillator detector with the resolution similar to the Jiangmen Underground Neutrino Observatory (JUNO) may register more than 5000 events. We have to rely predominantly on the elastic neutrino-electron scattering $\nu + e^- \to \nu + e^-$ and the elastic neutrino-proton scattering $\nu + p \to \nu + p$ for the spectra of $\nu^{}_e$ and $\nu^{}_x$, where $\nu$ denotes collectively neutrinos and antineutrinos of all three flavors and $\nu^{}_x$ for $\nu^{}_\mu$ and $\nu^{}_\tau$ as well as their antiparticles. To demonstrate the validity of our approach, we also attempt to reconstruct the neutrino spectra by using the time-integrated neutrino data from the latest numerical simulations of delayed neutrino-driven supernova explosions.


I. INTRODUCTION
The precious observations of the neutrino burst from Supernova (SN) 1987A in the Large Magellanic Cloud by Kamiokande-II [1], Irvine-Michigan-Brookhaven [2] and Baksan [3] have essentially demonstrated that the delayed neutrino-driven explosion mechanism works well for a core-collapse supernova [4][5][6][7][8]. Their experimental data are consistent with the inverse beta-decay ν e + p → e + + n (IBD) initiated by SN ν e of an average energy E ν e ≈ 12 MeV, and the total gravitational binding energy released mainly via neutrinos is about 3 × 10 53 erg [9,10] for a duration of 10 seconds. However, it is impossible to get more useful and robust information on SN neutrino energy spectra from two dozens of neutrino events totally observed in all those three experiments. See, e.g., Ref. [11], for an excellent overview on the statistical analysis of neutrino data from SN 1987A.
Thanks to recent tremendous progress in neutrino oscillation experiments, a few large water-Cherenkov (WC), liquid-scintillator (LS) and liquid-argon time projection chamber (LAr-TPC) neutrino detectors have already been in operation (e.g., Super-Kamiokande [12], Borexino [13] and KamLAND [14]) or will be available in the near future (e.g., Hyper-Kamiokande [15], JUNO [16] and DUNE [17]). In the WC and LS detectors, the IBD channel is dominant and most sensitive to SN ν e , while the charged-current reaction ν e + 40 Ar → e − + 40 K * in the LAr-TPC detectors is crucially important for probing the early-time neutronization burst of ν e [18]. As was first pointed out in Ref. [19], the elastic scattering of SN neutrinos of all flavors on the protons in the LS detectors can hopefully be implemented to extract the spectral information of ν µ and ν τ and their antiparticles (which will be collectively denoted as ν x ). Recently, this idea has been further studied in Ref. [20], where it has been shown that the energy spectrum of SN ν x can be partially reconstructed from the elastic neutrino-proton scattering (pES) events. Therefore, for a future galactic core-collapse SN, we will be able to achieve a high-statistics measurement of neutrino signals in multiple complementary channels.
In this paper, we study whether it is possible to accomplish a complete reconstruction of the energy spectra of SN neutrinos of all three distinct flavors, i.e., ν e , ν e and ν x , from a single large LS detector, such as JUNO [16] under construction and the proposed LENA [21].
For the JUNO detector, which is designed to have 20 kiloton LS with a proton fraction about 12%, there will be about 5000 IBD events for a galactic SN at a distance of 10 kpc. At the same time, about 1500 pES events will be registered above the visible-energy threshold of 0.2 MeV, while the number of elastic neutrino-electron (eES) events is approximately 400 [16,22]. For the LENA detector with the similar proton fraction and the energy resolution of 7% at 1 MeV [21], the corresponding statistics will be increased by a factor of about 2.5 due to the total target mass scaled upward by the same factor [21]. Apart from the target mass, the statistics will also be increased (or decreased) for a closer (or farther) SN. The motivation of a full reconstruction of SN neutrino energy spectra is two-fold. First, the formation of energy spectra of SN neutrinos is obviously associated with the production processes [23][24][25], so a direct extraction of neutrino energy spectra from future observations will be helpful in understanding the production of SN neutrinos and exploring the true explosion mechanism [26]. Second, besides the standard Mikheyev-Smirnov-Wolfenstein (MSW) matter effects [27,28], significant flavor conversions of SN neutrinos induced by the coherent neutrino-neutrino scattering are likely to take place in the SN environment where a dense neutrino gas is present [29][30][31][32][33][34]. See, e.g., Refs. [35][36][37], for recent reviews on collective neutrino oscillations and latest developments. To verify the spectral splits and swaps of SN neutrinos caused by their self-interactions, one has to determine SN neutrino spectra of individual flavors.
The present work differs from Ref. [20] at least in two aspects. First, the practically powerful approach of singular value decomposition (SVD) with a proper regularization is applied to unfold the energy spectra of ν e , ν e and ν x mainly from the eES, IBD and pES channels, respectively. The energy resolution of the future large LS detector is fully taken into account, while an approximate and analytical method is adopted only for ν x in Ref. [20].
In order to demonstrate the validity of the SVD approach, we also attempt to extract the SN neutrino spectra from the neutrino data produced by the numerical simulations of SN explosions. Second, for a single LS detector, we show that it is possible to accomplish a complete reconstruction of ν e , ν e and ν x spectra from the observations of a future galactic SN.
Certainly, the combination of neutrino data from different types of detectors (e.g., Hyper-Kamiokande, JUNO and DUNE) running at the same time will further improve greatly the determination of SN neutrino spectra.
The remaining part of our paper is structured as follows. In Sec. II, we explain how the neutrino events in different channels are calculated for a typical JUNO-like detector. Then, the strategy to reconstruct SN neutrino spectra is outlined in Sec. III, where both the basic idea in the ideal case and the realistic SVD unfolding approach are briefly introduced.
We show the numerical results of reconstruction in both cases of the parametrized neutrino spectra and the simulated ones in Sec. IV. Finally, we summarize in Sec. V our main conclusions.

II. SUPERNOVA NEUTRINO EVENTS
The relevant reactions for SN neutrinos in the LS detectors have been studied in details in Refs. [16,22,38]. Although the charged-and neutral-current interactions of SN neutrinos with the 12 C and 13 C nuclei are also available there, the corresponding numbers of neutrino events will be only sub-dominant [22]. For simplicity, we just focus on the dominant channels, i.e., IBD for ν e , pES for ν x and eES for ν e , respectively.
In this work we take a 20 kiloton LS detector with the resolution similar to JUNO as the working example to explore the reconstruction performance. The detector energy resolution is assumed to be 3%/ E o /(MeV) with E o being the event observed energy, and the energy threshold is simply taken as 0.2 MeV for both the pES and eES interaction channels. The calculations of neutrino events in these three channels are quite straightforward and will be briefly summarized below.

A. SN Neutrino Spectra
The mechanisms of SN neutrino production are distinct at three different stages of SN evolution, namely, the early-time neutronization burst, the accretion phase and the cooling phase. However, the differential neutrino fluences or time-integrated neutrino energy spectra can be perfectly described by the Keil-Raffelt-Janka (KRJ) parametrization [23] where the subscript α runs over three neutrino flavors ν e , ν e and ν x , and the flavor-dependent total neutrino energy ε α is given in units of 5 × 10 52 erg. In addition, the distance D to a galactic SN is normalized to a typical value of 10 kpc, the neutrino energy E α and average energy E α are measured in MeV, and the spectral index γ α = 3 will always be adopted in our calculations. Note that γ α = 2 corresponds to the Maxwell-Boltzmann distribution, while γ α = 2.3 to the Fermi-Dirac distribution with a zero chemical potential [39]. In the assumption of energy equipartition, the total gravitational binding energy E g = 3 × 10 53 erg released in a core-collapse SN within ten seconds is shared by ν e , ν µ , ν τ and their antiparticles, i.e., ε α ≈ 5 × 10 52 erg. In our calculations, the nominal values of neutrino average energies will be set to E ν e = 12 MeV, E ν e = 14 MeV and E ν x = 16 MeV.
For the numerical simulations, the neutrino data at a given time provided by the Garching group [40] are fitted by the KRJ parametrization, and the total energy ε α , the average energy E α and the spectral index γ α are tabulated, so one has to integrate the time-dependent neutrino number fluxes over the emission time, instead of using Eq. (1). In contrast, for the numerical simulation models from the Japan group [41], the neutrino data are provided for both time and energy distributions. In this case, we integrate over time to obtain the neutrino energy spectra for each simulation model.

IBD Events
Due to its large cross section and clear signal, the IBD is the primary reaction for SN neutrino detection in the LS detectors. For an average energy of E ν e = 14 MeV and a galactic SN at D = 10 kpc, a JUNO-like detector will register about 5000 IBD events, which are adequate for a precise determination of the SN ν e spectrum with one percent errors [22] for both the total energy ε ν e and the average energy E ν e . The time coincidence of prompt and delayed signals from e + and n in the final states, respectively, renders the IBD to be almost free of background. The IBD cross section can be calculated accurately, in particular for MeV neutrinos [42,43]. For illustration, we take the following simple formula where the positron energy is given by E e + ≈ E ν e − ∆ np with the neutron-proton mass difference ∆ np ≡ m n − m p ≈ 1.293 MeV, and the positron momentum is p e + = E 2 e + − m 2 e . The energy threshold for incident neutrinos is E th ν e = ∆ np + m e ≈ 1.804 MeV, while the visible energy in the detector is E v = E e + + m e after the electron-positron annihilation. Therefore, the event spectrum can be figured out as where E o stands for the observed energy, and G(E o ; E v , δ E ) for the Gaussian function of E o with E v and δ E = 3%/ E o /(MeV) being the expectation value and the standard deviation, respectively. Given the target mass of 20 kiloton and a proton fraction of 12% for the typical LS target, one find that the total number of protons is N p ≈ 1.4 × 10 33 .
The integration of the event spectrum in Eq. (3) over the observed energy E o gives rise to the total number of IBD events. Assuming that the gravitational binding energy takes the typical value of E g = 3 × 10 53 erg for a core-collapse supernova and it is equally shared by neutrinos and antineutrinos of three flavors, which are emitted within a duration of ten seconds, one can obtain the number of IBD events around 5000 for E ν e = 14 MeV and γ ν e = 3 [16]. In addition, the excellent energy resolution δ E = 3%/ E o /(MeV) and simply linear relationship between the neutrino energy and the visible energy, namely, E v ≈ E ν e − 0.782 MeV, assure an accurate enough determination of SN ν e spectrum.

pES Events
As an advantage of LS detectors, the pES events can also be observed due to their low threshold of visible energies, if radioactive backgrounds and dark noises of photomultiplier tubes are well in control [16]. The general formula for the cross section of the pES can be found in Ref. [44], and can be recast into the following form, which is particularly useful for supernova neutrinos [19] dσ where only the zeroth-order term in E ν /m p is retained, T p is the recoil energy of the finalstate proton. The recoil energy will be significantly quenched in the LS, according to the Birks' law [45], such that only part of the kinetic energy yields the scintillation light, i.e., where the Birks' constant k B ≈ 0.01 cm MeV −1 and dE/dx denotes the corresponding energy-loss rate of protons in the LS materials. Although the observation of the pES events depends crucially on the radioactive backgrounds, such as the beta decays of the nuclear isotope 14 C in the LS, and the dark noises, it is reasonable to assume that an energy threshold of E th o = 0.2 MeV for the single events is achievable in future LS detectors. In this case, the pES event spectrum can be calculated as which will be further convoluted with the energy resolution function G(E o ; T ′ p , δ E ) to give the event spectrum with respect to the observed energy E o .
Note that according to the kinematics of the pES, to produce a proton with a recoil energy T p , we require the neutrino energy to be above E min α = (T p m p /2) 1/2 . In consideration of the backgrounds, only the events with E o > E th o can be identified. The number of pES events is significant due to the contributions of both ν µ and ν τ and their antiparticles, whose average energy is higher than those of ν e and ν e . It is easy to find that the total number of pES events is about 1500 for E th o = 0.2 MeV, while about 2000 for E th o = 0.1 MeV. This observation indicates that the reduction of radioactive backgrounds and dark noises is crucially important for the reaction channels with singles.

eES Events
The eES is another very important channel to detect SN neutrinos in the sense that it is sensitive to the early-time ν e burst from neutronization, which has been found to be a robust feature of SN neutrinos from numerical simulations [26]. The cross sections of neutrinos scattering off electrons at the tree level are well known in the standard theory of electroweak interactions [46,47] where the Fermi constant is G F = 1.166 × 10 −5 GeV −2 , and the kinetic energy of the final- For electron neutrinos, the coefficients are given by ǫ − = −1/2 − sin 2 θ w and ǫ + = − sin 2 θ w ; while ǫ − = 1/2 − sin 2 θ w and ǫ + = − sin 2 θ w for muon and tau neutrinos. For the cross sections of antineutrinos, one has to perform the exchange of ǫ − ↔ ǫ + in Eq. (7). As the recoil energy of electrons will entirely be converted into the scintillation light, the observed event spectrum can be derived as follows where N e is the total number of electrons in the detector and the minimal neutrino energy E min α ≈ T e /2 + T e (T e + 2m e )/2 required for an electron with a recoil energy of T e can be figured out from the kinematics. The contributions from neutrinos and antineutrinos of all three flavors are included, although the dominant one comes from electron neutrinos.
Similar to the pES events, the eES events with energies below E th o = 0.2 MeV will be hidden in the radioactive backgrounds and dark noises.
For illustration, the IBD and pES event spectra at a JUNO-like LS detector for a galactic SN at a distance of 10 kpc have been shown in the first row of Fig. 1. In the second row, the eES events are displayed as well, where the number of events is given in both linear and logarithmic scales. In all these plots, the theoretical expectations of event numbers are denoted by black solid curves, while the toy Monte Carlo (MC) results by blue dots, together with the bin widths and statistical uncertainties being represented respectively by horizontal and vertical error bars. Note that the toy MC samples are randomly generated within the ROOT framework [48] according to the event spectra in Eqs. (3), (6) and (8), where an energy resolution of 3%/ E o /(MeV) is assumed.

III. STRATEGY FOR RECONSTRUCTION
As demonstrated in Ref. [22], the future large LS detectors could separately measure the average energies of SN ν e , ν e and ν x with reasonably good precisions, if the spectral indices are fixed at γ α = 3 (for α = ν e , ν e , ν x ). However, without any presumed information, is it possible to reconstruct neutrino spectra of all three flavors directly from experimental data?
In this work, we show that the answer is positive, and attempt to unfold all the SN neutrino spectra directly from simulated experimental neutrino data.

A. The Ideal Case
Ignoring the energy smearing, it is quite straightforward to obtain the SN ν e spectrum from the IBD events. In practice, we divide the whole energy range of the observed IBD events into N bins with the bin width , and the number of events in each bin is denoted by n i . The event distribution can be repre- According to the IBD event rate given in Eq. (3) and the corresponding relationship MeV between the neutrino energy E i and the observed energy E ′ i , we have where . Therefore, the initial ν e spectrum is reconstructed as f i at the energy E i , where the same pattern of binning for the initial neutrino energy is assumed.
Next, we proceed with the reconstruction of the SN ν x spectrum from the pES events.
The distribution of pES events can be denoted as n i in the energy bin of [ Since the observed recoil energy T ′ i of the final-state proton corresponds to T i due to the quenching effects, one has to implement the Birks' law in Eq. (5) to establish the relationship between T ′ i and T i , or equivalently to convert the observed event spectrum into that with respect to T i . To be explicit, we have where being the inverse function of that defined in Eq. (5) and is the bin width of the initial recoil energy. Furthermore, in order to produce the final-state proton with a recoil energy of T i , we need an initial-state neutrino of energies above the minimum E i min = (T i m p /2) 1/2 , and all the neutrinos with E ≥ E i min can give rise to the pES events in the bin at T i .
Unlike in the IBD case, we shall now categorize the initial neutrino energy into N different is fixed by the minimal neutrino energy to produce a recoil energy T i−1 for i = 1, 2, · · · , N. Since the average energy of ν x is larger than those of ν e and ν e , the pES receives dominant contributions from ν x . To illustrate the reconstruction method, we temporarily neglect the sub-leading contributions from ν e and ν e . All the contributions will be taken into account for the realistic case in the next section. Defining f j ≡ dF ν x /dE ν x | E j as the fluence at E j ≡ (E j + E j−1 )/2, one can establish the relationship where the bin width ∆E j ≡ E j − E j−1 and the matrix element K ij can be identified without any ambiguity. Note that the summation is performed over the bin index j corresponding to E j ≥ E i min = (T i m p /2) 1/2 , so the matrix element K ij vanishes for i > j. To be explicit, we have  One can choose the energy bins of the observed energy T ′ in such a way that the numbers of events in any bins are comparable in magnitude. Furthermore, the upper-triangular form of the matrix K ij leads to an analytical inversion where it is unnecessary to calculate K −1 directly. The reason is simply that one can first obtain f N = (n N /∆T ′ N )/K N N and then fix f i iteratively for i = N − 1, N − 2, · · · , 2, 1. Finally, let us consider the eES channel, for which the recoil energy of the final-state electron can be accurately measured and we assume that the possible backgrounds of other singles will be reduced by applying the approach of pulse shape discrimination. As before, the observed recoil energies of electrons can be divided into N bins [T 0 , T 1 ], [T 1 , T 2 ], · · · , [T N −1 , T N ], where the recoil energies are identified as the observed energies without any quenching effects. In the eES channel, since the cross section of ν e is much larger than those of ν e and ν x , we take ν e as an example to show the reconstruction method. Under this approximation, the eES event spectrum can be readily calculated as where n i stands for the event number in the observed energy bin [T i−1 , T i ] with the width In addition, the minimal neutrino energy E i min required to produce the recoil energy T i is given by E i min = T i /2 + T i (T i + 2m e )/2, indicating that E i min ≈ T i for a relatively large recoil energy T i ≫ 1 MeV. The initial neutrino energy bins [E j−1 , E j ] for j = 1, 2, · · · , N can be chosen as E j−1 = E j−1 min such that all the energy bins with indices j ≥ i contribute to the event bin i. Therefore, one can immediately realize that the approach for the pES channel can also be applied to the eES channel. The number of eES events, however, will be much smaller, leading to a worse precision on the extraction of the ν e spectrum.

B. The Realistic Case
In the ideal case, we have neglected the finite energy resolution and acceptance limitation of the detector. To incorporate them, we shall employ the unfolding technique to reconstruct the SN neutrino spectra. This technique is quite generic and frequently encountered in particle physics [49,50] to extract the true distribution from experimental data. Here we follow the strategy proposed in Ref. [51], where the SVD approach is clearly addressed and the necessity of an additional regularization is well explained. For completeness, we first briefly introduce the SVD method, and then apply it to the unfolding of SN neutrino spectra in the next section.
The spectral unfolding can be identified as a linear inverse problem ofÂx = b, whereÂ is an m × n response matrix characterizing the detector response effects (e.g., the energy resolution and the limited acceptance). Additionally, x denotes an n-dimensional vector of the true distribution, while b is an m-dimensional vector of the observed distribution. For the SN neutrino spectrum under consideration, x is just the predicted neutrino event spectrum, which is roughly the convolution of the SN neutrino spectrum with the cross section, and b is the observed event distribution. In fact,Â can be calculated from a large sample of MC simulations, and for each bin of the true energy, the summation of each column ofÂ (corresponding to the whole range of observed energies) is normalized to one.
To extract the true distribution, one can decompose the response matrix asÂ = USV T , where U and V are respectively the m × m and n × n orthogonal matrices while S is an m × n diagonal matrix with only non-negative and non-increasing diagonal elements. Now the unfolding problem can be solved by using the generalized inverse matrix ofÂ. However, due to some small singular values in S and the measurement errors in b, a direct extraction of x with the exact inversion of the response matrix may lead to an unacceptable oscillatory solution. The standard way to suppress the spurious oscillating components is to introduce a regulator, which is actually an extra penalty term in the least-squares function of the corresponding inverse problem with a regularization parameter characterizing the weight of the penalty. The optimal value of the regularization parameter can be obtained by balancing After preparing the detector response matrices for the IBD, pES and eES channels, we can apply the SVD unfolding method to the reconstruction of the SN neutrino spectra. In the unfolding processes, the binning schemes for the observed spectra depend on the event statistics and are carefully treated in order to have comparable numbers of neutrino events in each observed energy bin. For the true energy spectra, we employ the equal-size binning schemes but combine the bins near the boundaries due to the limited statistics. The actual binning scheme for each unfolding process can be readily read out from the energy spectra which will be presented in the next section. The practical realization of the SVD unfolding algorithm is based on the RooUnfold package [52] and the choice for the optimal value of the regularization parameter can be found in Ref. [51].

IV. FINAL RECONSTRUCTION RESULTS
In this section, we present the final results on the reconstruction of SN neutrino energy spectra. First of all, we simulate a large number of trials for SN neutrino data in the IBD, pES and eES channels by assuming the SN distance of 10 kpc and the detector parameters discussed in Sec. III. As we concentrate on the power of reconstruction, flavor conversions of SN neutrinos will temporarily be left aside. For the observed spectra in each trial, the SVD unfolding algorithm is applied to obtain the reconstructed neutrino spectra. To illustrate the dependence of the unfolding performance on the event statistics, we also simulate the neutrino data for the SN at the distances of 1.0 kpc and 0.2 kpc and reconstruct neutrino spectra in these two cases. To explore the impact of regularization, we also compare the reconstruction results of neutrino spectra in the ideal case and those with the SVD unfolding algorithm. Moreover, we use the method of bin-to-bin separation to extract the original neutrino spectra of three flavors in their common energy range. Finally, we reconstruct the SN neutrino spectra by using the numerical SN models from both the Garching group [40] and the Japan group [41].

A. Reconstructed Neutrino Energy Spectra
The reconstructed neutrino energy spectra for the IBD, pES and eES channels, which are actually the SN neutrino fluences weighted with the corresponding cross sections, are shown in Fig. 3. The simulated data are for a SN at 10 kpc, and the KRJ parametrization of the SN neutrino fluences with E ν e = 12 MeV, E ν e = 14 MeV and E ν x = 16 MeV has been adopted. The true SN neutrino energy spectra are shown as black histograms. The vertical error bars are the statistical uncertainties arising from the those of the observed SN neutrino spectra, while the horizontal ones show the bin widths. From the plots one can observe that the reconstructed neutrino energy spectra are nicely in accordance with the true ones and their differences are within the statistical fluctuations which indicates that the biases from the regularization are significantly smaller than statistical uncertainties in the corresponding energy bins.
To explore the impact of regularization, both the method in the ideal case and the SVD method are applied to the same observed spectra of the IBD and pES channels as given in Fig. 3. The reconstructed neutrino spectra via these two methods, where the results of the ideal case are denoted by black square points and those using the SVD method by blue circle points, are shown in Fig. 4. The error bars represent the statistical uncertainties propagated from the observed spectra. For the ideal case, the reconstructed spectra are calculated according to Eqs. (9) and (11) for the IBD and pES channels, respectively. From such a comparison, one can easily observe the advantages of the SVD unfolding method, namely, more stable reconstructed results and lower statistical errors.
In the left column of Fig. 5, the results of reconstructed neutrino spectra for the SN distances of 10 kpc, 1 kpc and 0.2 kpc are summarized. In each panel, the reconstructed spectra after unfolding for the IBD, pES and eES channels are represented by black, blue and red points, respectively. The histograms are the corresponding true neutrino energy spectra. We can see that the reconstructed and true energy spectra well agree with each other within the statistical uncertainty. For the SN distance at 10 kpc, the best precision of the IBD, pES and eES channels can reach the levels of 2%, 10% and 50%, respectively.
When the SN distances become closer and thus the statistics becomes larger, the precisions get better, indicating that the systematic bias of the unfolding method itself is negligible.

B. Original Neutrino Energy Spectra
The reconstructed neutrino energy spectra from the IBD, pES and eES channels are essentially the original neutrino fluences weighted by cross sections. To further extract the original spectra within their common energy region (e.g., E ν 20 MeV), we choose the same  Fig. 3 have been adopted. The results of the ideal case are denoted by black square points and those using the SVD method by blue circle points. The error bars are statistical uncertainties propagated from the observed spectra.
binning scheme for the reconstructed neutrino spectra in all three channels and implement the bin-to-bin separation method, which can be clearly interpreted as the following equations where n IBD (E i ν ), n pES (E i ν ) and n eES (E i ν ) stand for the reconstructed distributions in the neutrino energy bin E i ν for the IBD, pES and eES channels, and F α (for α = ν e , ν e , ν x ) are the original SN neutrino fluences. Here N p and N e are the numbers of free protons and electrons in the LS target, while σ IBD , σ αp , σ αe are the cross sections of the IBD, pES and eES reactions, respectively. By solving the above three equations for each neutrino energy bin E i ν , we achieve a complete reconstruction of the SN ν e , ν e and ν x energy spectra, as shown in the right column of Fig. 5 for the SN distances of 10 kpc, 1 kpc and 0.2 kpc.
From the upper-right plot of Fig. 5, where the reconstruction has been performed in the case of a SN distance of 10 kpc, it is clear that theν e energy spectra from the IBD channel can be determined with a precision better than 2% in the vicinity of E ν = 20 MeV. The high-energy tail of the ν x spectrum starting from 20 MeV are also well reconstructed, where the statistical uncertainty increases from 10% at E ν = 20 MeV to 30% at E ν = 30 MeV.
However, for just one trial of the SN ν e events at the distance of 10 kpc, it is difficult to provide an accurate reconstruction of the ν e spectrum because of the large statistical fluctuation, which can be at the level of 100%. On the other hand, with a much larger statistics when the SN distance is either 1 kpc or 0.2 kpc, one can observe an excellent reconstruction performance for all three neutrino flavors. This demonstrates the powerful capability of the SVD unfolding method. It is worthwhile to point out that SN ν e can also be observed in the LS detectors via the charged-current ν e -12 C interaction [22], thus the incorporation of the ν e -12 C channel will improve the precision of reconstruction.

C. Reconstruction with Numerical SN Models
In the previous discussions, the KRJ parametrization of SN neutrino fluences has been used for numerical simulations. In order to demonstrate the robustness of the SVD unfolding approach, we now consider the SN neutrino fluences directly from the numerical SN models, which have been provided by the Garching [40] and the Japan [41] groups. In total, 16 numerical models from the Garching group and 21 models from the Japan group are considered, covering a variety of progenitor star masses and different treatment of neutrino transport. For each numerical SN model we simulate 500 trials of the neutrino data and then apply the same reconstruction algorithm. For each trial, the bias is defined as the relative difference between the reconstructed and true neutrino spectra in each energy bin, which includes both the statistical fluctuation and the actual bias from the unfolding method. To reduce the statistical fluctuation, we further define the averaged bias among 500 trials as: where k is the index for the trial and N = 500 is the number of total trials for one numerical model. This definition is applicable to all the IBD, pES and eES channels.
In Fig. 6, the averaged bias of different numerical SN models for a SN distance of 10 kpc has been given for the reconstruction of SN ν e , ν x and ν e from the IBD, pES and eES channels, respectively. The left column is for 16 Garching numerical SN models, while the right column for the 21 Japan numerical SN models. From these plots, we can observe Comparing the left and right columns, one can see larger variations for the numerical SN models from the Garching group than those from the Japan group. This might be due to the fact that the Garching models cover an even wider range of core-collapse SNe, for which the basic properties and the main features of neutrino fluences differ greatly. Finally we conclude that the SVD unfolding method is a useful and robust method to reconstruct the neutrino energy spectra both in the IBD, pES and eES channels and for the three flavor neutrinos, at least for a particular region of neutrino energies.
Before closing this section, we make some remarks on the flavor conversions of SN neutrinos, which have so far been ignored in our discussions. Both the ordinary MSW matter effects [27,28] and the collective oscillations [29][30][31][32][33][34][35][36][37] lead to significant flavor conversions when neutrinos are propagating from the neutrino sphere to the surface of the SN. If only the MSW effects are taken into account, which become relevant in the SN mantle, there will be a partial or complete conversion between ν e and ν x (and between ν e and ν x ), depending on the neutrino mass ordering. Therefore, the final energy spectrum of a given neutrino flavor would be the mixture of the initial spectra of two neutrino flavors with distinct average energies. On the other hand, the collective neutrino oscillations, which are induced by the coherent forward neutrino-neutrino scattering, may result in spectral splits of SN neutrino spectra [36]. However, it remains an open question whether the collective neutrino oscillations do happen in the real SN environment [37].
The SVD unfolding method under discussions is independent of any specific parametrization of the neutrino fluxes as well as the spectra at the detector, no matter whether neutrinos flavor conversions are included. As we have shown before, the event statistics is one of the most important factors that affect the reconstruction performance. The reconstruction tests with different numerical simulation models from the Garching and Japan groups demonstrate the robustness of the unfolding method for different inputs of the neutrino fluxes and energy spectra. On the other hand, the reconstructed neutrino spectra of three flavor neutrinos at the detector would be very useful to test the paradigm of neutrino flavor conversions, which will be presented in a separated work in the future.
For a future galactic core-collapse SN, a high-statistics measurement of neutrino signals can definitely be achieved. In particular, SN neutrinos of all three flavors ν e , ν e and ν x can be detected in the IBD, eES and pES channels in a single large LS detector, such as JUNO. The precise determinations of SN neutrino spectra and flavor contents will be crucially important to verify the delayed neutrino-driven explosion mechanism, and to reveal the complicated pattern of SN neutrino flavor conversions [26].
In this paper, we apply the SVD approach with a proper regularization [51] to reconstruct SN neutrino energy spectra for all flavors in a 20 kiloton LS detector with the resolution similar to JUNO. First of all, we explain how to calculate the SN neutrino events, and introduce an analytical recipe to reconstruct neutrino spectra in the ideal case where the finite energy resolution is neglected. Then, taking account of the realistic energy resolution, we simulate the neutrino data by using the parametrized neutrino spectrum (e.g., the KRJ parametrization with E ν e = 12 MeV, E ν e = 14 MeV and and E ν x = 16 MeV). Starting with the reconstructed spectra, we also combine the IBD, eES and pES channels to extract the original neutrino spectra in a common region of neutrino energies. It turns out that for a SN at 10 kpc there is no problem to accurately reconstruct ν e and ν x spectra, but the reconstruction of ν e is limited by the relatively low statistics. Finally, instead of the parametrized SN fluences, the numerical SN models from the Garching and Japan groups are analyzed as well. In addition, the robustness and validity of the SVD approach for both analytical and numerical neutrino data are demonstrated to be excellent.
Even though a single large LS detector like JUNO has the great potential to detect SN neutrinos of all three flavors, the statistics in the eES channel is quite low for a SN at 10 kpc. In this case, it is very important to notice that the WC detectors (e.g., Hyper-Kamiokande) and the LArTPC detectors (e.g., DUNE) are exceptionally good at observing the SN ν e . Therefore, different types of future large detectors are actually complementary to each other. A combined analysis of the LS, WC and LArTPC detectors in reconstructing SN neutrino spectra will be very interesting and deserve dedicated studies in the near future.