A model-independent approach to the reconstruction of multi-flavor supernova neutrino energy spectra

The model-independent reconstruction of the energy spectra of $\overline{\nu}^{}_e$, $\nu^{}_e$ and $\nu^{}_x$ (i.e., $\nu^{}_\mu$, $\nu^{}_\tau$ and their antiparticles) from the future observation of a galactic core-collapse supernova (SN) is of crucial importance to understand the microscopic physics of SN explosions. To this end, we propose a practically useful method to combine the multi-channel detection of SN neutrinos in a large liquid-scintillator detector (e.g., JUNO), namely, the inverse beta decay $\overline{\nu}^{}_e + p \to e^+ + n$, the elastic neutrino-proton scattering $\nu + p \to \nu + p$ and the elastic neutrino-electron scattering $\nu + e^- \to \nu + e^-$, and reconstruct the energy spectra of $\overline{\nu}^{}_e$, $\nu^{}_e$ and $\nu^{}_x$ by making the best use of the observational data in those three channels. In addition, the neutrino energy spectra from the numerical simulations of the delayed neutrino-driven SN explosions are implemented to demonstrate the robustness of our method. Taking the ordinary matter effects into account, we also show how to extract the initial neutrino energy spectra in the presence of neutrino flavor conversions.

To diagnose the true pattern of neutrino flavor conversions and even pin down the SN explosion mechanism [23], we have to detect SN neutrinos of all flavors and fully reconstruct their energy spectra. However, the sparse data from SN 1987A do not suffice for such a purpose [24]. For a galactic core-collapse SN at a distance around 10 kpc [25], a few large water-Cherenkov (WC), liquid-scintillator (LS) and liquid-argon time projection chamber (LAr-TPC) neutrino detectors worldwide that are in operation (e.g., Super-Kamiokande [26], Borexino [27] and KamLAND [28] ) or under construction (e.g., Hyper-Kamiokande [29], JUNO [30], and DUNE [31]) will register a large number of neutrino events and have a great potential to provide complete flavor information on SN neutrinos. First of all, the energy spectrum of ν e can be unambiguously determined via the inverse beta decay channel ν e + p → e + + n (IBD) in both WC and LS detectors. Then, the ν e energy spectrum can be well measured in the LAr-TPC through the charge-current interaction ν e + 40 Ar → e − + 40 K * , and in the WC detector via the elastic neutrino-electron scattering ν + e − → ν + e − (eES), which receives the contributions from all neutrino flavors but is most sensitive to ν e because of the largest cross section. Finally, the energy spectrum of ν x , which collectively denotes ν µ , ν τ and their antiparticles, can be partially extracted from the elastic neutrino-proton scattering channel ν + p → ν + p (pES) in the LS detectors, as first proposed in Ref. [32] and further studied in Ref. [33].
In the previous work [34], we have already shown that it is possible to accomplish a complete reconstruction of the energy spectra of SN neutrinos ν e , ν e and ν x in a single large LS detector. One salient feature of the LS detector is its low energy threshold, which renders it capable of observing the recoiled protons resulting from the pES process that is mainly sensitive to SN ν x . Therefore, the energy spectra of ν e , ν e and ν x can essentially be extracted from the IBD, eES and pES events, respectively. In this paper, we have improved the previous study in Ref. [34] in several important aspects. First, while the SN ν e energy spectrum can be precisely determined via the IBD data, the energy spectrum of ν e cannot be accurately extracted from the eES data. The key point is that although the eES cross section of ν e is about six times larger than that of ν µ or ν τ (or their antiparticles), the total contribution to the eES events from the latter four neutrino flavors is obviously significant.
In Ref. [34], it has been assumed that ν e dominates over all other flavors in the eES channel and the detector response matrix between the initial neutrino energy and the observed event energy is approximately taken to be universal for all neutrinos. In the present work, this assumption is relaxed and the full detector response matrix will be implemented. Second, the flavor conversions of SN neutrinos, although only the ordinary MSW matter effects in the SN mantle are taken into account, serve as another complication for the spectrum unfolding of multi-flavor neutrinos. Third, the central idea of our reconstruction method is to treat ν e , ν e and ν x energy spectra on the same footing in all three reaction channels, and build the overall detector response matrix according to their individual interactions with the target particles in the LS. We stress that such a model-independent approach is also applicable to solar neutrinos (with two flavors ν e and ν µ/τ ) and ultrahigh-energy cosmic neutrinos (with three flavors ν e /ν e , ν µ /ν µ and ν τ /ν τ ), when the statistics is sufficiently large in the relevant next-generation experiments.
The remaining part of our paper is organized as follows. After a brief description of SN neutrino detection in the LS detector in Sec. II, we introduce our reconstruction method in Sec. III. Then, we present in Sec. IV the unfolding results of SN neutrino energy spectra, and investigate the impact of the energy threshold of the detector and the dependence on the numerical models of SN neutrinos. Moreover, the flavor conversions of SN neutrinos in the presence of MSW matter effects are discussed. Finally, we summarize our main results and conclude in Sec. V.

II. SUPERNOVA NEUTRINO EVENTS
In a core-collapse SN, the gravitational potential energy about 3×10 53 erg is released and 99% of it is carried away by neutrinos. Neutrinos of all flavors with energies of a few tens MeV are emitted. The duration of neutrino emission lasts for about ten seconds in three distinct phases, namely, the early-time neutronization burst, the accretion phase and the cooling phase. In this work, however, we focus on the time-integrated SN neutrino energy spectra and their reconstruction from simulated experimental data.
In LS detectors, the relevant reactions for SN neutrinos have been studied in detail in Refs. [30,35,36]. Although the charged-and neutral-current interactions with the 12 C and 13 C nuclei are also available therein, the corresponding numbers of neutrino events are subdominant [35]. For simplicity, we consider only three dominant channels, i.e., the IBD, eES, and pES. The neutrino event rates in these channels for a JUNO-like detector have been calculated in Ref. [34] and will be recapped below for completeness.

A. SN Neutrino Spectra
The differential fluences or time-integrated energy spectra of SN neutrinos can reasonably be described by the Keil-Raffelt-Janka (KRJ) parametrization [9] dF α dE α = 3.5 × 10 13 where the subscript α runs over three neutrino flavors ν e , ν e and ν x , and ε α is the total neutrino energy for each flavor 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. As an example for the analytical model of SN neutrino energy spectra, the spectral index γ α = 3 will always be adopted. In the assumption of energy equipartition, the total gravitational binding energy E g = 3 × 10 53 erg released in the SN explosion within ten seconds is shared by ν e , ν µ , ν τ and their antiparticles, i.e., ε α ≈ 5 × 10 52 erg. Whenever the analytical model is referred to, 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 numerical models of SN neutrino spectra, we make use of the simulation results from the Japan group [37], where the time and energy distributions of all neutrino flavors are provided for a wide range of progenitor star masses and different values of metallicity. In this case, we numerically integrate the two-dimensional distribution over time to obtain the neutrino energy spectra for each simulation model.
It is worthwhile to mention that all the neutrino energy spectra from either the analytical model with the KRJ parametrization or the numerical models from the simulations of SN explosions are implemented as the original ones in the following discussions. Only in Sec. IV C do we investigate the impact of neutrino flavor conversions on the reconstruction strategy.

B. SN Neutrino Events
In Ref [34], The IBD channel is the golden channel for SN neutrino detection in LS detectors due to the time coincidence of the prompt and slow signals and the large cross section of the IBD reaction [38,39]. This channel is solely sensitive to ν e . The other two elastic scattering channels, pES and eES, receive contributions from neutrinos of all flavors ν e , ν e and ν x . For SN neutrinos with a few tens of MeV, the cross section for pES is almost the same for all flavors [32,40]. In the eES channel, the cross sections of electron flavor neutrinos ν e are larger than that of ν x [41,42] due to the charged-current interaction of the former. The signals of recoiled protons and electrons in the LS detectors, from the pES and eES channel respectively, can be distinguished by utilizing the technique of pulse shape discrimination.
In the following discussions, all the events in those three channels are treated ideally without considering the uncertainties from background signals and detection efficiencies. The effects   of 12 C-related interaction channels on the SN neutrino spectra reconstruction will be studied in a future separated work.
For illustration, the IBD and pES event spectra in 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 event spectrum is displayed as well, where the number of events has been given in both linear (left panel) and logarithmic (right panel) scales. In all these plots, the theoretical predictions for the event numbers are denoted by black solid curves, while the toy MC results by blue dots, for which the widths of energy bins and the statistical uncertainties are respectively represented by horizontal and vertical error bars. Note that the toy MC samples are randomly generated within the ROOT framework [43] according to the event spectra given in Eqs. (3), (6) and (8) of Ref. [34], where an energy resolution of 3%/ E o /(MeV) is always adopted and the analytical SN neutrino flux model with E ν e = 12 MeV, E ν e = 14 MeV and E ν x = 16 MeV is assumed.

III. STRATEGY FOR RECONSTRUCTION
As demonstrated in Ref. [34], one can reconstruct the SN neutrino spectra of all flavors with three unfolded spectra of IBD, pES and eES channels in one single large LS detector via a simple bin-to-bin separation method. However, as the number of events for ν e dominates over those for all other flavors in the eES channel, it has been assumed in Ref. [34] that the detector response matrix between the initial neutrino energy and the observed event energy is universal for all neutrinos. In this section, we shall improve the strategy for reconstruction by relaxing this assumption, and extract the energy spectra of all flavor neutrinos directly from the observed spectra of these three channels by using the overall response matrix.

A. The separated analysis
For comparison, we first summarize the basic idea of Ref. [34] for the spectral reconstruction. For the IBD channel, the observed event spectrum can be calculated as follows where N p is the number of protons in the LS detector, F ν e stands for the spectrum of the SN ν e , and S IBD denotes the observed event spectrum. In addition, σ IBD ν e is the cross section matrix for the IBD reaction, which maps the initial neutrino energy E ν into the visible energy of the final-state particle E v . In Eq. (2), D IBD is the probability matrix, which accounts for the detector effects and converts E v into the observed energy E o . Note that those matrices and vectors have been constructed by discretizing the relevant continuous energy spectra into finite energy bins. To be explicit, we can recast Eq. (2) into the matrix form where n o , n v and n ν are the number of bins for the observed energy, the visible energy of the final-state particle and the initial energy of SN neutrinos, respectively. Comparing between Eq.
(2) and Eq. (3), one can easily identify the definitions of the corresponding matrices.
For practical purposes, we further normalize the cross-section matrix σ IBD νe and the initial spectrum F νe in the following way  . . .
where σ i for i = 1, 2, · · · , n ν stands for the total cross section for the incident neutrino with the central energy E i ν of the i-th energy bin. Eq. (4) can be rewritten in a more compact form as R IBD F IBD ν e = S IBD with R IBD ≡ D IBD σ IBD , where the definitions of relevant matrices are self-evident. Therefore, the reconstruction of the cross section weighted neutrino spectrum F IBD ν e from the observed event spectrum S IBD can be regarded as a linear inverse problem, which can be routinely solved with an unfolding method [34]. Here R IBD is exactly the detector response matrix in Ref. [34], which can be built with a large number of simulated events and is independent of SN neutrino models.
As for the pES and eES channels, the observed event spectra receive contributions from neutrinos and antineutrinos of all flavors, and can be calculated in a similar way where N p(e) denotes the number of protons for pES (electrons for eES), and α refers to different neutrino flavors. In analogy to the IBD channel, we can deal with the pES channel by introducing the detector response matrix R pES ≡ D pES σ pES and the cross section weighted spectrum α F pES α . Then, it is straightforward to extract α F pES α from the observed event spectrum R pES α F pES α = S pES . Such a treatment is quite reasonable in the pES channel, since neutrinos of all flavors interact with protons via the neutral-current interaction and the corresponding cross sections are almost the same. However, this is obviously not the case for the eES channel, for which the cross section of ν e -e − scattering is about two times larger than that of ν e -e − scattering and six times larger than that of ν x -e − scattering. In Ref. [34], it has been assumed that ν e dominates over all other flavors in eES channel and a universal response matrix R eES ≡ α D eES σ α eES is then used to achieve the reconstruction of the weighted true spectrum α F eES α . This approximation is only valid for the one flavor dominated case. Finally, the energy spectra for different flavor neutrinos can be simply separated bin-by-bin from the unfolded spectra F IBD ν e , α F pES α and α F eES α . Although the cross section of the eES for ν e is much larger than that for ν x and thus ν e gives rise to most of the eES events, the summation of the contributions from four flavors of ν x with high energies is not negligible. Due to the quenching effects on recoiled protons, only the high-energy part of the ν x spectrum (i.e., above 20 MeV) can be really reconstructed from experimental data. The eES channel will be able to provide useful information about the low-energy part of the ν x spectrum. Therefore, it is interesting to have a further look at the reconstruction of neutrino spectra in the eES channel.

B. The combined analysis
Now we put forward a combined analysis of all three observed spectra from IBD, pES and eES by grouping the multi-flavor neutrino spectra into an overall neutrino spectrum, and likewise for the event spectra. More explicitly, we have where S c is the whole event spectrum with (n IBD o + n pES o + n eES o ) bins, and F c is the overall neutrino spectrum with (n νe + n ν e + n νx ) bins. Therefore, the whole event spectrum S c , as observed in the LS detector, can be described in one single equation  where the (n IBD o + n pES o + n eES o ) × (n νe + n ν e + n νx ) matrix on the left-hand side is just the detector response matrix R c for the combined analysis, and the sum in its third column is running over ν µ , ν τ and their antiparticles. Some explanations for the structure of R c are necessary. First of all, the block matrices D IBD σ IBD νe and D IBD σ IBD νx in the first row of R c are actually vanishing, since the IBD reaction takes place only for ν e . Second, for the pES and eES channels, the response matrices have been determined by simulating the interactions of ν e , ν e and ν x with the target particles in the detector. In particular, in the eES channel, the elastic scattering of ν e , ν e and ν x with electrons has been investigated individually to produce three matrices in the last row of R c . Third, it is worthwhile to emphasize that such a combined analysis treats ν e , ν e and ν x on the same footing and it is independent of SN neutrino models. Moreover, this method can be easily extended to include the observations from the WC and LAr-TPC detectors, leading to a global analysis of all SN neutrino data. This can be realized by adding a new row into R c , which is determined by the specified interaction channel in a given detector, and accordingly a new row in S c .
To demonstrate how the combined analysis works, we concentrate on the LS detector.
Given Eq. (7), one can immediately apply the spectral unfolding approach to extract F c .
See, e.g., Refs. [44,45], for a general review on the unfolding problem in particle physics.
In the following, we implement the Singular Value Decomposition (SVD) method with a proper regularization scheme, as proposed in Ref. [46], to reconstruct SN neutrino spectra.
The regularization is important to suppress the spurious oscillating components in the final results. The main strategy for reconstruction is outlined as below: • The detector response matrix for the combined analysis should be constructed, and the results are now depicted in Fig. 2. Comparing Fig. 2 with the matrix R c in Eq. Afterwards, all the involved bins will be multiplied by the corresponding total cross section σ at E ν and the number of target particles in the specified channel. Finally, all nine histograms generated in this way are shown in Fig. 2. Note that the three block matrices for eES channels are different, which demonstrate the limitation of the separated analysis as discussed in the previous subsection.
• After preparing the response matrix, we then apply the SVD unfolding method to reconstruct of the SN neutrino spectra. In the unfolding process, the binning scheme for the observed event spectrum in each channel depends on the event statistics and should be carefully handled in order to guarantee a comparable number of neutrino events in each bin of the observed energy, as illustrated in Fig. 1 for a SN distance at 10 kpc. For each true neutrino energy spectrum, we employ the equal-size binning scheme but combine the bins at the boundaries due to the limited statistics. The actual binning scheme for the realistic 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 TSVDUnfold in ROOT, where a proper regularization parameter is set for this work.
As we have mentioned, the response matrix has been built from the simulated neutrino events in each detection channel for a given detector. Therefore, it depends only on the experimental setup and neutrino interactions with the target particles, implying that the combined analysis can be applied to any SN neutrino model, namely, both analytical and numerical ones, and the realistic SN explosion. In the next section, we shall present the final results of the reconstructed SN neutrino spectra.

IV. RECONSTRUCTION OF NEUTRINO SPECTRA
Following the approach of combined analysis in the previous section, we reconstruct the energy spectra of all flavor SN neutrinos by using the simulated events described in Section II. For comparison, the observed spectra are simulated for a SN distance at 10 kpc, 1 kpc, and 0.2 kpc. In addition, we investigate the impact of the energy threshold and the regularization parameter on the reconstruction. To quantitatively assess the model dependence, we also calculate the bias distributions of the reconstructed spectra by repeating the reconstructions for several numerical SN models from the Japan group. In the end, we explain how to reconstruct the original spectra of SN neutrinos in the presence of neutrino flavor conversions.

A. Results for analytical models
To generate the SN neutrino event spectra, we consider the KRJ parametrization of SN MeV. The observed spectra in the pES and eES channels are divided into equal-size bins in a logarithmic scale, while that for IBD is equally binned in a linear scale. The number of bins for IBD, pES and eES is simply taken to be the same, i.e., 20, 30 and 40 bins for the SN at 10 kpc, 1 kpc, and 0.2 kpc, respectively, which can be read out from Fig. 3. Then, the SVD unfolding method is implemented to extract the true neutrino spectra F c from the whole observed spectra S c .
The reconstructed SN neutrino spectra are shown in Fig. 3, where three rows correspond to the case of SNe at 10 kpc (upper panels), 1 kpc (middle panels), and 0.2 kpc (lower panels). The plots in the left column of Fig. 3 focus on the high-energy parts (i.e., above about 20 MeV) of the neutrino spectra, while those in the right column represent the full spectra. In all the plots, the solid histograms stand for the true SN neutrino energy spectra while the points denote the reconstructed spectra. The vertical error bars attached to the points indicate the statistical uncertainties, arising from those of the observed SN neutrino spectra, and the horizontal ones show the bin widths. In the cases of the three SN distances, the results are in general better than those obtained in Ref. [34] via a simple bin-by-bin separation method. It is straightforward to understand why the combined analysis works better. First, the cross sections of different flavor neutrinos in eES channel are treated accurately while there is an approximation of ν e -domination in Ref. [34]. Second, taking account of the correlation among different reaction channels, the combined analysis is able to reduce the meaningless fluctuations in the simple bin-by-bin separation procedure. For the SN at 1 kpc, the best precision for the ν e , ν e and ν x spectra can reach the level of 1%, 20% and 6%, respectively. It is obvious that the precision gets better when the SN distance becomes smaller and thus the statistics turns out to be larger.
It is worthwhile to note that the reconstructed ν e and ν x spectra in the right column of Fig. 3 deviates significantly from the true spectra below 20 MeV and they are anti-correlated between each other. The main reason for this behavior is that such low-energy neutrinos can only produce protons with low recoil energies, which after the quenching effects will be lying below the threshold E th o = 0.2 MeV of the observed energy. Therefore, the pES channel at the LS detector is only sensitive to SN neutrinos with energies above 20 MeV or so. Then it is clear that the deviation is a systematic bias of the energy threshold and the ν e and ν x spectra below 20 MeV can be constrained but cannot reconstructed accurately through one single observed spectrum in the eES channel. We shall give a detailed discussion on the effect of different energy thresholds in the next subsection.

The impact of energy threshold
To check the impact of E th o on the reconstruction, we generate another trial of events with E th o = 0.01 MeV, which is just taken for illustration and will certainly be impossible to achieve in the present and next-generation LS detectors. For the SN at 1 kpc, the unfolded neutrino spectra are shown with orange points in Fig. 4. MeV is about 6.7 MeV, which can also be seen from the spectra of ν e and ν x of Fig. 4. The reconstructed spectra of ν e and ν x are well consistent with the true spectra for bins above 6.7 MeV. The comparison between these two cases of different values of energy threshold clarifies that a lower threshold of the observed energy can help to extract the ν x spectrum accurately in a wider range of energies.

The impact of regularization parameters
As usual, the regularization parameter is introduced in the unfolding procedure in order to suppress the spurious oscillatory components which arise from the direct computation of the inverse of the response matrix. In doing so, the bias will be unavoidably brought into the final results at the same time. To optimize the value of the regularization parameter, one should make a balance between the reduction of spurious oscillatory components and the bias. Given the KRJ parameterization, we simulate 500 trials of neutrino events with the SN at 1kpc. Then the combined method is applied to reconstruct the neutrino spectra with a high, medium and low values of the regularization parameter respectively. 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 reconstruction bias from the unfolding method. In order to reduce the statistical fluctuation, we introduce the mean bias b M over 500 trials and its standard deviation σ M , namely, where k is the index for the trial and N = 500 is the number of total trials for the KRJ parameterization.
In parameters. As the regularization parameter becomes larger, the spurious oscillations in the unfolded spectra are getting more suppressed, but the resultant biases are also more sizable.
In the current study, we do not attempt to obtain the optimal value of the regularization parameter for each trial, but directly choose a suitable value to keep the statistical errors of the unfolded spectra smaller than or comparable to the induced bias.

B. Results for numerical models
To demonstrate the robustness and the model-independence of the combined method, we implement the SN neutrino fluences from twelve numerical models simulated by the Japan group [37]. In addition, the same calculation is performed again but with a response matrix constructed by using the SN neutrino spectra from this numerical model. The strategy for reconstruction is the same as that adopted for the analytical models. Therefore, for each trial in a specified numerical model, we have the results of unfolded neutrino spectra.
In Fig. 6, the distributions of the mean bias as defined in Eq. 8 for ν e , ν e and ν x are shown in the upper, middle and lower row, respectively and the dashed horizontal lines are shown   for the best precision of the reconstructed spectra of different flavor neutrinos for the SN at 1 kpc in Fig. 3. The left column summarizes the results with the fixed response matrix from the flat model, while the right column is for results with the response matrix from the corresponding numerical model. The same value of the regularization parameter is adopted as that for the SN at 1kpc in Fig. 3. Some important conclusions can be drawn. First, the bias distributions for different numerical models, as denoted by the colored curves in Fig. 6, are quite similar to each other, implying that this combined analysis is robust and modelindependent. Second, the distributions in left and right column are also well consistent.
Although the bias of ν e at a few low-energy bins seems model dependent, the reconstruction for spectra above 20 MeV is robust and not affected by the initial neutrino spectra used to build the response matrix. Third, as seen from those plots in Fig. 6, the averaged bias for the central part of energies in all different flavors is smaller than the statistical errors shown by the dashed horizontal lines. This observation indicates that the adopted regularization parameter is very reasonable.

C. Neutrino flavor conversions
In the previous discussions, the flavor conversions of SN neutrinos have been completely ignored. When propagating outward from the neutrino sphere, SN neutrinos may experience spectral splits or swaps from the collective neutrino oscillations [12][13][14][15][16][17][18][19][20] due to the dense neutrino background. Several recent studies [47][48][49] show that the self-induced fast oscillations may take place close to the neutrino sphere. However, it remains an open question whether the collective neutrino oscillation do happen in the real SN environment [20]. In the mantle of the SN, the MSW matter effects will play an important role, resulting in a partial or complete conversion between ν e and ν x (or between ν e and ν x ) depending on the neutrino mass ordering. On the way to the Earth, although SN neutrinos have lost quantum coherence [50,51], there will be regeneration effects due to the Earth matter [52][53][54][55][56][57][58]. Therefore, the final SN neutrino spectrum of a given flavor entering into the detectors is actually a mixture of the initial spectra of different flavors.
To explain how to deal with neutrino flavor conversions in the reconstruction of SN neutrino spectra, we take into account the MSW matter effects [52], for which the overall picture is clear and well understood. Once the collective oscillations of SN neutrinos are established, it will be straightforward to incorporate them into our analysis. According to Ref. [52] and the latest neutrino oscillation data [59], we can find for the normal neutrino or for the inverted neutrino mass ordering (IO):  where F α stands for the initial neutrino spectra as given by the KRJ parametrization while F ′ α is the spectra after flavor conversions. The conversion matrix between the initial and final spectra will be denoted as C, in which the neutrino mixing angle θ 12 ≈ 33 • is taken from Ref. [59].
In consideration of neutrino flavor conversions, the observed event spectra S c in the LS detector shall be written as where F ′ c = C ·F c with F c being the initial neutrino spectra. Since the pES channel is subject only to the neutral-current interaction, the event spectrum is not affected by neutrino flavor conversions. However, the event spectra in the IBD and eES channels will differ from those in the scenario without flavor conversions. To extract the initial SN neutrino spectra, we can apply the combined method as well if the response matrix R c is convolved with the flavor conversion matrix C. The results of such a convolution are shown in Fig. 7.
Two trials of SN neutrino events in the IBD, pES and eES channels for a JUNO-like detector are simulated, including neutrino flavor conversions as in Eq. (10) and Eq. (11) for NO and IO, respectively. The initial neutrino spectra are described by the KRJ parametrization with E ν e = 12 MeV, E ν e = 14 MeV and E ν x = 16 MeV for a SN distance of 1 kpc. Then, following the combined analysis, we reconstruct the initial SN neutrino spectra and show the final results in Fig. 8. For the IO case, the ν x spectrum can be well reconstructed because F ′ ν e = F ν x and ν e at the detector is precisely measured in the IBD channel. While for the NO case, the IBD channel constrains the spectra of both ν e and ν x because F ′ ν e = cos 2 θ 12 F ν e + sin 2 θ 12 F ν x . Therefore the spectral accuracies of ν e and ν x are not as good as that of ν x in the IO case. When a more complicated scenario of neutrino flavor conversions is considered, one can replace the flavor conversion matrix C with the new one and repeat the analysis to extract the initial SN neutrino spectra. Therefore, our strategy for reconstruction is useful to test the pattern of SN neutrino flavor conversion. Note that the current study of neutrino oscillation effects is only included in the time-integrated neutrino spectra reconstruction. The scenario of neutrino flavor conversions might be different for different phases of the SN neutrino burst. We want to stress that the method proposed here can still be applicable for the reconstruction of the time-dependent neutrino energy spectra, but one needs to suffer from the relatively lower event statistics.

V. SUMMARY
For a future galactic core-collapse SN, we have proposed a model-independent approach to reconstruct all flavor SN neutrino energy spectra by performing a combined analysis of the IBD, pES and eES detection channels in a 20 kiloton JUNO-like LS detector. First of all, we briefly recap the calculation of SN neutrino events in the LS detector and the separated method used in Ref. [34] to reconstruct SN neutrino spectra, where however the response matrix for different neutrino flavors in eES channel is not fully considered. Then, the combined method is introduced to treat all neutrino flavors ν e , ν e and ν x in three channels on the same footing, and applied in the spectral reconstruction with the simulated SN neutrino events. Similar calculations have been carried out for the SN at different distances (i.e., 10 kpc, 1 kpc, and 0.2 kpc). In addition, we investigate the impact of the threshold energy of the detector and the regularization parameter of the unfolding method on the spectral reconstruction. The combined method is demonstrated to be robust and modelindependent via the analysis of both analytical and numerical neutrino data. Finally, taking account of neutrino flavor conversions under the MSW matter effects in the SN mantle, we explain how to implement the combined analysis to extract the initial neutrino spectra in the presence of flavor conversions.
Although we have concentrated on the spectral reconstruction with a LS detector, the detections from the large WC and LAr-TPC detectors should be utilized to probe the SN neutrino spectra globally [60]. It is intuitively convenient for the combined method to accomplish this task. What one has to do is just to extend the the response matrix and the observed event spectra with the information from other detectors and the relevant detection channels. In this case, the low statistics in the eES channel of the LS detectors can be compensated by the WC and LAr-TPC detectors. On the other hand, the advantage of the low energy threshold of LS detectors is maintained to reconstruct the ν x spectrum.
The combined analysis of the LS, WC and LAr-TPC detectors in the reconstruction of SN neutrino spectra is very interesting and deserves another dedicated study. Moreover, this method can also be used to reconstruct the spectra of solar neutrinos and ultrahigh-energy cosmic neutrinos when a large statistics in the multi-flavor detection is accumulated.