Developing an end-to-end simulation framework of supernova neutrino detection

Massive stars can explode as supernovae at the end of their life cycle, releasing neutrinos whose total energy reaches $10^{53}$ erg. Moreover, neutrinos play key roles in supernovae, heating and reviving the shock wave as well as cooling the resulting protoneutron star. Therefore, neutrino detectors are waiting to observe the next galactic supernova and several theoretical simulations of supernova neutrinos are underway. While these simulation concentrate mainly on only the first one second after the supernova bounce, the only observation of a supernova with neutrinos, SN 1987A, revealed that neutrino emission lasts for more than 10 seconds. For this reason, long-time simulation and analysis tools are needed to compare theories with the next observation. Our study is to develop an integrated supernova analysis framework to prepare an analysis pipeline for treating galactic supernovae observations in the near future. This framework deals with the core-collapse, bounce and proto-neutron star cooling processes, as well as with neutrino detection on earth in a consistent manner. We have developed a new long-time supernova simulation in one dimension that explodes successfully and computes the neutrino emission for up to 20 seconds. Using this model we estimate the resulting neutrino signal in the Super-Kamiokande detector to be about 1,800 events for an explosion at 10 kpc and discuss its implications in this paper. We compare this result with the SN 1987A observation to test its reliability.


Introduction
A massive star leads to a spectacular explosion called a supernova (SN) at the end of its life. The supernova releases 10 51 erg in the form of ejected matter. It also emits 10 53 erg, which amounts to 10% of the solar rest mass energy from gravitational binding energy, in the form of neutrinos within just a few tens of seconds. Supernova explosions are among the most energetic phenomena in the universe, releasing into space both elements that were synthesized over the entire life of the parent star as well as heavy elements synthesized during the explosion itself. For this reason, supernovae provide important information about the chemical evolution of the universe.
In order to study the supernova explosion mechanism itself, it is important to study the neutrinos emitted during the evolution of core-collapse supernovae. The standard scenario of core-collapse supernovae is believed to be the so-called neutrino-driven explosion, in which the interaction of neutrinos and matter plays an essential role [1][2][3][4]. Neutrinos are involved in all phases of the explosion through trapping, emission and absorption in the collapse, as well as during the bounce and explosion resulting in the formation of a neutron star. This general scenario was shown to be viable by the detection of neutrinos from SN 1987A [5][6][7]. However, the detailed mechanism has not been clarified because only about 20 neutrino events were observed from SN 1987A, which was only able to provide information on the scale of the total energy budget and pointed to neutrino diffusion over long time scales from a compact object. Since that time neutrino detectors have been waiting to detect the next galactic supernova. For example, Super-Kamiokande [8], which is a large water Cherenkov detector, is expected to detect thousands of events from a supernova occurring in our Galaxy. Once such neutrinos are detected with high statistics, they will provide valuable information that can elucidate the full dynamics of core-collapse supernovae and the subsequent formation of neutron stars.
It is therefore important to prepare various sets of theoretical models and their predicted signals to allow for rapid and thorough exploration of the possible physical scenarios using any observed supernova neutrinos. Studies of supernova neutrinos covering their evolution from the explosion to the neutron star birth have so far been limited to a few explosion models. Although there are many numerical studies of supernovae, most simulations aim to understand the supernova mechanism focusing on only the first second of the process, which is critical to determine whether a supernova explodes at all. Accordingly, there are only a few studies of future supernova neutrino detection and they typically adopt a particular model which successfully explodes and forms a neutron star [9][10][11][12][13][14][15]. There are also other studies on the number of expected events based on simulations with prescribed explosions [16,17] as well as approximate analytic solutions [18]. Without sufficient simulations covering the entire evolution of the supernova explosion, neutron star formation, and the corresponding neutrino event prediction in a terrestrial detector, it is difficult to compare theoretical predictions with a real observation to obtain constraints on model parameters.
To overcome these difficulties, recent supernova studies predict neutrino signals in terrestrial detectors [17,[19][20][21] and to the same end we are developing a framework for 2/23 supernova-neutrino analysis (Fig. 1). Our framework consists of a supernova simulator, a detector simulator and a supernova analyzer. The supernova simulator uses a package of simulation software to provide theoretical predictions of the neutrino emission over long time scales, covering the core-collapse, bounce, explosion, and cooling phase of the proto-neutron star (PNS). The detector simulator provides mock data representing detected neutrino events using the output of the supernova simulator. The analyzer aims to provide methods to compare the mock sample and the observational data to quickly analyze the properties of a supernova based on the results from the supernova and detector simulators. The purpose of our framework is to connect the numerical simulations of the stellar collapse and explosion with the prediction and observation of neutrino events at detectors on earth. We want to bridge the gap from the simulation to the detector so that one can analyze a burst in a systematic manner and to eventually extract model information from the supernova explosion.
The first steps of the framework development are the creation of a supernova simulator that performs long-time simulations of over 20 seconds and a detector simulator that generates mock data by Monte-Carlo sampling from the output of the supernova simulator. We provide the supernova analyzer to make a rapid analysis of a real supernova burst and to find connections between such an observation and theory using the evolution of the number of events and their energy as a function of time. Note that previous studies have also attempted to produce the late time neutrino spectra [21]. In the present work we additionally discuss the mean and variance of the number of events and their energies in order to provide robust information even when statistics are small, which may be beneficial immediately following a supernova observation.
In this paper, we demonstrate the flow of our framework by showing results of the supernova simulation described below and mock data generated therewith. We describe the new long-time simulation in §2 and §3, which correspond to the "SN simulator" in Fig. 1. We show examples of the construction of mock data sampled from our simulation and analyses of event rates at Super-Kamiokande in §4, corresponding to the "Detector simulator" and "Analyzer" in Fig. 1, respectively.

Hydrodynamics simulation Methods
We employ GR1D [22,23] for hydrodynamic simulations in which general relativistic hydrodynamics equations, as well as multi-energy neutrino radiation transport equations, are solved in a spherically symmetric geometry. 1 Neutrino transport is solved using the truncated momentum formalism of Ref. [24].
Since GR1D is open source it has the advantage of allowing the community to reproduce and extend our results, which is not necessarily the case with other custom codes. 2 We employ the nuclear equation of state (EOS) based on the density-dependent relativistic 1 The code is publicly available at https://www.GR1Dcode.org. 2 The truncated momentum formalism used in GR1D is likely to produce qualitatively reasonable results and we expect that our conclusions in this work will not change significantly even if we employ more sophisticated numerical schemes for computing the details of neutrino transfer. In particular GR1D addresses three flavor neutrino transport while recently calculations using six flavors have also been performed [25,26]. The additional flavors are likely to influence the initial collapse and may also affect the late stage neutrino spectra and will therefore be considered in a future work.  Fig. 1 Schematic diagram of the integrated supernova framework. This framework comprises a long-time supernova simulator, a detector simulator and an event analyzer. Before detection of a supernova, we provide mock samples to allow for practice analysis. After detection, we can quickly analyze the supernova to provide feedback to various supernova models.
mean-field (DD2) model [27]. Neutrino interactions are calculated using a numerical table made with NuLib 3 [23]. We employ a 9.6 M zero-metallicity progenitor provided by A. Heger (2016, private communication), which has been used in the previous works and shown to explode even assuming spherical symmetry [28,29]. This allows us to perform a long-time simulation from the initial core-collapse through the PNS cooling without introducing other phenomenological modeling.
Multi-dimensional effects such as the standing accretion-shock instability [30,31], rotation [32] and convection [33] also play important roles in supernovae and are expected to influence the neutrino signal somewhat. However in this study, we do not consider multi-dimensional effects due to their computational expense. Indeed, in the case of a light progenitor around 9.0 M , multi-dimensional effects are expected to have little impact on the neutrino signal [30] so we ignore them here. Note that there are methods under development that approximately incorporate multi-dimensional effects into spherically symmetric simulations [20,34].
In the following we explain our modifications of GR1D for long-term simulations in §2.1, the grid settings in §2.2, and the neutrino reactions in §2.3.

Modifications for long-term simulations
GR1D is publicly available software for the simulation of core-collapse supernova explosions that is well-suited to modelling the accretion phase, which typically occurs up to one second after the bounce. Since we are interested in the PNS cooling phase in this paper, which happens later, we have modified the original GR1D in two ways. First, the original grid profiles of GR1D are optimized to resolve the supernova shock evolution, which is not suitable for the PNS cooling phase. To perform the PNS cooling simulation we have adjusted the grids so as to resolve the steep density gradient at the surface of the PNS. Details are presented in the next section. We also change the CFL number from the original setting of 0.5 to 0.25 starting 8 s after the bounce. Second, we find that the long-term evolution causes some thermodynamic quantities such as the density, temperature, electron fraction, and η = µ/k B T , with µ being the chemical potential of electrons, to exceed the bounds of GR1D's numerical tables. To avoid such overflows, we have reproduced and extended those tables to cover a wider range of parameters. Further, we have chosen to fix any values that exceed the limits of the new tables to their closest extremum. Details are given in §2.3.

Grid settings
During PNS cooling, the surface density gradient is extremely steep, decreasing from ∼ 10 14 g cm −3 to ∼ 10 8 g cm −3 within 20 km (see the red line in Fig. 2). Therefore we have employed fine grids to resolve this gradient. For this purpose we use custom2, which allows the 30 innermost zones to have zone widths that decrease logarithmically from 1 km to 0.1 km and allows the intermediate 78 zones to have a constant width (0.1 km) with an otherwise logarithmic progression that extends up to 5000 km (see the blue line in Fig. 2). There are 300 total grid points used throughout the simulation. As shown in Fig. 2 the finest resolution grid is assigned for the steepest density gradient at the PNS surface.

Neutrino reaction table
In GR1D the neutrino reactions are given by a numerical table that is computed in advance using NuLib and DD2 as an EOS input. Since the original table does not cover the whole thermodynamic range necessary for the PNS cooling phase, we have produced an expanded 5/23 table. It spans ρ = 10 6−15.5 g cm −3 with 82 logarithmically-spaced points T = 0.05 − 150 MeV with 65 logarithmically-spaced points, Y e = 0.015 − 0.55 with 82 linearly-spaced points, and η = 0.1 − 100 with 61 logarithmically-spaced points.
Neutrino interactions with matter are an integral part of the explosion's evolution. In this work we consider the emission and absorption reactions where ν e ,ν e , p, n, e − , e + , and (A, Z) represent electron-type neutrinos, electron-type antineutrinos, protons, neutrons, electrons, positrons, and a nucleus with mass number A and atomic number Z, respectively. These reactions are calculated based on [35] with weak-magnetism and recoil corrections from [36] taken into account. Neutrino absorption on heavy nuclei is also based on [35,37]. Similarly we consider elastic scattering via the following reactions where α is the helium nucleus, ν indicates that the reaction is insensitive to the neutrino flavor, and ν i indicates that the reaction depends on flavor. These are also based on [35,37]. Inelastic scattering has been calculated following [37]. We consider the following thermal processes, which have been calculated by [35,37]. Note that electron-positron annihilation is considered only for ν x throughout the simulation, where ν x refers to non-electron-type flavors. Including electron flavors results in the simulation stopping prematurely around the bounce. Note though that this reaction has little influence on the late phase neutrino spectra [23]. In addition, the rate of electron-positron annihilation is a few orders of magnitude lower than that of nucleon-nucleon bremsstrahlung at energies lower than 50 MeV [38]. Though the original GR1D includes nucleon-nucleon bremsstrahlung for ν x only, we find that without including electron-type neutrinos the average energies of these flavors remain constant at late times, though they are expected to decrease based on physical considerations. Accordingly, we have added them to the bremsstrahlung process to resolve this issue. Note that we do not consider medium suppression for simplicity [38].

Hydrodynamic simulation results
In this section we present results from our simulation, summarizing the dynamical features in section 3.1 before discussing neutrino luminosities and average energies in section 3.2. 6/23 Figure 3 shows the trajectories of several mass shells and the shock wave in red and black, respectively. The shock trajectory clearly shows a successful shock propagation. The innermost red line corresponds to the mass coordinate 1.20 M , while the outermost thick line indicates that of 1.37 M . Here the thick red lines are separated by 10 −2 M . In order to see the mass shells outside the PNS in detail, thin red lines for 1.330 M to up to 1.374 M are shown in intervals of 10 −3 M . At the beginning of the simulation the mass shells fall into the center gradually and then rapidly just before the core bounce. The initial central density of the progenitor is 10 9 g cm −3 , and the core collapse continues until the central density exceeds the nuclear saturation density ∼ 10 14 g cm −3 . The core bounce of the z9.6 progenitor occurs 0.254 s after the start of the simulation.

Dynamical properties
After the core bounce, the mass shells inside 1.36 M accrete onto the PNS surface, while others are ejected after the passage of the shock. Hence the baryonic mass of the PNS remnant of this model is 1.36 M . The shock wave continues to propagate outward until it reaches the simulation boundary of 5000 km 0.35 s after the core bounce.
Here we define the commonly used diagnostic explosion energy [39] as where is the binding energy, α = √ −g 00 is the lapse function (see O'Connor et al. (2010) [22] for details), ρ is the rest-mass density, is the specific internal energy, P is the pressure, W is the Lorentz factor, and dV is the three-volume element for the curved space-time metric. With this definition the diagnostic explosion energy of this model converges to 4.2 × 10 49 erg, which is significantly smaller than the observed value of 10 51 erg [40]. However, this value is consistent with that of previous studies that used the same progenitor model [28,29]. Observations of SNe indicate that they have a rather broad distribution of explosion energies, from ∼ 10 50 erg to 10 52 erg [41]. The explosion energy is naively determined using the binding energy of progenitor layers in the vicinity of the final remnant mass of the compact object so that stars with small binding energies are possible candidates for weak explosions (see e.g., [42]).
Radial profiles of the entropy and Y e at different times are shown in Fig. 4. Before the bounce at −100 ms (blue dotted line), the entropy is as low as 1.2 k B /baryon and is constant with respect to the mass coordinates. The electron fraction at that time is constant and as high as 0. Comparing our results with the simulations in Refs. [28] and [29], we note that the shock wave of our simulation propagates outward without slowing down, while the shock waves in previous studies slow slightly at radii around 2 × 10 7 cm (Ref. [28] indicates stronger deceleration than Ref. [29]). The diagnostic explosion energy (∼ 4.2 × 10 49 erg) observed in our simulation is close to those of previous studies (∼ 2 × 10 49 erg in Ref. [28] and ∼ 1-4×10 49 erg in Ref. [29] for 1D simulations) and similarly our resultant PNS baryonic mass, 1.36M , is also close to that of Ref. [29] (1.38M ). Thus, we conclude that overall behaviors of the simulations are consistent with one another. We speculate that small differences are due to the differences in the adopted microphysics modeling, such as the EOS or neutrino interaction model.

Neutrino properties
We simulated the time evolution of neutrinos in the z9.6 model. The luminosity L ν , average energy E ν , and the root-mean-square (RMS) energy E 2 ν are shown in Figs. 5, 6, and 7, respectively. In this paper, the definitions of these energies are where dN (E)/dE is the neutrino number density per unit energy. In Fig. 5, the electron neutrino (ν e ) luminosity begins to rise at ∼ 20 ms before the core bounce, and then falls off temporarily around the bounce due to neutrino trapping. Note that there are few neutrinos of the other types before the core bounce since the electrons are degenerate at the lower temperatures at this stage in the evolution. During this period the dominant reaction in detectors on the ground would therefore be electron scattering, which is sensitive to the direction of the incoming neutrino and hence the SN location on the sky. Shortly after the core bounce, the ν e luminosity rapidly increases to 6.5 × 10 53 erg s −1 due to the neutronization burst and then drops to ∼ 1.0 × 10 53 erg s −1 . Note that the peak of the ν e luminosity from GR1D is likely to be higher than that of other simulations according to O'Connor et al. (2018) [43]. The anti-electron neutrinoν e luminosity gradually increases to the same level as the ν e luminosity after the core bounce. For non-electron-type neutrinos, ν x , which represents ν µ , ν τ , and their anti-particles collectively, the luminosity quickly jumps to 1.2 × 10 53 erg s −1 . Each of the flavors represented by ν x has the luminosity shown in the figure. After 100 ms, all luminosities gradually decrease and this trend continues during the PNS cooling phase. Finally, all luminosities have nearly the same value at 20 s. Figure 6 shows that the average ν e energy is 8 MeV initially, increases to 10 MeV, and then drops slightly during neutrino trapping. At the core bounce the average ν e energy reaches a peak value of 15 MeV. After the core bounce the average ν e energy is almost constant around 10 MeV and lasts for several hundreds of milliseconds. The averageν e energy after the core bounce is also constant but higher than that of ν e becauseν e interacts with matter more weakly than ν e and theν e neutrinosphere is at a smaller radius than that of ν e . The 9/23 average ν x energy is between 5 MeV and 3 MeV initially, rises to 17 MeV rapidly at the core bounce, and then slightly decreases to 15 MeV. Afterward, the average energies for all species gradually decrease to 6 MeV. In particular, the average energies ofν e and ν x energies almost overlap at late times. In comparison to the simulation of Ref. [44], the average energy of ν e in our simulation is higher by 0.5 MeV and that of ν x is lower by 0.5 MeV in the early phase. These small differences are likely due to differences in the microphysics as mentioned above §3.1. However these differences are not expected to have a large impact on response of the detector.
As shown in Fig. 7, the behavior of the RMS energies is similar to that of the average energies for all types of neutrinos. The RMS energies are related to the pinching parameter of the neutrino spectra. Both the average and RMS energies are translated to the observed distributions in a terrestrial detector, indicating observations will provide feedback on their initial distributions at the time of the SN. Note that there is a small ripple at 110 ms for all species because the shock wave crosses the point where we calculate the neutrino luminosities, 500 km, at that time. Figure 8 shows the total energy of the emitted neutrinos, that is the luminosity in

Detection simulations
In this section we demonstrate the detector simulator part of our analysis framework (Fig.  1). So as to produce mock data samples we perform Monte Carlo simulations of the neutrino 10 Super-Kamiokande is the largest underground water Cherenkov detector, consisting of a stainless tank whose height is 41.4 m and whose radius is 39.3 m, and is filled with 50 kton of ultra-pure water. The tank is separated into two regions, an inner and an outer detector. The inner detector views the water target with 11,129 photomultiplier tubes and the nominal fiducial volume is a 22.5 kton volume therein, which is defined to reduce backgrounds from the detector walls and ensure proper event reconstruction. However, for supernova neutrino bursts, many neutrino events will be detected over a time interval short enough that the 11  contribution from background events is expected to be negligible. Thus, hereafter we assume the entire inner detector volume of 32.5 kton can be used for observation. Since this study is the first step in the construction of our analysis framework, we do not include the details of the detector response for simplicity. Although in the supernova simulations we do not include neutrino oscillations, for the detector simulation we modify the fluxes to account for neutrino flavor oscillations in vacuum and in matter. Following Refs. [45,46], we mix neutrinos as, for the normal hierarchy and for the inverted hierarchy, where F ν and F ν are neutrino fluxes after and before neutrino oscillations, respectively, and p is 0.69. Note that we employ adiabatic flavor conversion for the following reason. The density of the progenitor in this study decreases rapidly outside the core and the region where the resonant (non-adiabatic) conversion can occur is in the vicinity of the iron core ( 10 9 cm). Although the progenitor density scale height is large enough to satisfy the adiabatic condition, the shock wave is expected to change the density structure and lead to non-adiabatic conversion once it reaches the resonant region. With a typical shock velocity of ∼ 10 9 cm s −1 , the shock passes through the resonant regime within 12/23 a few seconds. Thus, resonant conversion does not change the following discussion, which focuses on later timescales.

Detection properties
Here we focus on the number of events in SK, their angular distribution, and the time evolution of the event rate. Events are simulated assuming supernova distances of 1 kpc, 5 kpc, 10 kpc and 50 kpc. Our Monte Carlo simulations are performed as follows: (i) calculate the expected number of events in each time interval, which are chosen to be 0.01 s (t < 0.744 s, where t is the time measured relative to the bounce) and 0.1 s (t > 0.744 s), (ii) determine the time of each event in the time interval at random, (iii) evaluate the event spectrum in the time interval and (iv) determine the event energy according to that spectrum for each event.
We consider the inverse beta decay (IBD) reaction, and the electron scattering (ES) reaction, The IBD reaction has no sensitivity to the SN direction but has a reaction rate in water that is two orders of magnitude larger than that of electron scattering. We adopt the IBD cross section of Strumia and Vissaini [47]. 4 On the other hand, electron scattering has a much smaller reaction rate but the strong correlation of the neutrino and lepton directions provides sensitivity to the SN direction. The cross section for electron scattering is obtained from field theory as where the light speed c is 1, G F = 1.166 × 10 −11 MeV 2 is the Fermi coupling constant, θ is the angle between the directions of the neutrino and the scattered electron and the constants A, B and C are summarized in Table 1 and depend upon the neutrino flavor [49]. Incidentally, we do not include neutrino interactions on oxygen [50], which have no direction sensitivity and smaller reaction rates than IBD. We set the detection threshold energy to be 5 MeV for both positrons and electrons.
ν: Table 1 Constant parameters for electron scattering. Here g V = −0.5 + sin θ W , where θ W is the Weinberg angle sin 2 θ W ≈ 0.23, and g A = −0.5 The IBD event rate is shown in Fig. 9 and the number of IBD events in several time windows is shown in Table 2 for the case of a supernova at 10 kpc. The total number of events in our model is of the same order as those of the previous studies [9,16,17,51]. As seen in Table 2, half of the events come in the first second, which is a region in which the cooling phase calculation is as important as the early phase calculation.
The IBD rates rapidly increase up to 2800 Hz for no oscillation, 2700 Hz for the normal hierarchy and 2600 Hz for the inverted hierarchy before gradually decreasing to 10 Hz at 20 seconds (Fig. 9). Three sharp peaks in the ES rate can be seen around the bounce, which is the only point at which the ES rate dominates the IBD rate for the no oscillations case. Supposing neutrino oscillations, these peaks are reduced from 1300 Hz to 300 Hz for the normal hierarchy and 600 Hz for the inverted hierarchy. The peaks correspond to the peaks of the luminosity (Fig. 5) and average energy (Figs. 6 and 7) of electron neutrinos. The ES rates fall to 100 Hz at around 20 milliseconds and gradually decrease to 0.3 Hz at 20 seconds. In the right panel of Fig. 9 there are no differences in either of the IBD or ES rates for neutrino oscillations. As a result, neutrino oscillations have little impact on the observation except for the peak rate of the ES events. 14/23  Table 2 Number of events at SK for a supernova at 10 kpc. For each reaction N tot is the total number of events, N (t min ≤ t ≤ t max ) is the number of events in the time interval between t min and t max , and the number in the brackets is the ratio relative to N tot .
We show a scatter plot of the event time and lepton energy in Fig. 10. As seen from this figure, the fraction of events with energy greater than 30 MeV declines with time as is expected from the decline in the average neutrino energy shown in Figs. 6 and 7.
The distribution of true event directions is shown in Fig. 11. Here the supernova is assumed to have occurred at the center of the plot and the horizontal (vertical) direction represents event latitude (longitude). Blue points are IBD events and the red points are electron scattering events. Note that IBD events are distributed throughout the plot, since the outgoing lepton does not preserve the direction of the incoming neutrino, while the electron scattering events concentrate at the center. As a result, the latter can be used to point back to the supernova. Figure 12 shows the time evolution of the mean energy of IBD events for no oscillation observed at each assumed supernova distance. In this figure the red curves show the theoretical expectation and the blue points show examples from the mock data sets. For the latter the width of the time bins is taken to be 1 second, except for the last two bins of the 50 kpc model. Here the error bars E err are defined as where N bin is the number of events, E i is the positron energy of the i-th event andĒ is the average energy of events in the time bin. The mock sample for the supernova at 1 kpc has mean energies that are almost the same as the theoretical curve. For the mock samples with from supernovae at 5 kpc and 10 kpc, the difference in the mean energies from the theoretical curve is less than 1 MeV before 5 seconds, but is larger after 5 seconds. However, the time evolution of the mean energies reproduces the theoretical curve well. For the mock sample at 50 kpc the time evolution does not track the expectation. Note that here we only consider only statistical uncertainties, which are symmetric, while the neutrino energy distribution in each time bin is asymmetric. In this sense the mock data distribution does not match the expectation well due to insufficient statistics. We will consider asymmetric errors including the effect of the shape of the distribution in the future. Figure 13 shows the charged particle energy and cos θ of individual events in the mock samples, where θ is the angle between the neutrino and the final state charged particle. The left panel shows the distributions for the model at 5 kpc, the middle is at 10 kpc, which 15/23 is the distance to the galactic center, and the right is at 50 kpc, which is the distance to the Large Magellanic Cloud. The top panel assumes no oscillations, the middle is the for the normal hierarchy and the bottom is for the inverted hierarchy. The electron scattering events (red) are shown stacked on the IBD events (blue).
From Fig. 13 we can see that the total event number decreases with the supernova distance, as expected. The IBD events are flat in the cos θ distribution. On the other hand, electron scattering events have peaks in the cos θ histograms toward the supernova direction for the mock samples for the supernovae at 5 kpc and 10 kpc. At 50 kpc, the same peak cannot be resolved clearly.

SK-Gd
The next stage of SK operations, known as SK-Gd [52], started in 2020. For the SK-Gd period a gadolinium compound has been dissolved into the pure water to improve the detector's ability to tag neutrons from IBD. SK-Gd mainly aims to detect the diffuse supernova neutrino background, however, it is also expected to improve the ability to determine the location of a supernova burst. Indeed, neutron tagging can be used to remove IBD events, leaving only electron scattering events for a precise determination of the supernova direction. Figure 14 shows the angular distributions of IBD and electron scattering events assuming no neutron tagging, 50%, and 90% tagging efficiency for a supernova at 10 kpc. It is clear that improved neutron tagging enhances the forward scattering peak.

Comparison with SN 1987A
To demonstrate the supernova analyzer in our framework, we compare our model with Kamiokande's observation of SN 1987A's neutrinos [5]. For this purpose we assume that the supernova distance is 51.4 kpc, the detector fiducial volume is 2.14 kton, and the detection threshold is 7.5 MeV in this subsection. We realize 100 Monte Carlo simulations using only IBD interactions. The left panel of Fig. 15 shows a histogram of the number of events obtained in these simulations. Though the average expected observation in our model is four events, 11 events were observed from SN 1987A. The right histogram shows the distribution of mean energies, which can be compared with 15.4 MeV observed at Kamiokande from Hirata et al. (1987) [5]. In comparison, the average energy of our simulation is 16.4 MeV.
As described in §3 our model corresponds to a weak explosion the number of events on the left side of Fig. 15 is within our expectations, though it is lower than the observations of SN 1987A. With this in mind, we expect the time evolution of the event accumulation to provide a comparison that is less dependent on the total number of events. Figure 16 shows a comparison of the cumulative number of events (CDF) over time. Though there are only 11 events from SN 1987A, their accumulation trend is similar to that of our simulation. To quantitatively compare with our model we perform a study based on the Kolmogorov-Smirnov test (KS test) in Fig. 17. Here the difference in the Kamiokande CDF and the model CDF at a time t from the first event is written as D(t). The maximum distance |D| max between the two is a measure of the compatibility of the two CDFs and is 0.238. that adopting different progenitor simulations may address our model's underestimation of the total number of events.

Summary and discussion
We have developed a framework to analyze neutrino events for studies of supernova physics and have shown predictions for neutrino emission over 20 seconds. This model follows neutrino emission from the beginning of collapse to the formation of a PNS by making extensions to the GR1D code. Using this simulation we describe a detector simulator to predict the neutrino signal at Super-Kamiokande. Mock data samples from our framework provide basic information to compare the simulation with the future observation of supernova burst. Assuming a supernova at 10 kpc about 2000 events are expected at SK, which is consistent with the model of Nakazato et al. (2013) [51]. In these simulations we find that electron scattering events from the neutronization burst can be distinguished from the more abundant IBD events, yielding information on the direction to the supernova. The gadolinium-loaded SK-Gd detector will improve this discrimination. Comparison of our model with the observation of SN 1987A indicates a consistent mean energy but our predicted event number of four is below the observed 11 events. Work is needed to increase the total flux of neutrinos without affecting the average neutrino energy. In the next work we will adopt different explosion models using the progenitor generation method of Ref. [53] to further explore the ability of long-time observations to constrain the physics of the burst.  Fig. 13 Distribution of cos θ and energy of charged particles from IBD (blue) and electron scattering (red) reactions. Electron scattering events (red) are shown stacked on IBD events (blue). The left, center and right panels are for the models with the supernova distances of 5 kpc, 10 kpc, and 50 kpc (bottom), respectively. The top, middle and bottom panels are for no oscillation, the normal hierarchy and the inverted hierarchy, respectively. The numbers in the legends indicate the total number of events.