The Neutrino Signal From Pair Instability Supernovae

TL;DR

This paper models the neutrino signals from pair-instability supernovae, highlighting their detectability in current and future neutrino detectors, which can help confirm the nature of these rare stellar explosions.

Contribution

It provides detailed predictions of neutrino fluxes from pair-instability supernovae across different progenitor masses and detector types, aiding future observational efforts.

Findings

Detectable neutrino signals from massive progenitors in current detectors at 10 kpc.

HyperKamiokande can detect signals from supernovae up to 50 kpc, including the Magellanic Clouds.

Light progenitors produce fewer neutrino events, often below detection thresholds.

Abstract

A very massive star with a carbon-oxygen core in the range of $64$ M$_{\odot}<M_{\mathrm{CO}}<133$ M$_{\odot}$ is expected to undergo a very different kind of explosion known as a pair instability supernova. Pair instability supernovae are candidates for superluminous supernovae due to the prodigious amounts of radioactive elements they create. While the basic mechanism for the explosion is understood, how a star reaches a state is not, thus observations of a nearby pair-instability supernova would allow us to test current models of stellar evolution at the extreme of stellar masses. Much will be sought within the electromagnetic radiation we detect from such a supernova but we should not forget that the neutrinos from a pair-instability supernova contain unique signatures of the event that unambiguously identify this type of explosion. We calculate the expected neutrino flux at Earth from two, one-dimensional pair-instability supernova simulations which bracket the mass range of stars which explode by this mechanism taking into account the full time and energy dependence of the neutrino emission and the flavor evolution through the outer layers of the star. We calculate the neutrino signals in five different detectors chosen to represent present or near future designs. We find the more massive progenitors explode as pair-instability supernova which can easily be detected in multiple different neutrino detectors at the `standard' supernova distance of $10\;{\rm kpc}$ producing several events in DUNE, JUNE and SuperKamiokande, while the lightest progenitors only produce a handful of events (if any) in the same detectors. The proposed HyperKamiokande detector would detect neutrinos from a large pair-instability supernova as far as $\sim 50\;{\rm kpc}$ allowing it to reach the Megallanic Clouds and the several very high mass stars known to exist there.

Figures

Tables

Table 1. Table 1: A summary of the detector types. ∗ See SNOwGLoBES documentation for discussion on phototube coverage. † For IceCube, the mass given is the ‘effective’ mass.

Detector	Type	Mass (kt)
Super-Kamiokande type: 30% phototube coverageFukuda et al. (2003)∗	Water Cherenkov	50
Hyper-Kamiokande typeHyp (2016)	Water Cherenkov	374
DUNE type detectorAcciarri et al. (2016)	Liquid Ar	40
JUNO type detectorAn et al. (2016)	Scintillator	20
IceCubeAbbasi et al. (2011)	Water Cherenkov	3500†

Table 2. Table 2: Numbers of interactions per detector for each mass ordering and a PISN at 10 kpc. These event counts are for the whole neutrino burst. The last two columns represents the number of interactions observed when neutrino oscillations are not taken into account.

Mass	Detector	NMO		IMO		Unoscillated
		Helm	SFHo	Helm	SFHo	Helm	SFHo
P150	Hyper-Kamiokande	$1.77$	$1.78$	$1.74$	$1.75$	$3.02$	$3.05$
	Super-Kamiokande	$0.24$	$0.24$	$0.23$	$0.23$	$0.40$	$0.41$
	DUNE	$0.14$	$0.14$	$0.15$	$0.15$	$0.25$	$0.25$
	JUNO	$0.10$	$0.10$	$0.10$	$0.10$	$0.17$	$0.17$
P250	Hyper-Kamiokande	$52.23$	$50.08$	$43.32$	$41.98$	$85.70$	$84.19$
	Super-Kamiokande	$6.98$	$6.69$	$5.79$	$5.61$	$11.46$	$11.26$
	DUNE	$2.95$	$2.78$	$3.17$	$3.06$	$5.30$	$5.20$
	JUNO	$3.13$	$3.00$	$2.48$	$2.40$	$5.06$	$4.97$

Equations

$$ı\frac{d{S}^{(XX)}}{dr}=H^{(X)}{S}^{(XX)},$$

$$U_V=\left(\begin{array}{ccc}1& 0& 0\\ 0& c_{23}& s_{23}\\ 0& -s_{23}& c_{23}\end{array}\right)\left(\begin{array}{ccc}{c}_{13}& 0& s_{13}{e}^{-i{\delta }_{CP}}\\ 0& 1& 0\\ -s_{13}{e}^{i{\delta }_{CP}}& 0& c_{13}\end{array}\right)\left(\begin{array}{ccc}{c}_{12}& s_{12}& 0\\ -s_{12}& c_{12}& 0\\ 0& 0& 1\end{array}\right),$$

$$H_V^{(f)}=\frac{1}{2E}{U}_V\left(\begin{array}{ccc}{m}_1^2& 0& 0\\ 0& m_2^2& 0\\ 0& 0& m_3^2\end{array}\right)U_V^{†},$$

$$\left(m_2^2-m_1^2,|m_3^2-m_2^2|,{\theta }_{12},{\theta }_{13},{\theta }_{23},{\delta }_{CP}\right)=\left(7.5\times 10^{-5}{\text{eV}}^2,2.32\times 10^{-3}{\text{eV}}^2,33.9^{\circ },9^{\circ },45^{\circ },0\right).$$

$${\rho }^{MSW}=\frac{m_N}{\sqrt{2}{G}_F{Y}_e}\left|\frac{\delta m^2\cos 2\theta }{2E}\right|$$

$$d{\mu }_{\alpha ,i,\theta }=\sqrt{2}{G}_{\text{F}}\frac{(1-\cos \Theta )}{4\pi c{\ell }^2}\left(\frac{l_{\alpha ,i}}{⟨E_{\alpha }⟩}-\frac{l_{\bar{\alpha },i}}{⟨E_{\bar{\alpha }}⟩}\right)2\pi \sin \theta r_i^2\delta r_{i}d\theta$$

$${\ell }^2=r_i^2+R^2-2r_{i}R\cos \theta$$

$$\cos \Theta$$

$$p_{\alpha \beta }=\sum _i|U_{V,\alpha i}{|}^2{P}_{i\beta }^{\text{(mf)}}\left(R_{\ast },R_0\right),$$

$$F_{\alpha }=\frac{1}{4\pi d^2}\sum _{\beta }{p}_{\alpha \beta }(E){\Phi }_{\beta }(E)$$

Full text

I Introduction

Pair instability supernovae (PISNe) are the explosive end point of very massive stars. After central carbon burning, the cores of these stars reach conditions that result in a dynamical instability due to the formation of electron-positron pairs (Barkat et al., 1967; Rakavy and Shaviv, 1967; Fraley, 1968). The conversion of photons to electron-positron pairs softens the equation of state and leads to a contraction which increases the density and temperature triggering explosive burning of the oxygen. The energy release from the burning is sufficient to reverse the contraction and unbind the star in the form of a PISN explosion (Heger and Woosley, 2002; Chatzopoulos and Wheeler, 2012). The existence of such supernovae would be important from both a chemical evolution Karlsson et al. (2008); Wise et al. (2012) and galaxy formation (Whalen et al., 2008) perspective. In addition to the large amounts of 56Ni they produce, the predicted nucleosynthesis pattern from PISNe has a strong odd-even staggering (Heger and Woosley, 2002). Such a pattern has not yet been found conclusively in any metal-poor stars Ji et al. (2015) but one tentative case has been identified Aoki et al. (2014).

For a star to experience the pair instability (PI), it needs to be massive enough to form a carbon-oxygen core in the range of $64$ M${}_{\odot}<M_{\mathrm{CO}}<133$ M*⊙* (Heger and Woosley, 2002). The corresponding zero age main sequence (ZAMS) mass is highly dependent on the details of stellar evolution and in particular on the mass loss history, which depends on the metallicity and rotation, but is not well understood (Vink, 2015). The long-standing expectation was that only non-rotating stars of zero metallicity (which are thought to have very little mass loss) born with ZAMS masses in the range of 140 M${}_{\odot}<M_{\mathrm{ZAMS}}<260$ M*⊙* explode this way (Heger and Woosley, 2002). Stars with higher ZAMS mass will collapse directly to a black hole; stars with an initial mass less than 140 M*⊙* may undergo a series of pair instability pulses (pulsational PISN) but ultimately explode as core-collapse supernovae. However it has been found that including rotation in stellar models can facilitate a chemically more homogeneous evolution. This leads to larger carbon-oxygen cores for the same initial mass, and hence shifts the mass range of those stars which undergo PISN to a lower ZAMS mass limit. In simulations at zero metallicity and with initial rotations rates of 80% of Keplerian velocity, stars with an initial masses as low as 40 M*⊙* are found to undergo a pulsational PISN and stars above 65 M*⊙* undergo a full PISN Chatzopoulos and Wheeler (2012). Also more recently, the viewpoint that the star had to be metal free to explode as a PISN has been challenged. Langer et al. (Langer et al., 2007) find that PISN may occur at metallicities as high as Z${}_{\odot}/3$ and simulations at solar metallicity find that a magnetic field at the surface of the massive star can reduce the mass-loss rate so that even stars at near-solar metallicity can explode as PISN Georgy et al. (2017). Thus it appears there is still a great deal we do not understand about exactly which stars explode as PISN. Observations of PISN would be of great value in answering the numerous questions about the progenitors and, if PISN can occur at non-zero metallicites, then it may be possible to observe a nearby event from which we could hope to obtain high quality information.

The theoretical expectation for the rate of PISN is very uncertain because of the aforementioned uncertain mass loss of massive stars. The estimate by Langer et al. Langer et al. (2007) is that there is one PISN event for every 1000 SNe in the local universe. To make this estimate the authors considered only the mass range 140 M${}_{\odot}<M_{\mathrm{ZAMS}}<260$ M*⊙* and thus do not include the previously mentioned extension down to lower ZAMS masses in the presence of rapid rotation. Additionally, while they include the effects of magnetic torques on mixing, they do not consider a surface dipolar magnetic field and its affect on the mass-loss rate, as in Georgy et al. (2017). Due to this uncertainty, we adopt a lower limit for the minimum ZAMS mass of PISN progenitors of 100 M*⊙. If that is the case, one can find many candidates of potential PISN progenitors within the Milky Way and nearby galaxies. Several stars with masses exceeding 100 M⊙* have been found in the Magellanic Clouds (Kudritzki R.-P., Lennon D.J., Haser S.M., Puls J., Pauldrach A.W.A., ). Crowthers et al. (Crowther et al., 2016) identified almost two dozen stars in R136/30 Dor whose initial masses exceed $\sim 100$ M*⊙*. Very massive stars are also known in NGC3603 (Schnurr et al., 2008) and the Arches cluster (Martins et al., 2008). Some of these stars exceed the previously assumed upper mass limit of $\sim 150\;{\rm M_{\odot}}$ for stars in metal-rich galaxies (Figer, 2005; Koen, 2006). If one of these stars were to explode as a PISN, one wonders how such an event might be identified.

There are several observational features which might be used to identify a PISN. In recent years, a number of authors have revisited the pair-instability mechanism computing expected lightcurves for various PISN models from zero metallicity to $10^{-3}\,$Z*⊙* Kasen et al. (2011); Whalen et al. (2013, 2014a); Kozyreva et al. (2017). These studies have found the peak luminosity in the models can be comparable to the peak luminosity of observed superluminous supernovae. Indeed a few of the observed superluminous supernovae with slower decay times have been discussed as possibly being of PISN origin: for example, the energy required to power the lightcurve of SN 2006gy (at a distance of $73\;{\rm Mpc}$) could be provided by the decay of 22 M*⊙* of 56Ni Smith et al. (2007); SN 2007bi has been suggested to be a pair-instability supernova by Gal-Yam et al. Gal-Yam et al. (2009), and another candidate was found by Cooke et al. Cooke et al. (2012). Superluminous supernova PS1-14bj also has a slowly evolving light curve, consistent with PISN models Lunnan et al. (2016). Of course, in addition to the energy from the decay of the large of amount of 56Ni, other processes, such as interactions of the ejecta with a circum-steller medium, could also contribute to the high luminosity. However we note that the calculated PISN lightcurves at non-zero metallicity are generally dimmer than the lightcurves at zero metallicity so PISN in the local Universe may become hidden among other classes of supernovae, see e.g. Whalen et al. (2014b) or Kozyreva et al. (2014). Should a nearby PISN occur, one would expect the increased quality of the observations might be able to find enough differences between a PISN and other supernova types to distinguish them but, as always, extracting information from the electromagnetic signal requires a high degree of modeling of the ejecta material which frequently introduces some uncertainty.

A much cleaner way to unambiguously distinguish PISN from other types of stellar explosions is the neutrino signal because the neutrino emission mechanism for PISN is the same as for Type Ia supernovae Kunugise and Iwamoto (2007); Odrzywolek and Plewa (2011); Seitenzahl et al. (2015); Wright et al. (2016, 2017) but on a much larger scale. The aim of this paper is to compute the expected neutrino signal from a PISN and determine its detectability in current and near-future neutrino detectors in order to answer the question of whether the neutrinos can be added to the suite of tools one can use to identify a PISN. We take into account the full time and energy dependence of the emission and the follow the neutrino flavor transformation in the envelope with a complete 3-flavor oscillation code.

Our paper is arranged as follows. In Sec. II we briefly describe the simulation of the pair-instability supernova used to compute the neutrino signal. In Sec. III we outline the algorithm for computing the neutrino spectrum and luminosity from the simulation and then in Sec. IV we present the calculation of the flavor transformation through the envelope. We then consider (Sec. V) a suite of neutrino detectors in order to determine how well, and to what distance, a pair-instability supernova could be detected. Section VI discusses the availability of nearby high mass stars and we conclude in Sec. VII.

II Supernova Simulation

For the PISN simulation we use two GENEC (Ekström et al., 2012; Yusof et al., 2013) progenitor models, P250 and P150 (Kozyreva et al., 2017; Gilmer et al., 2017). As their names suggest, the initial masses of these models are 250 M*⊙* and 150 M*⊙, respectively. Throughout the rest of our paper we shall distinguish our calculations by the progenitor model. The models have initial metallicity Z = 10-3* and were evolved without rotation. The Schwarzschild criterion was used for convection and the outer convective zone was treated according to the mixing length theory with $\alpha_{MLT}=1.0$ ($\alpha_{MLT}=l/H_{\rho}$, where $l$ is the mixing length and $H_{\rho}$ is the density scale height). Nearing the end of core carbon burning (after 25,000–30,00 time steps; 300–800 mass zones), when a sufficiently large part of the core has become unstable, the models are ready to explode. At this point the carbon-oxygen core mass of P250 is 126.7 M*⊙* which is near the upper end of the PISN mass range. The P150 model has a carbon-oxygen core mass of 65.7 M*⊙* which is at the lower end of the PISN mass range. Thus the neutrino emission we calculate from these two models should bracket what we can expect from an actual PISN.

To model the explosion phase, we used version 4.3 of the FLASH multiphysics code (Fryxell et al., 2000; Dubey et al., 2009). We used a subset of standard code modules suitable for the present application including the Helmholz EOS and the Aprox19 nuclear reaction network (Timmes and Swesty, 2000), the directionally unsplit hydrodynamics solver (Lee and Deane, 2009), the improved multipole solver for self-gravity (Couch et al., 2013), and its adaptive meshing capability. We run on a 1D grid of spherical geometry with an effective mesh resolution of 1.3 $\times$ 108 cm. The grid has the freedom to refine once (a local doubling of resolution) according to the density and temperature gradients. We set the refinement cutoff values to half of their default values. We consider only the innermost 5 $\times$ 1010 cm of the GENEC progenitor models because the outer material will not contribute to either the explosion dynamics, neutrino emission or neutrino flavor oscillation. The boundary conditions used are ‘reflect’ and ‘diode’ (a modification of ‘outflow’ in which matter is prevented from flowing into the grid) for the inner and outer boundaries, respectively. We then evolve the models through the contraction and explosion phases until nucleosynthesis is complete. Additional details on the progenitor models and the explosion calculation can be found in Gilmer et al. (2017). The thermodynamic trajectory of the core then serves as an input for the neutrino emission calculation.

III Neutrino Production

The strategy used to calculate the neutrino emission from our PISN simulations is similar to the one used in Wright et al. (2016) and Wright et al. (2017) for the case of Type Ia SN. We use the software package NuLib Fuller et al. (1982); Langanke and Martínez-Pinedo (2000); Langanke et al. (2003); Burrows et al. (2006); Oda et al. (1994); Steiner et al. (2013); O’Connor and Ott (2010); Sullivan et al. (2016) to calculate the emission as a function of time, neutrino energy and flavor and due to the different emission processes. The weak processes included are electron and positron capture on neutrons, protons and nuclei. The only thermal process included is neutrino pair creation due to electrons annihilating with positrons. This emission calculation is done by post-processing the simulations described in the previous section. Density, temperature, electron fraction and isotopic composition are extracted from each simulation for each time-slice and radial zone. These quantities are used to setup an Equation of State (EOS) which is then used by NuLib to compute neutrino emissivity. For comparison, two different EOSs have been considered. The first is the Helmholtz equation of state (based on Itoh et al. (1996)), the second is the SFHo EOS Steiner et al. (2013). The strength of the Helmholtz EOS is consistency, because it is the same EOS used in the original FLASH simulation. However, the simulation (with the Helmholtz EOS and the Aprox19 nuclear reaction network (Timmes and Swesty, 2000)) does not track all isotopes and thus neutrino emission from weak processes on isotopes not included, is missing. Nonetheless, given that the dominant isotopes are tracked, it is expected that these missing nuclei would not account for the bulk of neutrino emission from weak processes but it might miss spectral features especially towards higher neutrino emission energies. To test this hypothesis we consider a second EOS, SFHo, which assumes Nuclear Statistical Equilibrium (NSE) thus allowing us to a include emission from a wider range of nuclei. This EOS is an accurate description of the composition of material in a core-collapse supernova (CCSN) when $T>0.5\text{ MeV}\approx 5.8\text{ GK}$. However, our PISN simulation does not get as hot or dense as a CCSN, and thus our strategy, when using the SFHo EOS, is to only compute the weak emission from computational zones with $T>3\text{ GK}$ but when we calculate the thermal emission with this EOS we include zones with $T>1.2\text{ GK}$, which is approximately the lower bound of the SFHo EOS (see Steiner et al. (2013)). Our choice for the cutoff of $T>3\text{ GK}$ for zones from which we calculate the emission due to weak processes was so that we only include zones hot and dense enough to be close to NSE, while avoiding using a NSE EOS in zones that should not be in equilibrium. This strategy was proven to be successful when applied to the case of Type Ia SN Wright et al. (2016) and Wright et al. (2017).

Figure 1 shows the results of the neutrino emission calculations using NuLib. The time denoted as $t=0$ s is the time of maximum compression which is determined by locating the minimum in the gravitational energy. It is no surprise that the time of maximum compression is also the time of maximum neutrino emission. This is because increasing compression generally causes an increase in temperature which itself causes greater thermal activity, including neutrino emission. It is worth emphasizing that the neutrino signal from a PISN has a long duration, $\sim 30$ s, much longer than that of Type Ia supernovae Odrzywolek and Plewa (2011); Wright et al. (2016, 2017) and even core-collapse supernovae Totani et al. (1998); Fischer et al. (2010). In this figure, pair production (all six neutrino flavors), electron capture on nuclei ($\nu_{e}$), and electron and positron capture on nucleons ($\nu_{e}$ and $\bar{\nu}_{e}$) are given in red, blue and green respectively (positron capture on nuclei is very subordinate). The results for the P150 (left) and the P250 (right) models are presented as well as the results for the Helmholtz EOS and the SFHo EOS. Additionally, results when the Helmholtz EOS is assumed and a $T>3$ GK cutoff is imposed, are also displayed. For P250, there is no thin red line representing SFHo pair emission with $T>1.2$ GK because for this SN, a $T>3$ GK cut-off still captures the vast majority of the signal (unlike in the P150 case). In all cases we see that thermal (pair) emission dominates. For P150, pair emission using SFHo with a $T>1.2$ GK cut-off agree with the Helmholtz results with no cut-off. Also in agreement are the SFHo and Helmholtz pair emission results when both use a $T>3$ GK cut-off. For P250 at times during the peak of the neutrino emission, pair emission using SFHo with a $T>3$ GK cut-off agree with the Helmholtz results with no cut-off and those with a $T>3$ GK cut-off. Thus we conclude that EOS choice has little effect on neutrino pair emission, which is the dominant emission process. We can also conclude that the use of a temperature cut-off needs to be done carefully: using a cut-off temperature which is too high can lead to substantial error in the emission calculations especially for the for the cooler (lower mass) PISN.

For both electron capture on nuclei and electron and positron capture on nucleons, we again see that the temperature cut-off has a much small effect. This is clear because the thick blue and green lines match very well with their respective dashed lines, especially around times of peak emission for both the P150 simulation and the P250 simulation. However, the choice of EOS does make a large difference for this emission processes. We observe that the thin blue and green lines are above the thick or dashed blue and green lines, especially at early times. The difference indicates the SFHo EOS, which assumes NSE, over-predicts the amount of weak neutrino emission. Given that the P150 simulation barely reaches temperatures above the $T>3$ cut-off, the large over prediction for an EOS which assumes equilibrium has been reached is not surprising. However for both simulations, it is clear that the differences between the EOSs decrease as the explosions progress. If the mismatch at early times can be attributed to the SFHo predicting neutrino emission from nuclei that have not yet been synthesized then, later on when those nuclei have indeed been formed, the assumption of NSE is closer to reality and the SFHo EOS is more accurate. While these differences in the weak process emission are interesting, we stress that thermal emission is the dominant emission process in both simulations and because this emission can be calculated accurately, the mismatches in the weak emission only contribute to the uncertainty in sub-leading features of the spectrum.

That pair production is the dominant thermal process for a PISN is shown in Fig. 2 which shows the all-flavor neutrino luminosity for a number of thermal processes computed for the P250 simulation. The thermal processes included in the calculation are pair, photo-, plasma, bremsstrahlung and recombination neutrino processes and were calculated using the code sneut4 which is based on Itoh et al. (1996) which can be found at http://cococubed.asu.edu/code_pages/nuloss.shtml. These calculations have no spectral information but allow us to determine which processes are important. Figure 2 shows that pair production is by far the most important thermal process at all times. Also, the pair results from NuLib very nearly overlap with the pair results from sneut4, giving further confidence to our results. A similar analysis was done for the P150 simulation and the conclusion that pair production is by far the most important thermal process remains true.

The advantage the NuLib calculations have is that they contain spectral information. Figure 3 shows the neutrino spectrum as a function of neutrino energy for all six flavors and for the time slice around maximum emission. The red (blue) lines represent the P250 (P150) simulation and the solid and dashed lines represent the results using the Helmholtz and SFHo EOSs respectively. The figure shows, that for all flavors other than $\nu_{e}$, the spectra have the same expected thermal shape with a peak at $\sim 1.5$ MeV for the P250 simulation and $\sim 1$ MeV for the P150 simulation. The $\nu_{e}$ spectra has a strong contribution from weak emission which boosts the overall rate and additionally introduces a spectral feature at 10 MeV when we use the SFHo EOS. Given the results from Wright et al. (2016) and Wright et al. (2017), we recognize this 10 MeV peak as weak emission mostly from electron capture on copper. We find that this 10 MeV peak is present at all time slices (for the SFHo results). However, given the aforementioned mismatch between the weak process emission using the two EOSs, our confidence in the presence of this feature in the spectrum is low initially and increases as the simulation progresses (especially for the P250 simulation).

Finally, the time evolution of these spectra are simply described as brightening then fading with little change of the shape except for the spectral feature around 10 MeV in $\nu_{e}$ which becomes relatively more pronounced later in the simulation. This behavior is demonstrated by the gray region which is the spectra from the SFHo, P250 results at $t=12.6$ s. Another interesting point is that the thermal peak energy does not change much as the supernova progresses. For the Type Ia supernova analysis previously performed with this same strategy, the thermal peak shifted to lower energies as the explosion progressed implying that each zone in the Type Ia supernova cooled over time. The different evolution for the PISN indicates that the hot zones remain hot for a much longer period of time than in the Type Ia simulations.

IV Neutrino Flavor Transformation

In order to accurately determine the neutrino signal from a PISN that reaches Earth, the effects of neutrino oscillation need to be accounted for. This includes both the neutrino flavor oscillations that take place during neutrino propagation through the stellar mantle and the decoherence that arises as the neutrino propagates through the vacuum between the supernova and Earth. The calculation strategy we have used for the flavor transformations is the same as the one used in Wright et al. (2016) and Wright et al. (2017).

Neutrino oscillation phenomena are calculated by solving the Schrödinger equation. The Schrödinger equation in some basis $(X)$ is

[TABLE]

where $H^{(X)}$ is the Hamiltonian in that basis and the evolution matrix $S^{(YX)}(r_{2},r_{1})$ connects the neutrino states in the $(X)$ basis at some initial position $r_{1}$ to the states in a possibly different basis $(Y)$ at $r_{2}$. Note that when we are not referring to a specific basis, we shall drop the superscript. The evolution of the antineutrinos is governed by an evolution matrix $\bar{S}$ which evolves according to the same equation but with a Hamiltonian $\bar{H}$. From the evolution matrix we define the transition probability $P^{(YX)}_{yx}=\left|S^{(YX)}_{yx}\right|^{2}$ as the probability for some state $x$ in the $(X)$ basis at $r_{1}$ to be in state $y$ in the $(Y)$ basis at $r_{2}$. Antineutrino transition probabilities will be denoted with an overbar. The two bases commonly used in the literature are the flavor basis – which has basis states $\nu_{e}$, $\nu_{\mu}$, and $\nu_{\tau}$ – and the matter basis (often called the mass basis when in vacuum) – which has basis states $\nu_{1}$, $\nu_{2}$, and $\nu_{3}$ Bethe (1986). For future reference, we shall use Greek letters $\alpha$ and $\beta$ to denote generic flavor states and the Roman symbols $i$ and $j$ to denote generic matter basis states. Again, we denote the antineutrino states in the two bases by an overbar. Using this notation, the neutrino transition probability for a neutrino that was initially state $j$ in the matter basis to be detected as state $i$ in the matter basis is denoted by $P^{\text{(mm)}}_{ij}=P(\nu_{j}\rightarrow\nu_{i})$.

The mixing matrix $U_{V}$ which describes the transformation of the flavor basis to the mass basis is parameterized by three angles $\theta_{12}$, $\theta_{23}$ & $\theta_{13}$, and a CP-violating phase $\delta_{CP}$, and can be written as

[TABLE]

with $c_{ij}=\cos\theta_{ij}$ and $s_{ij}=\sin\theta_{ij}$.

IV.1 Vacuum Hamiltonian

In order to solve the Schrödinger equation we need to define the Hamiltonian $H$ ($\bar{H}$ for the antineutrinos). In the vacuum, $H$ is given by a single matrix $H_{V}$ whose exact structure depends upon the basis one is using. In the mass basis the three neutrino states have masses $m_{1}$, $m_{2}$ and $m_{3}$ and the vacuum Hamiltonian is diagonal. The vacuum Hamiltonian in the flavor basis, denoted by $H_{V}^{(f)}$, is

[TABLE]

with $E$ the neutrino energy. The mixing matrix angles and phases, together with the squared differences between neutrino masses, are part of what determines the nature of the oscillation phenomenology. In this paper we adopt the following values for all results

[TABLE]

where $m^{2}_{3}-m^{2}_{2}>0$ is for normal mass ordering (NMO) and $m^{2}_{3}-m^{2}_{2}<0$ is for inverted mass ordering (IMO). For the antineutrinos the vacuum Hamiltonian $\bar{H}_{V}$ is simply the complex conjugate of $H_{V}$.

IV.2 Matter Hamiltonian

In matter we must add to $H_{V}$ an additional term often known as the ‘matter Hamiltonian’, $H_{M}$, so that $H=H_{V}+H_{M}$. In the flavor basis $H_{M}^{(f)}$ is given by $H^{(f)}_{M}=\sqrt{2}\,G_{\text{F}}\,n_{e}\,\text{Diag}\left(1,0,0\right)$, where $n_{e}$ is the neutrino density and $\text{Diag}\left(1,0,0\right)$ is a 3-by-3 matrix where the only non-zero entry is unity in the first diagonal element. For the antineutrinos the matter Hamiltonian $\bar{H}_{M}$ is related to $H_{M}$ by $\bar{H}_{M}=-H_{M}$.

Thus we see neutrino flavor transformation depends on the density of the material and its electron fraction $Y_{e}$. For the PISN simulations considered here, the electron fraction is very nearly $Y_{e}=0.5$ for the entire duration of significant neutrino emission. The stellar density profile is shown in Fig. 4. The time slices chosen are the first, peak, last and the first time slice with a shock present. We see that the star develops a shock at $\sim 7\times 10^{9}$ cm at $t\approx-5.4$ s and $\sim 10^{4}$ g/ccm for the P250 simulation and similarly at $\sim 10^{10}$ cm at $t\approx 6.8$ s and $\sim 10^{2}$ g/ccm for the P150 simulation. This shock then propagates outward into regions of lower density. The densities at which MSW resonances occur for 1 MeV neutrinos are indicated by the dashed lines. To calculate the density of the high-density MSW resonance, $\rho^{MSW}_{High}$, we use the two-flavor approximation

[TABLE]

and $\delta m_{32}^{2}$ for the mass splitting $\delta m^{2}$ in the formula and $\theta_{13}$ for the mixing angle $\theta$, while for the low-density MSW resonance $\rho^{MSW}_{Low}$, we use $\delta m_{12}^{2}$ for the mass splitting and $\theta_{12}$ for the mixing angle. In Eq. 4 $m_{N}$ is the nucleon mass.

IV.3 Neutrino self interaction potential

Another possible contribution to neutrino flavor oscillation arises due to neutrino self-interaction, i.e. we must add a third term to $H$ (see Duan and Kneller (2009); Duan et al. (2010); Mirizzi et al. (2015) for reviews where neutrino self interaction is discussed in core-collapse supernovae). Naively, one might not anticipate this contribution to be important in a PISN because the neutrino densities are quite low. However, the matter density - shown in Fig. 4 - is also much lower than one finds in core-collapse supernovae and thus it is not immediately obvious whether or not neutrino self-interactions can be ignored.

In order to check whether self-interaction effects might occur, we compute the strength of the neutrino self-interaction upon a radially emitted test neutrino. One important detail we need to include in our estimate is that the matter in a PISN is nowhere optically thick to neutrinos. Given that our simulation is one dimensional, our first step is to map the density and neutrino luminosity onto a sequence of spherical shells labeled by an integer $i$ and with a thickness $\delta r_{i}$ given by the radial grid spacing. For shell $i$, the neutrino luminosity emitted in flavor $\alpha$ per unit volume is $l_{\alpha,i}=L_{\alpha,i}/V_{i}$ where $V_{i}$ is the shell volume and $L_{\alpha,i}$ the luminosity of the shell in neutrino flavor $\alpha$. Because the matter densities are much too low to trap neutrinos, the neutrinos are emitted isotropically from each volume element of a shell. Thus, for a test neutrino emitted radially at a distance $R$ from the center of the star, the contribution to the neutrino self interaction potential (not including oscillations) from neutrinos emitted from a volume element which is at a distance $r_{i}$ from the center of the star and at an angle $\theta$ with respect to the ray of the test neutrino, is

[TABLE]

where $G_{\text{F}}$ is the Fermi constant, $\Theta$ is the angle between the ray of the test neutrino and the neutrinos emitted from the volume element, $\ell$ is the distance between the volume element and the test neutrino, $\left\langle E_{\alpha}\right\rangle$ is the average neutrino energy and $\bar{\alpha}$ represents the antineutrino of flavor $\alpha$. The distance $\ell$ is simply

[TABLE]

and the angle $\Theta$ is

[TABLE]

A particular shell’s contribution to the total neutrino luminosity is assumed to scale as $L_{\alpha,i}\propto T_{i}^{9}$ in accordance with the results derived in Burrows, Reddy, & Thompson Burrows et al. (2006) for pair production neutrinos. We determine the proportionality constant by constraining the sum over shells to give the total luminosity we have previously calculated. The contribution from an entire shell can be obtained by performing the integration over the angle $\theta$ and the contribution from the whole star can be obtained by summing over all shells. In this way we calculate the total strength $\mu_{\alpha}$ of the neutrino self-interaction.

Our temperature cut-off strategy when computing the emission means that if a shell has a temperature below the cut-ff then its neutrino luminosity is zero. Thus, for each time-slice, only shells within a certain radius emit neutrinos. This outermost neutrino-emitting shell is not a neutrino-sphere because the neutrinos are not all emitted from its surface (which is the case for a neutrino-sphere in CCSN). Nonetheless, it is instructive to treat the outermost shell as a neutrino-sphere and compute the self-interaction potential via the usual equations in the neutrino BULB model (i.e. via Eqn. 40c taken from Duan et al. (2006)) to compare to our shell calculation. These pure BULB model results are an upper limit for our shell model calculations because, by assuming all neutrino luminosity comes from the outermost emitting shell, there is greater flux of neutrinos from large angles $\Theta$ and thus the self-interaction is artificially increased via the factor $\left(1-\cos\Theta\right)$ in Eqn. 5.

Figure 5 displays the self-interaction potential strengths as a function of radius for the P250 simulation. The six subplots represent time-slices around the peak neutrino emission and each subplot displays the $\left(L_{\alpha},\left\langle E_{\alpha}\right\rangle\right)$ data valid for that time-slice as an inset table. The black line represents the matter potential which is $\sqrt{2}\,G_{\text{F}}\,Y_{e}\,\rho/m_{N}$, the blue and red lines are the electron and muon-or-tau self-interaction potentials respectively. The dashed blue and red lines represent the corresponding self-interaction potentials calculated via the BULB model. The gray horizontal lines represent the vacuum potentials with $H_{\text{V,ij}}=\Delta m_{ij}^{2}/(2\left\langle E\right\rangle)$ and the average energy being the mean across all flavors. From Fig. 5 we conclude the self-interaction potentials are many orders of magnitude below the matter potential or vacuum Hamiltonian and indicates we need not consider neutrino self-interaction effects when calculating the the neutrino flavor evolution. As expected, the self-interaction potentials calculated via our shell model are lower than those calculated by the BULB model. Additionally, Fig. 5 demonstrates how the shock affects the high density MSW resonance before affecting the low density MSW resonance. A similar analysis was performed for the P150 simulation and the neutrino self-interactions potential was even weaker than in the P250 simulation.

IV.4 Matter basis transformation probabilities

The matter basis is defined to be the basis where the eigenvalues of the total Hamiltonian $H=H_{V}+H_{M}$ appear on the diagonal. The eigenvalues are arranged so that they reflect the ordering of the masses appearing in the vacuum Hamiltonian. The matter basis for the antineutrinos is defined similarly. The significant advantage to using the matter basis over others is that for adiabatic evolution the evolution matrix in this basis is close to diagonal and transition probabilities $P_{ij}$ are constant.

Indeed we find adiabatic evolution is almost exactly the case at early times because the density profiles of the simulations are not steeply falling functions of the radius $r$. This means the matter basis neutrino transition probabilities between states $\nu_{i}$ and $\nu_{j}$ are approximately unity for $i=j$ and zero for $i\neq j$. The same is also true for the antineutrinos. However, adiabatic evolution cannot continue for all epochs because the adiabaticity depends upon the density gradient and the density profiles from our simulation have moving shocks. The presence of shocks changes the adiabaticity and the change is most noticeable at those epochs where the shock is in the vicinity of the high (H) and/or low (L) Mikheyev-Smirnov-Wolfenstein (MSW) resonances Dighe and Smirnov (2000). At these times there is near-total ‘hopping’ from one matter eigenstate to the other Parke (1986) which translates into the neutrino remaining in whichever flavor state it was in before the shock. Thus we expect significant changes in the flavor transition probabilities oscillations as the shock moves through the mantle of the supernova which will, in turn, lead to a change in the number of events we detect here on Earth. For a normal mass ordering both the H and L resonances affect the neutrinos because the L resonance is in the 12 mixing channel (by ij mixing channel we mean mixing is between matter states $\nu_{i}$ and $\nu_{j}$), and the H resonance is in the 23 channel. For an inverted mass ordering the L resonance remains in the 12 neutrino mixing channel but the H resonance moves to the 13 antineutrino mixing channel Dighe and Smirnov (2000).

Figure 6 shows the neutrino matter basis transition probabilities $P^{\text{(mm)}}_{ij}\left(E_{\nu}\right)$ for both mass orderings and as a function of energy for a number of time-slices around peak neutrino emission for the P250 simulation. For the top panel, which represents normal mass ordering, we see energy-dependent departures from adiabatic behavior in the 12 and 23 mixing channels. This partial diabatic evolution appears at $t\sim-5.4$ s for neutrinos with energy around 0.2 MeV in the 12 mixing channel and around 3 MeV in the 23 mixing channel, and moves to higher energies at later times. The evolution of the transition probabilities are easily understood. As the SN evolves, the shock moves from high to low densities and the MSW resonances are inversely proportional to the neutrino energy. As the figures shows, the diabaticity in the 23 mixing channel ($\sim 80\%$) is much higher that in the 12 mixing channel ($\sim 20\%$) due to the smaller value of $\theta_{13}$ compared to $\theta_{12}$.

The bottom subplot of Fig. 6 represents neutrino matter basis transition probabilities for inverted mass ordering. Here we only expect diabatic behavior in the 12 mixing channel because only the L resonance occurs for neutrinos in inverted mass ordering. Indeed that is what is is visible in the bottom subplot and the energy and time dependence of this Low MSW resonance is the same as for the the top subplot.

Similarly, Fig. 7 represents the antineutrino matter basis transition probabilities $\bar{P}^{\text{(mm)}}_{ij}\left(E_{\bar{\nu}}\right)$ as a function of energy for a number of time-slices around peak antineutrino emission for both mass orderings for the P250 simulation. Figure 7 shows similar behavior as seen in Fig. 6. For antineutrinos in the inverted mass ordering, the diabatic behavior is seen in the 13 mixing channel as was seen in the 23 mixing channel for neutrinos in the normal mass ordering. This happens because the H resonance switches to the antineutrinos for the inverted mass ordering and is between antineutrino matter states $\bar{\nu}_{1}$ and $\bar{\nu}_{3}$. A small amount of mixing is seen in the 12 antineutrino mixing channel in both the normal and inverted mass ordering but with much smaller amplitude. For the energies shown, the scale of the diabatic behavior for antineutrinos (in either mass ordering) in the 12 mixing channel is small ($<10\%$).

IV.5 Flavor basis transition probabilities at Earth

Once the neutrinos have propagated through the star and start their long voyage to Earth, we need to account for the decoherence of the neutrino wavepacket Giunti et al. (1998). Accounting for this effect, the probability $p_{\alpha\beta}$ that a neutrino emitted as flavor $\beta$ in the supernova will be detected as flavor $\alpha$ at Earth is given by

[TABLE]

where $R_{0}$ represents the radius of the neutrino production point (near the center of the supernova which we take to be the origin) and $R_{*}$ represents the radius of the outer edge of the supernova. One defines antineutrino flavor basis survival probabilities $\bar{p}_{\alpha\beta}$ by using the matrix $\bar{P}^{\text{(mf)}}_{i\beta}\left(R_{*},R_{0}\right)$. Note these probabilities $p$ and $\bar{p}$ are different from the transition probabilities we discussed in the previous section because they do not come from an evolution matrix. In Fig. 8 the electron neutrino and antineutrino survival probabilities $p_{ee}$ and $\bar{p}_{ee}$ as a function of energy. Neutrino (antineutrino) survival probabilities are displayed in the top (bottom) row for each PISN simulation and the left (right) column represents the normal (inverted) mass ordering. It is clear the general trend is that flavor oscillations hamper $\nu_{e}$ or $\bar{\nu}_{e}$ survival. Only for $\bar{\nu}_{e}$ in NMO does greater that 50% of the electron flavor survive. For the P250 simulation, the ‘humps’ which appear for electron neutrinos in the NMO and electron antineutrino in the IMO are mostly caused by the passage of the shock through the H resonance. Without the shock the probability that an electron neutrino or electron antineutrino would be detected in the same state at Earth is just a few percent. As seen in the electron neutrinos in the IMO, the passage of the shock through the L resonance produces a much smaller effect and at later times. For the P150 simulation, it is clear that there are very few energy / time dependent oscillation effects. In this case the neutrino flavor evolution is almost entirely adiabatic at all times because the shock forms quite late and at low densities. In fact, for the P150 simulation, the shock never passed through the H resonance because it formed beyond that radius, and moved so slowly that it barely passes through the L resonance before the neutrino emission subsides.

IV.6 The neutrino flux at Earth

By combining the flavor basis oscillation probabilities and the neutrino emission spectra, we calculate the neutrino flux seen by neutrino detectors at Earth. These fluxes are given by

[TABLE]

where $d$ is the distance to the supernova and $\Phi_{\beta}(E)$ is the differential spectrum of flavor $\beta$ at the point of emission. In what follows we set the distance to the PISN to be $d=10\;{\rm kpc}$ but will comment about the event rates at other distances. Figure 9 shows the flux as a function of neutrino energy including oscillations (at peak emission time). By comparing Fig. 9 and Fig. 3 the effects of neutrino oscillation can be easily seen. Once again, the red (blue) lines represent the P250 (P150) simulation and the solid and dashed lines represent the results using the Helmholtz and SFHo EOSs respectively. The gray region is the oscillated spectra from the SFHo, P250 results at $t=12.6$ s. In the SFHo results, the 10 MeV spectral feature mentioned in Sec. III, which prior to oscillation was only present in the $\nu_{e}$ flux, is now present in all three neutrino flavors. Furthermore, in all cases, even though the unoscillated electron flavor flux dominates over the other flavors, the oscillated flux spectrum shows that much of this electron flavor has oscillated into muon and tau flavor in both the neutrino and the antineutrino case. The effects of the shock are difficult to see on the logarithmic scales used to make these figures, but they exist between $t=-5.4$ s and $t=0$ s for the P250 simulation. At these times the changes, due to diabatic evolution induced by the shock, occur at energies sufficiently close to the energy of peak emission to make an appreciable difference. All other diabatic effects of the shock occur at energies where the luminosity was low enough that any variation is unnoticeable.

V Neutrino Detection

Ultimately, the goal is to be able to measure a neutrino signal in a neutrino detector and by using that signal, be able to identify source characteristics. In order to determine what the neutrino signal from our PISN simulations would be, we use the software package SNOwGLoBES to simulate the expected event rates in a variety of detectors. Table 1 lists the detectors considered together with their type and mass. The detector description includes the work "type" to indicate that the detector model used in SNOwGLoBES is only approximate. This is because, for the existing detectors, more accurate descriptions of the detectors are not openly published and for the future detectors, the detector specifications are not finalized. Furthermore, realistic detector thresholds are not included here because they are not well established for the future detectors and thus we do not include them for the existing detectors for consistency. All results include a computational 0.5 MeV threshold. We also assume perfect detection efficiency, that is, the incoming neutrino’s energy is perfectly reconstructed from the detected particle’s energy. This assumption is again made because efficiencies are not known for all the detectors considered and thus the events calculated here are labeled interaction events. These assumptions are overly optimistic, but they allow more consistent comparison between detectors given the information available. At the end of this section, an analysis of the effect of detector thresholds and efficiency (also called smearing) is presented.

Table 2 shows the expected interaction event count for all the detectors listed in Table 1, except IceCube, for a PISN occurring at 10 kpc. The results are reported for both the P150 and P250 simulations, both EOSs and both mass orderings as well as the case of no neutrino oscillations for comparison. These interaction event counts are the total from the whole $\sim 30$ s neutrino burst as well as across all energies. For reference, the expected event rates in these same detectors for a Type Ia supernova at 10 kpc are roughly two orders of magnitude smaller Odrzywolek and Plewa (2011); Wright et al. (2016, 2017) than the P250 case, while the event rates for a CCSN are approximately three orders of magnitude larger Mirizzi et al. (2015). Firstly, we note that the effect of neutrino oscillations is to lower the event count in all detectors and all cases. This is because charged current interactions usually produce the strongest signals and the overall trend of oscillations is to convert electron flavor into muon and tau flavor which are detected via neutral current. Next, we note that for Super-Kamiokande (SK), Hyper-Kamiokande (HK), and JUNO there are more events for NMO than IMO whereas for DUNE, there are more events for IMO than for NMO. This is because SK, HK and JUNO are sensitive to the inverse beta decay (IBD) channel and DUNE is not. As shown in Figure 9, the $\bar{\nu}_{e}$ flux is about an order of magnitude greater in NMO than in IMO. Our investigations (detailed below) reveal that it is the IBD channel that is responsible for why the NMO event count is greater than the IMO event count in the detectors that are sensitive to it. Next, Table 2 shows that SK, DUNE and JUNO would detect several and HK would detect several tens of interaction events over the span of the neutrino signal for a PISN at the upper end of the progenitor mass range. Such a detection would indicate high-temperature nuclear burning for an extended period of time. This allows us to easily distinguish a PISN from other types of supernovae: the neutrino signal of a CCSN reaches a peak luminosity during the neutronization burst which occurs $\sim 100\;{\rm ms}$ after core bounce Fischer et al. (2010), the entire duration of the neutrino emission from a Type Ia is $\sim 5\;{\rm s}$ at most Odrzywolek and Plewa (2011); Wright et al. (2016, 2017), while the PISN neutrino signal lats for $\sim 30$ seconds and takes $\sim 15$ seconds to reach peak luminosity. Even if other discriminators are not available, this temporal evolution alone would be a smoking-gun signature of a PISN and that some stars must end their lives in this way. The table also reveals that the differences between the predicted number of events by the two EOSs is minor. Lastly, it is clear that the P150 simulation would just barely be observable at $d=10\;{\rm kpc}$ and that the P150 simulation produces around 25 times fewer events than the P250 simulation.

By going beyond total event counts and considering the time and energy structure of the signal we can investigate the observability of the various temporal and spectral features in the neutrino signal. Figure 10 displays the interaction event rates expected to be observed in the detectors under consideration from a PISN at 10 kpc. The top four panels show that, in log scale, NMO and IMO give roughly the same event time-profile at all detectors however, the bottom four panels, using a linear scale, reveal the aforementioned differences between the rates from NMO and IMO. For both orderings our general expectation is for a Gaussian-like burst of neutrinos over a period of $\sim 30$ s. Solid and dashed lines represent results from the Helmholtz and SFHo EOS respectively and it is clear that the choice of EOS has little impact. If we look very carefully in the P250 simulation we find in the IMO case a slight increase between $t=-5$ s and $t=0$ s due to shock induced diabatic evolution. This burst structure shows that at the peak, HK would be seeing a little more than 5 interaction events per second for the NMO. However, for SK, DUNE and JUNO, bigger time bins (or a nearer SN) would be required to confidently see activity per bin. As expected, the P150 simulation has a much wider and dimmer signal.

Figure 11 shows the energy structure of the neutrino burst released from a 10 kpc PISN as seen in the detectors under consideration. From all 8 plots, it is clear that the majority of the signal is below 5 MeV in all cases. The top plots reveals that even though the 10 MeV spectral feature from electron capture on copper is seen in both mass orderings using the SFHo EOS, it is a very small part of the signal and is thus less likely to be detected. The bottom 4 plots show a slight bump that is visible around 2-3 MeV in the NMO differential spectrum. This bump is from the contributions of the IBD channel which are more significant in NMO than in IMO. This is a very interesting spectral feature that, if ever detected at high enough statistical significance, could place constraints on the neutrino mass ordering.

We now examine the spectral event rate per channel in order to quantify and illustrate the previous points about the strength of the IBD in NMO over IMO. Figure 12 shows the interaction event count per 0.5 MeV energy bins for the full $\sim 30$ s duration of a P250 PISN at 10 kpc (using the SFHo EOS). The reason that the IBD channel is present for HK but not for DUNE is because it is dominant for detectors built from materials that are composed of hydrogen e.g. water and scintillator. As previously noted, Fig. 9 reveals that the $\bar{\nu}_{e}$ flux is about an order of magnitude greater in NMO than IMO. This is the reason why the IBD contribution is shown by Fig. 12 to be much larger in NMO than in IMO. This gives rise to the strong double peak seen in Fig. 11 for the water and scintillator detectors in NMO. Figure 12 also shows that the neutrino elastic scattering on free electrons form a significant contribution to the total event count. For HK, contributions from interactions with oxygen are quite minimal but for DUNE, contributions from interactions with argon are important, especially at energies above 4 MeV. Below $\sim 8$ MeV, where the bulk of the signal is, Fig. 12 is unchanged by choice of EOS. Above $\sim 8$ MeV, Helmholtz spectra continue to decline according to the thermal trends visible below $\sim 8$ MeV.

Finally, we examine the effects of detector thresholds and energy smearing which, considering that the bulk of the signal is below 5 MeV, are expected to be an important consideration. Figure 13 shows what percentage of events would survive given a particular detector event threshold. This analysis was for interaction events (left column), for energy-smeared events (right column), and for both mass orderings. The solid and dashed lines are for the P250 and P150 simulation respectively. The figure reveals that a threshold of 5 MeV would reduce the detectable signal in all cases essentially to zero. However, a threshold of 2 MeV means that, for the P250 simulation, Hyper-K, Super-K and JUNO would retain more than 50% of their interaction events while DUNE would retain closer to 40%. The case is worsened when energy-smearing is added but for a threshold of 2 MeV we see Hyper-K and Super-K still retain more than 50% of their events, JUNO drops to 30-40% (dependent on mass ordering) and DUNE drops to a bit more than 10% retention. In short, energy smearing and detector thresholds do reduce the detectable signal for the P250 simulation, but given a threshold of around 2 MeV, much of the signal can still be detected, especially in the case of Hyper-K and Super-K. We also note that the spectral peak between 2-3 MeV which appeared for the NMO would still be visible for a 2 MeV threshold and, if observed, would still suggest the mass ordering is normal. The case of the P150 simulation however, shows that the Hyper-K and Super-K would still see 50% of their smeared events, while for DUNE and JUNO, the signal is significantly attenuated.

V.1 IceCube Event Rate

The PISN neutrino signal in IceCube must be treated differently than the other detectors listed in Table 1. The IceCube neutrino detector is located inside the ice that covers Antarctica and, while it is not primarily designed to investigate low energy neutrinos, it nonetheless is sensitive to them. These ‘detections’ will consist of an increase of the background hum of low-energy events in the individual Photo-Multiplier Tubes. The events IceCube measures from a PISN are mostly from IBD but there is some contribution from elastic scattering off of electrons in the ice. However no direction or energy information of individual neutrino interactions will be extracted so the two event types cannot be distinguished. What can be determined is the overall flux of the neutrinos and the time structure. In order to be detected the PISN neutrino signal needs to be sufficiently above statistical fluctuations of the usual background event rate. Thus, even though the PISN will produce many events in IceCube, the metric of detectability is different from the other detector types considered.

Figure 14 shows the detection confidence of a P250 PISN at IceCube for a given PISN distance. As expected, given the abundance of $\bar{\nu}_{e}$ and hence IBD events, detection is easier in NMO than in IMO. However, in either mass ordering, it is clear that a PISN would have to be within one kiloparsec in order to be confidently detected above background at IceCube. The existence of stars with mass high enough to explode as a PISN at these relatively nearby distances seems quite unlikely. The P150 PISN case is even less likely to be detected by IceCube (3$\sigma$ requires the PISN to be within $\sim 70$ pc).

VI Astronomical context

Thus far, we assumed that the PISN was located at 10 kpc. This was chosen because it is roughly the distance to the Galactic center and the distance used in the SN community as the most probable distance for the next Galactic SN. If we adopt a conservative lower limit for the ZAMS mass of a PISN of $100$ M*⊙* then this proximal distance is not unrealistically close: within the Milky Way there are several stars with masses around $100$ M*⊙* or above. These include:

several stars in the Arches cluster Nagata et al. (1995); Cotera et al. (1996), close to the Galactic center at a distance of $d\approx 8\;{\rm kpc}$ Reid (1993); Malkin (2012); Camarillo et al. (2017), are estimated to have masses above $100$ M*⊙* Gräfener et al. (2011), 2. 2.

the OB association Cygnus OB2 Massey and Thompson (1991) at a distance of $d=1.7\;{\rm kpc}$ Schuster et al. (2009) also contains stars with masses approaching $100$ M*⊙*, 3. 3.

the mass of the primary in the HD 15558 system, at a distance of $d=2.3\;{\rm kpc}$, was estimated to be $150\pm 50$ M*⊙* De Becker et al. (2006) (see also Garmany and Massey (1981)), 4. 4.

the primary star in the binary $\eta$ Carina which is, coincidentally, also at a distance of $d=2.3\;{\rm kpc}$, has a mass that is also modeled as being greater than $100$ M*⊙* Kashi and Soker (2010); Clementel et al. (2015).

Whether these Galactic very high mass stars are too metal rich to explode as a PISN will have to be determined when the event occurs. The neutrino signal will easily distinguish the explosion type. A number of very high mass stars are also known within the Large Magellanic Cloud (LMC) - which has a lower metallicity of $0.43\,$Z*⊙* 111Computed from $\log(Z/Z_{\odot})=[\mathrm{Fe}/\mathrm{H}]$ and the values of [Fe/H] found in Choudhury et al. (2016) - and which is at a distance of $d=49.97\;{\rm kpc}$ Pietrzyński et al. (2013):

The R136 open cluster contains nine stars with masses greater than $100$ M*⊙* at 1-sigma Crowther et al. (2016). 2. 2.

NGC 3603-A1 with the more massive component of the binary having a mass $116\pm 31$ M*⊙* Schnurr et al. (2008). 3. 3.

WR21a whose more massive component has a mass of $103.6\pm 10.2$ M*⊙* Tramper et al. (2016).

For stars at this distance the event rates given in table 2 need to be scaled downward by a factor of $25$. Thus we see Hyper-Kamiokande is the only detector at the present time or in the near future capable of detecting events from a large PISN at LMC distances. Again, the neutrinos from the explosion of any of these stars can be used to distinguish a PISN explosion from a core-collapse supernova.

VII Conclusion

Pair-Instability Supernovae represent a intriguing option for the conclusion of stellar evolution in the case of very massive stars. They could also potentially be the source of some observed super-luminous supernovae and early generations of such stars would greatly change the chemical evolution and dynamics of galaxies. Much about PISN and their progenitors remains uncertain and observations of PISN would greatly improve our understanding. Distinguishing a PISN from other supernova types may not be as straight-forward as expected if the lightcurves of PISN at higher metallicities are not super-luminous but rather similar to the light curves of other supernovae. The goal of this paper is to determine whether the neutrino signal from a PISN in the Milky Way or a nearby galaxy could be detected and used to identify the explosion mechanism independent of the electromagnetic emission.

The neutrino emission from a PISN simulation using P150 and P250 progenitors was calculated using the NuLib code. We determined that the dominant emission process was the thermal process of electron positron annihilation into neutrino pairs of all flavors. Also of great importance was the emission via weak processes which produced a 10 MeV peak in the spectrum that is familiar from SNe Ia investigations as coming predominantly from electron capture on copper. The neutrino emission is significant for a duration of $\sim 30$ s, peaks at the time of maximum stellar compression, and the average energy is around 1-2 MeV.

We calculated the neutrino flavor oscillations through the stellar envelope and accounted for decoherence between the SN and Earth. The overall effect of flavor oscillations is to convert most of the electron neutrinos into muon and tau neutrinos. Electron antineutrinos have a greater probability of surviving, especially in NMO. The presence of the shocks in the SN density profile caused time/energy dependent diabatic evolution, which at certain times and energies, had a non-trivial impact on the oscillated flux.

The oscillated flux was then used as input for the code SNOwGLoBES in order to calculate the interaction event rate at various detectors. For a P250 PISN at 10 kpc, we find that HK could measure several tens of events and SK, DUNE and JUNO would measure several events. For the P150 PISN at 10 kpc only HK would detect enough events to observe the explosion. These predictions are not sensitive to the choice of EOS. Thus present and near-future neutrinos detectors can identify a PISN at the Galactic center providing a useful discriminator of the explosion type. The spectral distribution of the events reveals that most events would be below 5 MeV. However, because of the IBD contribution, the spectral distribution of events for the NMO has a double peak structure while the spectral distribution of events for the IMO does not.

Our conclusion is that the gross features of the neutrino signal from a PISN are well-understood, and that the signal contains distinct signatures which are potentially detectable with present neutrino detectors should such a supernova occur in the Milky Way, and with Hyper-Kamiokande if a high mass PISN is located in the Large Magellanic Cloud. Additionally, the signal has spectral features that could help determine the neutrino mass ordering if the supernova is sufficiently close. Further refinement of the model and consideration of the small neglected effects will reduce the uncertainty of the predictions and allow for better extraction of quantitative information should a nearby PISN occur.

VIII Acknowledgments

We are grateful to Evan O’Connor for his help with the implementation of NuLib. This work was supported at NC State by DOE grants DE-FG02-02ER41216 and SC0010263.