The Neutrino Signal From Pair Instability Supernovae
Warren P. Wright, Matthew S. Gilmer, Carla Fröhlich, James P. Kneller
TThe Neutrino Signal From Pair Instability Supernovae
Warren P. Wright, ∗ Matthew S. Gilmer, † Carla Fröhlich, ‡ and James P. Kneller § Department of Physics, North Carolina State University, Raleigh, North Carolina 27695, USA (Dated: July 19, 2018)A very massive star with a carbon-oxygen core in the range of 64 M (cid:12) < M CO <
133 M (cid:12) is expected to undergo a very different kind of explosion known as a pair instability super-nova. Pair instability supernovae are candidates for superluminous supernovae due to theprodigious amounts of radioactive elements they create. While the basic mechanism for theexplosion is understood, how a star reaches a state is not, thus observations of a nearbypair-instability supernova would allow us to test current models of stellar evolution at theextreme of stellar masses. Much will be sought within the electromagnetic radiation we de-tect from such a supernova but we should not forget that the neutrinos from a pair-instabilitysupernova contain unique signatures of the event that unambiguously identify this type ofexplosion. 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 bythis mechanism taking into account the full time and energy dependence of the neutrinoemission and the flavor evolution through the outer layers of the star. We calculate theneutrino 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 easilybe detected in multiple different neutrino detectors at the ‘standard’ supernova distance of10 kpc producing several events in DUNE, JUNE and SuperKamiokande, while the lightestprogenitors only produce a handful of events (if any) in the same detectors. The proposedHyperKamiokande detector would detect neutrinos from a large pair-instability supernovaas far as ∼
50 kpc allowing it to reach the Megallanic Clouds and the several very high massstars known to exist there.
I. INTRODUCTION
Pair instability supernovae (PISNe) are the explosive end point of very massive stars. After central carbonburning, the cores of these stars reach conditions that result in a dynamical instability due to the formationof electron-positron pairs [1–3]. The conversion of photons to electron-positron pairs softens the equation ofstate and leads to a contraction which increases the density and temperature triggering explosive burning ofthe oxygen. The energy release from the burning is sufficient to reverse the contraction and unbind the star in ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] a r X i v : . [ a s t r o - ph . H E ] O c t the form of a PISN explosion [4, 5]. The existence of such supernovae would be important from both a chemicalevolution [6, 7] and galaxy formation [8] perspective. In addition to the large amounts of Ni they produce,the predicted nucleosynthesis pattern from PISNe has a strong odd-even staggering [4]. Such a pattern has notyet been found conclusively in any metal-poor stars [9] but one tentative case has been identified [10].For a star to experience the pair instability (PI), it needs to be massive enough to form a carbon-oxygencore in the range of 64 M (cid:12) < M CO <
133 M (cid:12) [4]. The corresponding zero age main sequence (ZAMS) mass ishighly dependent on the details of stellar evolution and in particular on the mass loss history, which dependson the metallicity and rotation, but is not well understood [11]. The long-standing expectation was that onlynon-rotating stars of zero metallicity (which are thought to have very little mass loss) born with ZAMS massesin the range of 140 M (cid:12) < M
ZAMS <
260 M (cid:12) explode this way [4]. Stars with higher ZAMS mass will collapsedirectly to a black hole; stars with an initial mass less than 140 M (cid:12) may undergo a series of pair instabilitypulses (pulsational PISN) but ultimately explode as core-collapse supernovae. However it has been found thatincluding rotation in stellar models can facilitate a chemically more homogeneous evolution. This leads to largercarbon-oxygen cores for the same initial mass, and hence shifts the mass range of those stars which undergoPISN to a lower ZAMS mass limit. In simulations at zero metallicity and with initial rotations rates of 80% ofKeplerian velocity, stars with an initial masses as low as 40 M (cid:12) are found to undergo a pulsational PISN andstars above 65 M (cid:12) undergo a full PISN [5]. Also more recently, the viewpoint that the star had to be metal freeto explode as a PISN has been challenged. Langer et al. [12] find that PISN may occur at metallicities as highas Z (cid:12) / et al. [12] is that there is one PISN event for every1000 SNe in the local universe. To make this estimate the authors considered only the mass range 140 M (cid:12) 260 M (cid:12) and thus do not include the previously mentioned extension down to lower ZAMS masses inthe presence of rapid rotation. Additionally, while they include the effects of magnetic torques on mixing, theydo not consider a surface dipolar magnetic field and its affect on the mass-loss rate, as in Georgy et al. [13].Due to this uncertainty, we adopt a lower limit for the minimum ZAMS mass of PISN progenitors of 100 M (cid:12) .If that is the case, one can find many candidates of potential PISN progenitors within the Milky Way andnearby galaxies. Several stars with masses exceeding 100 M (cid:12) have been found in the Magellanic Clouds [14].Crowthers et al. [15] identified almost two dozen stars in R136/30 Dor whose initial masses exceed ∼ 100 M (cid:12) .Very massive stars are also known in NGC3603 [16] and the Arches cluster [17]. Some of these stars exceed thepreviously assumed upper mass limit of ∼ 150 M (cid:12) for stars in metal-rich galaxies [18, 19]. If one of these starswere 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 ofauthors have revisited the pair-instability mechanism computing expected lightcurves for various PISN modelsfrom zero metallicity to 10 − Z (cid:12) [20–23]. These studies have found the peak luminosity in the models canbe comparable to the peak luminosity of observed superluminous supernovae. Indeed a few of the observedsuperluminous supernovae with slower decay times have been discussed as possibly being of PISN origin: forexample, the energy required to power the lightcurve of SN 2006gy (at a distance of 73 Mpc) could be providedby the decay of 22 M (cid:12) of Ni [24]; SN 2007bi has been suggested to be a pair-instability supernova by Gal-Yam et al. [25], and another candidate was found by Cooke et al. [26]. Superluminous supernova PS1-14bj also hasa slowly evolving light curve, consistent with PISN models [27]. Of course, in addition to the energy from thedecay of the large of amount of Ni, other processes, such as interactions of the ejecta with a circum-stellermedium, could also contribute to the high luminosity. However we note that the calculated PISN lightcurves atnon-zero metallicity are generally dimmer than the lightcurves at zero metallicity so PISN in the local Universemay become hidden among other classes of supernovae, see e.g. [28] or [29]. Should a nearby PISN occur, onewould expect the increased quality of the observations might be able to find enough differences between a PISNand other supernova types to distinguish them but, as always, extracting information from the electromagneticsignal 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 neutrinosignal because the neutrino emission mechanism for PISN is the same as for Type Ia supernovae [30–34] buton a much larger scale. The aim of this paper is to compute the expected neutrino signal from a PISN anddetermine its detectability in current and near-future neutrino detectors in order to answer the question ofwhether the neutrinos can be added to the suite of tools one can use to identify a PISN. We take into accountthe full time and energy dependence of the emission and the follow the neutrino flavor transformation in theenvelope 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 supernovaused to compute the neutrino signal. In Sec. III we outline the algorithm for computing the neutrino spectrumand luminosity from the simulation and then in Sec. IV we present the calculation of the flavor transformationthrough 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 ofnearby high mass stars and we conclude in Sec. VII. II. SUPERNOVA SIMULATION For the PISN simulation we use two GENEC [35, 36] progenitor models, P250 and P150 [23, 37]. As theirnames suggest, the initial masses of these models are 250 M (cid:12) and 150 M (cid:12) , respectively. Throughout the restof our paper we shall distinguish our calculations by the progenitor model. The models have initial metallicityZ = 10 − and were evolved without rotation. The Schwarzschild criterion was used for convection and theouter convective zone was treated according to the mixing length theory with α MLT = 1 . α MLT = l/H ρ ,where l is the mixing length and H ρ is the density scale height). Nearing the end of core carbon burning (after25,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 (cid:12) which is nearthe upper end of the PISN mass range. The P150 model has a carbon-oxygen core mass of 65.7 M (cid:12) which is atthe lower end of the PISN mass range. Thus the neutrino emission we calculate from these two models shouldbracket what we can expect from an actual PISN.To model the explosion phase, we used version 4.3 of the FLASH multiphysics code [38, 39]. We used a subsetof standard code modules suitable for the present application including the Helmholz EOS and the Aprox19nuclear reaction network [40], the directionally unsplit hydrodynamics solver [41], the improved multipole solverfor self-gravity [42], and its adaptive meshing capability. We run on a 1D grid of spherical geometry with aneffective mesh resolution of 1.3 × 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 defaultvalues. We consider only the innermost 5 × cm of the GENEC progenitor models because the outermaterial 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 preventedfrom flowing into the grid) for the inner and outer boundaries, respectively. We then evolve the models throughthe contraction and explosion phases until nucleosynthesis is complete. Additional details on the progenitormodels and the explosion calculation can be found in Gilmer et al. [37]. The thermodynamic trajectory of thecore 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[33] and [34] for the case of Type Ia SN. We use the software package NuLib [43–50] to calculate the emission asa function of time, neutrino energy and flavor and due to the different emission processes. The weak processesincluded are electron and positron capture on neutrons, protons and nuclei. The only thermal process includedis neutrino pair creation due to electrons annihilating with positrons. This emission calculation is done bypost-processing the simulations described in the previous section. Density, temperature, electron fraction andisotopic composition are extracted from each simulation for each time-slice and radial zone. These quantitiesare 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 (basedon [51]), the second is the SFHo EOS [48]. The strength of the Helmholtz EOS is consistency, because it isthe same EOS used in the original FLASH simulation. However, the simulation (with the Helmholtz EOS andthe Aprox19 nuclear reaction network [40]) does not track all isotopes and thus neutrino emission from weakprocesses on isotopes not included, is missing. Nonetheless, given that the dominant isotopes are tracked, itis expected that these missing nuclei would not account for the bulk of neutrino emission from weak processesbut it might miss spectral features especially towards higher neutrino emission energies. To test this hypothesiswe consider a second EOS, SFHo, which assumes Nuclear Statistical Equilibrium (NSE) thus allowing us toa include emission from a wider range of nuclei. This EOS is an accurate description of the composition ofmaterial in a core-collapse supernova (CCSN) when T > . ≈ . T > T > . T > NuLib . The time denoted as t = 0 sis 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 isbecause increasing compression generally causes an increase in temperature which itself causes greater thermalactivity, including neutrino emission. It is worth emphasizing that the neutrino signal from a PISN has a longduration, ∼ 30 s, much longer than that of Type Ia supernovae [31, 33, 34] and even core-collapse supernovae[52, 53]. In this figure, pair production (all six neutrino flavors), electron capture on nuclei ( ν e ), and electronand positron capture on nucleons ( ν e and ¯ ν e ) are given in red, blue and green respectively (positron captureon nuclei is very subordinate). The results for the P150 (left) and the P250 (right) models are presented aswell as the results for the Helmholtz EOS and the SFHo EOS. Additionally, results when the Helmholtz EOS isassumed and a T > T > . T > T > . T > FIG. 1. PISN total neutrino luminosity as a function of time arising from the various neutrino emission processesconsidered. The results for the P150 simulation are shown on the left, the P250 on the right. Emission due to pairproduction are shown as red and purple lines, electron capture on nuclei are the blue lines, and and electron and positroncapture on nucleons are the green lines. The various temperature cuts used are given in the legend. We have also calculatedthe weak and thermal emission using the Helmholtz EOS using temperature cutoffs in order to make a comparison althoughthis EOS does not fail at lower temperatures in the same way as the SFHo EOS. times during the peak of the neutrino emission, pair emission using SFHo with a T > T > T > FIG. 2. PISN total neutrino luminosity for thermal processes for the P250 simulation. The different thermal processesare color coded as per the legend and the code used to generate each is given in parentheses. The calculation assumes theHelmholtz EOS for all processes. The "Pair (Nulib)" line is a reproduction of the "Helm: Pair" line from Fig. 1. 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 thermalprocesses included in the calculation are pair, photo-, plasma, bremsstrahlung and recombination neutrinoprocesses and were calculated using the code sneut4 which is based on [51] which can be found at http://cococubed.asu.edu/code_pages/nuloss.shtml . These calculations have no spectral information but allowus to determine which processes are important. Figure 2 shows that pair production is by far the most importantthermal 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 theconclusion 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 theneutrino spectrum as a function of neutrino energy for all six flavors and for the time slice around maximumemission. The red (blue) lines represent the P250 (P150) simulation and the solid and dashed lines representthe results using the Helmholtz and SFHo EOSs respectively. The figure shows, that for all flavors other than ν e , the spectra have the same expected thermal shape with a peak at ∼ . ∼ ν e spectra has a strong contribution from weak emission which booststhe overall rate and additionally introduces a spectral feature at 10 MeV when we use the SFHo EOS. Given FIG. 3. PISN neutrino flux spectra. Each curve is the sum of all considered weak and thermal processes at the timesliceof maximum emission. The gray region is the spectra from the SFHo, P250 results at t = 12 . the results from [33] and [34], we recognize this 10 MeV peak as weak emission mostly from electron captureon copper. We find that this 10 MeV peak is present at all time slices (for the SFHo results). However, giventhe aforementioned mismatch between the weak process emission using the two EOSs, our confidence in thepresence of this feature in the spectrum is low initially and increases as the simulation progresses (especiallyfor the P250 simulation).Finally, the time evolution of these spectra are simply described as brightening then fading with little changeof the shape except for the spectral feature around 10 MeV in ν e which becomes relatively more pronouncedlater in the simulation. This behavior is demonstrated by the gray region which is the spectra from the SFHo,P250 results at t = 12 . IV. NEUTRINO FLAVOR TRANSFORMATION In order to accurately determine the neutrino signal from a PISN that reaches Earth, the effects of neutrinooscillation need to be accounted for. This includes both the neutrino flavor oscillations that take place duringneutrino propagation through the stellar mantle and the decoherence that arises as the neutrino propagatesthrough the vacuum between the supernova and Earth. The calculation strategy we have used for the flavortransformations is the same as the one used in [33] and [34].Neutrino oscillation phenomena are calculated by solving the Schrödinger equation. The Schrödinger equa-tion in some basis ( X ) is ı dS ( XX ) dr = H ( X ) S ( XX ) , (1)where H ( X ) is the Hamiltonian in that basis and the evolution matrix S ( Y X ) ( r , r ) connects the neutrino statesin the ( X ) basis at some initial position r to the states in a possibly different basis ( Y ) at r . Note that whenwe are not referring to a specific basis, we shall drop the superscript. The evolution of the antineutrinos isgoverned by an evolution matrix ¯ S which evolves according to the same equation but with a Hamiltonian ¯ H .From the evolution matrix we define the transition probability P ( Y X ) yx = (cid:12)(cid:12)(cid:12) S ( Y X ) yx (cid:12)(cid:12)(cid:12) as the probability for somestate x in the ( X ) basis at r to be in state y in the ( Y ) basis at r . Antineutrino transition probabilitieswill be denoted with an overbar. The two bases commonly used in the literature are the flavor basis – whichhas basis states ν e , ν µ , and ν τ – and the matter basis (often called the mass basis when in vacuum) – whichhas basis states ν , ν , and ν [54]. For future reference, we shall use Greek letters α and β to denote genericflavor states and the Roman symbols i and j to denote generic matter basis states. Again, we denote theantineutrino states in the two bases by an overbar. Using this notation, the neutrino transition probability fora neutrino that was initially state j in the matter basis to be detected as state i in the matter basis is denotedby P (mm) ij = P ( ν j → ν i ).The mixing matrix U V which describes the transformation of the flavor basis to the mass basis is parame-terized by three angles θ , θ & θ , and a CP-violating phase δ CP , and can be written as U V = c s − s c c s e − iδ CP − s e iδ CP c c s − s c 00 0 1 , with c ij = cos θ ij and s ij = sin θ ij . A. Vacuum Hamiltonian In order to solve the Schrödinger equation we need to define the Hamiltonian H ( ¯ 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. In0the mass basis the three neutrino states have masses m , m and m and the vacuum Hamiltonian is diagonal.The vacuum Hamiltonian in the flavor basis, denoted by H ( f ) V , is H ( f ) V = 12 E U V m m 00 0 m U † V , (2)with E the neutrino energy. The mixing matrix angles and phases, together with the squared differences betweenneutrino masses, are part of what determines the nature of the oscillation phenomenology. In this paper weadopt the following values for all results (cid:16) m − m , | m − m | , θ , θ , θ , δ CP (cid:17) = (cid:16) . × − eV , . × − eV , . ◦ , ◦ , ◦ , (cid:17) . (3)where m − m > m − m < H V is simply the complex conjugate of H V . B. 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 ( f ) M is given by H ( f ) M = √ G F n e Diag (1 , , n e is the neutrinodensity and Diag (1 , , 0) 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 ¯ H M is related to H M by ¯ H M = − H M . FIG. 4. Density profiles for various time slices for P150 (left) and for P250 (right). The two horizontal dashed lines arethe 2-flavor MSW densities for 1 MeV neutrinos. The density of the high-density MSW resonance, ρ MSWHigh , we use δm for the mass splitting and θ for the mixing angle, while for the low-density MSW resonance ρ MSWLow , we use δm for themass splitting δm and θ for the mixing angle θ . 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 . ∼ × cm at t ≈ − . ∼ g/ccm for the P250 simulation and similarly at ∼ cm at t ≈ . ∼ g/ccm for the P150 simulation. This shock then propagates outward into regions of lower density. Thedensities at which MSW resonances occur for 1 MeV neutrinos are indicated by the dashed lines. To calculatethe density of the high-density MSW resonance, ρ MSWHigh , we use the two-flavor approximation ρ MSW = m N √ G F Y e (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) δm cos 2 θ E (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (4)and δm for the mass splitting δm in the formula and θ for the mixing angle θ , while for the low-densityMSW resonance ρ MSWLow , we use δm for the mass splitting and θ for the mixing angle. In Eq. 4 m N is thenucleon mass. C. Neutrino self interaction potential Another possible contribution to neutrino flavor oscillation arises due to neutrino self-interaction, i.e. wemust add a third term to H (see [55–57] for reviews where neutrino self interaction is discussed in core-collapsesupernovae). Naively, one might not anticipate this contribution to be important in a PISN because the neutrinodensities are quite low. However, the matter density - shown in Fig. 4 - is also much lower than one finds incore-collapse supernovae and thus it is not immediately obvious whether or not neutrino self-interactions can beignored. In order to check whether self-interaction effects might occur, we compute the strength of the neutrinoself-interaction upon a radially emitted test neutrino. One important detail we need to include in our estimate isthat 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 aninteger i and with a thickness δr i given by the radial grid spacing. For shell i , the neutrino luminosity emittedin flavor α per unit volume is l α,i = L α,i /V i where V i is the shell volume and L α,i the luminosity of the shell inneutrino flavor α . Because the matter densities are much too low to trap neutrinos, the neutrinos are emittedisotropically 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 anangle θ with respect to the ray of the test neutrino, is dµ α,i,θ = √ G F (1 − cos Θ)4 πc‘ (cid:18) l α,i h E α i − l ¯ α,i h E ¯ α i (cid:19) π sin θ r i δr i dθ (5)where G F is the Fermi constant, Θ is the angle between the ray of the test neutrino and the neutrinos emittedfrom the volume element, ‘ is the distance between the volume element and the test neutrino, h E α i is the2 FIG. 5. The various Hamiltonian contributions for time slices around peak emission. Black denotes the matter Hamil-tonian, gray horizontal lines are the vacuum potentials, blue are the electron component of the self-interaction potentialand red the muon or tau self-interaction potential. The results are for the P250 simulation. average neutrino energy and ¯ α represents the antineutrino of flavor α . The distance ‘ is simply ‘ = r i + R − r i R cos θ (6)and the angle Θ is cos Θ = R − r i cos θ‘ (7)A particular shell’s contribution to the total neutrino luminosity is assumed to scale as L α,i ∝ T i in accordancewith the results derived in Burrows, Reddy, & Thompson [46] for pair production neutrinos. We determine the3proportionality constant by constraining the sum over shells to give the total luminosity we have previouslycalculated. The contribution from an entire shell can be obtained by performing the integration over the angle θ and the contribution from the whole star can be obtained by summing over all shells. In this way we calculatethe total strength µ α of the neutrino self-interaction.Our temperature cut-off strategy when computing the emission means that if a shell has a temperature belowthe cut-ff then its neutrino luminosity is zero. Thus, for each time-slice, only shells within a certain radius emitneutrinos. This outermost neutrino-emitting shell is not a neutrino-sphere because the neutrinos are not allemitted from its surface (which is the case for a neutrino-sphere in CCSN). Nonetheless, it is instructive totreat the outermost shell as a neutrino-sphere and compute the self-interaction potential via the usual equationsin the neutrino BULB model (i.e. via Eqn. 40c taken from [58]) to compare to our shell calculation. Thesepure BULB model results are an upper limit for our shell model calculations because, by assuming all neutrinoluminosity comes from the outermost emitting shell, there is greater flux of neutrinos from large angles Θ andthus the self-interaction is artificially increased via the factor (1 − cos Θ) in Eqn. 5.Figure 5 displays the self-interaction potential strengths as a function of radius for the P250 simulation. Thesix subplots represent time-slices around the peak neutrino emission and each subplot displays the ( L α , h E α i )data valid for that time-slice as an inset table. The black line represents the matter potential which is √ G F Y e ρ/m N , the blue and red lines are the electron and muon-or-tau self-interaction potentials respec-tively. The dashed blue and red lines represent the corresponding self-interaction potentials calculated viathe BULB model. The gray horizontal lines represent the vacuum potentials with H V,ij = ∆ m ij / (2 h E i ) andthe average energy being the mean across all flavors. From Fig. 5 we conclude the self-interaction potentialsare many orders of magnitude below the matter potential or vacuum Hamiltonian and indicates we need notconsider neutrino self-interaction effects when calculating the the neutrino flavor evolution. As expected, theself-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 affectingthe low density MSW resonance. A similar analysis was performed for the P150 simulation and the neutrinoself-interactions potential was even weaker than in the P250 simulation. D. 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 inthe vacuum Hamiltonian. The matter basis for the antineutrinos is defined similarly. The significant advantageto using the matter basis over others is that for adiabatic evolution the evolution matrix in this basis is closeto diagonal and transition probabilities P ij are constant.4Indeed we find adiabatic evolution is almost exactly the case at early times because the density profiles ofthe simulations are not steeply falling functions of the radius r . This means the matter basis neutrino transitionprobabilities between states ν i and ν j are approximately unity for i = j and zero for i = j . The same is alsotrue for the antineutrinos. However, adiabatic evolution cannot continue for all epochs because the adiabaticitydepends upon the density gradient and the density profiles from our simulation have moving shocks. Thepresence of shocks changes the adiabaticity and the change is most noticeable at those epochs where the shockis in the vicinity of the high (H) and/or low (L) Mikheyev-Smirnov-Wolfenstein (MSW) resonances [59]. Atthese times there is near-total ‘hopping’ from one matter eigenstate to the other [60] which translates into theneutrino remaining in whichever flavor state it was in before the shock. Thus we expect significant changes inthe flavor transition probabilities oscillations as the shock moves through the mantle of the supernova whichwill, in turn, lead to a change in the number of events we detect here on Earth. For a normal mass orderingboth the H and L resonances affect the neutrinos because the L resonance is in the 12 mixing channel (by ijmixing channel we mean mixing is between matter states ν i and ν 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 resonancemoves to the 13 antineutrino mixing channel [59].Figure 6 shows the neutrino matter basis transition probabilities P (mm) ij ( E ν ) for both mass orderings and asa function of energy for a number of time-slices around peak neutrino emission for the P250 simulation. For thetop panel, which represents normal mass ordering, we see energy-dependent departures from adiabatic behaviorin the 12 and 23 mixing channels. This partial diabatic evolution appears at t ∼ − . ∼ ∼ θ compared to θ .The bottom subplot of Fig. 6 represents neutrino matter basis transition probabilities for inverted massordering. Here we only expect diabatic behavior in the 12 mixing channel because only the L resonance occursfor neutrinos in inverted mass ordering. Indeed that is what is is visible in the bottom subplot and the energyand 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 ¯ P (mm) ij ( E ¯ ν ) as a functionof energy for a number of time-slices around peak antineutrino emission for both mass orderings for the P250simulation. 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 inthe normal mass ordering. This happens because the H resonance switches to the antineutrinos for the invertedmass ordering and is between antineutrino matter states ¯ ν and ¯ ν . A small amount of mixing is seen in the 125 FIG. 6. The neutrino matter basis transition probabilities P (mm) ij ( E ν ) as a function of energy. In both the subplots thetop row, from left to right, shows P (mm)11 , P (mm)12 and P (mm)13 . The middle row from left to right shows P (mm)21 , P (mm)22 and P (mm)23 and the bottom row from left to right shows P (mm)31 , P (mm)32 and P (mm)33 . The mass ordering is normal for the topsubplot, inverted for the bottom subplot and the time is given in the legend. Results are for the P250 simulation. 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 12mixing channel is small ( < FIG. 7. The antineutrino matter basis transition probabilities ¯ P (mm) ij ( E ν ) as a function of energy. The layout is the sameas Fig. 6. Results are for the P250 simulation. E. Flavor basis transition probabilities at Earth Once the neutrinos have propagated through the star and start their long voyage to Earth, we need toaccount for the decoherence of the neutrino wavepacket [61]. Accounting for this effect, the probability p αβ thata neutrino emitted as flavor β in the supernova will be detected as flavor α at Earth is given by p αβ = X i | U V,αi | P (mf) iβ ( R ∗ , R ) , (8)where R represents the radius of the neutrino production point (near the center of the supernova which wetake to be the origin) and R ∗ represents the radius of the outer edge of the supernova. One defines antineutrino7 FIG. 8. p ee for neutrinos and ¯ p ee antineutrinos and both mass orderings. The time slices chosen are representative ofwhen oscillation features occur. The top four subplots are for the P150 simulation and the bottom four subplots are forthe P250 simulation. flavor basis survival probabilities ¯ p αβ by using the matrix ¯ P (mf) iβ ( R ∗ , R ). Note these probabilities p and ¯ p aredifferent from the transition probabilities we discussed in the previous section because they do not come froman evolution matrix. In Fig. 8 the electron neutrino and antineutrino survival probabilities p ee and ¯ p ee as afunction of energy. Neutrino (antineutrino) survival probabilities are displayed in the top (bottom) row foreach PISN simulation and the left (right) column represents the normal (inverted) mass ordering. It is clear8the general trend is that flavor oscillations hamper ν e or ¯ ν e survival. Only for ¯ ν e in NMO does greater that50% of the electron flavor survive. For the P250 simulation, the ‘humps’ which appear for electron neutrinosin the NMO and electron antineutrino in the IMO are mostly caused by the passage of the shock throughthe H resonance. Without the shock the probability that an electron neutrino or electron antineutrino wouldbe 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 theP150 simulation, it is clear that there are very few energy / time dependent oscillation effects. In this case theneutrino flavor evolution is almost entirely adiabatic at all times because the shock forms quite late and at lowdensities. In fact, for the P150 simulation, the shock never passed through the H resonance because it formedbeyond that radius, and moved so slowly that it barely passes through the L resonance before the neutrinoemission subsides. F. The neutrino flux at Earth By combining the flavor basis oscillation probabilities and the neutrino emission spectra, we calculate theneutrino flux seen by neutrino detectors at Earth. These fluxes are given by F α = 14 πd X β p αβ ( E )Φ β ( E ) (9)where d is the distance to the supernova and Φ β ( E ) is the differential spectrum of flavor β at the point ofemission. In what follows we set the distance to the PISN to be d = 10 kpc but will comment about the eventrates at other distances. Figure 9 shows the flux as a function of neutrino energy including oscillations (at peakemission time). By comparing Fig. 9 and Fig. 3 the effects of neutrino oscillation can be easily seen. Onceagain, the red (blue) lines represent the P250 (P150) simulation and the solid and dashed lines represent theresults using the Helmholtz and SFHo EOSs respectively. The gray region is the oscillated spectra from theSFHo, P250 results at t = 12 . ν e flux, is now present in all three neutrino flavors. Furthermore, inall cases, even though the unoscillated electron flavor flux dominates over the other flavors, the oscillated fluxspectrum shows that much of this electron flavor has oscillated into muon and tau flavor in both the neutrinoand the antineutrino case. The effects of the shock are difficult to see on the logarithmic scales used to makethese figures, but they exist between t = − . t = 0 s for the P250 simulation. At these times thechanges, due to diabatic evolution induced by the shock, occur at energies sufficiently close to the energy ofpeak emission to make an appreciable difference. All other diabatic effects of the shock occur at energies wherethe luminosity was low enough that any variation is unnoticeable.9 FIG. 9. Total oscillated neutrino flux from a PISN at 10 kpc. Each curve is the sum of all considered weak and thermalprocesses at the time slice of maximum emission. The top subplot is for NMO and the bottom subplot represents IMO.The gray region is the spectra from the SFHo, P250 results at t = 12 . V. NEUTRINO DETECTION Ultimately, the goal is to be able to measure a neutrino signal in a neutrino detector and by using thatsignal, be able to identify source characteristics. In order to determine what the neutrino signal from our0PISN simulations would be, we use the software package SNOwGLoBES to simulate the expected event ratesin a variety of detectors. Table I lists the detectors considered together with their type and mass. Thedetector description includes the work "type" to indicate that the detector model used in SNOwGLoBES isonly approximate. This is because, for the existing detectors, more accurate descriptions of the detectors arenot 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 detectorsand thus we do not include them for the existing detectors for consistency. All results include a computational0.5 MeV threshold. We also assume perfect detection efficiency, that is, the incoming neutrino’s energy isperfectly reconstructed from the detected particle’s energy. This assumption is again made because efficienciesare 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 theinformation available. At the end of this section, an analysis of the effect of detector thresholds and efficiency(also called smearing) is presented. Detector Type Mass (kt) Super-Kamiokande type: 30% phototube coverage[62] ∗ Water Cherenkov 50Hyper-Kamiokande type[63] Water Cherenkov 374DUNE type detector[64] Liquid Ar 40JUNO type detector[65] Scintillator 20IceCube[66] Water Cherenkov 3500 † TABLE I. A summary of the detector types. ∗ See SNOwGLoBES documentation for discussion on phototube coverage. † For IceCube, the mass given is the ‘effective’ mass. Table II shows the expected interaction event count for all the detectors listed in Table I, except IceCube,for a PISN occurring at 10 kpc. The results are reported for both the P150 and P250 simulations, both EOSsand both mass orderings as well as the case of no neutrino oscillations for comparison. These interaction eventcounts are the total from the whole ∼ 30 s neutrino burst as well as across all energies. For reference, theexpected event rates in these same detectors for a Type Ia supernova at 10 kpc are roughly two orders ofmagnitude smaller [31, 33, 34] than the P250 case, while the event rates for a CCSN are approximately threeorders of magnitude larger [57]. Firstly, we note that the effect of neutrino oscillations is to lower the eventcount in all detectors and all cases. This is because charged current interactions usually produce the strongestsignals and the overall trend of oscillations is to convert electron flavor into muon and tau flavor which aredetected via neutral current. Next, we note that for Super-Kamiokande (SK), Hyper-Kamiokande (HK), andJUNO there are more events for NMO than IMO whereas for DUNE, there are more events for IMO than forNMO. This is because SK, HK and JUNO are sensitive to the inverse beta decay (IBD) channel and DUNEis not. As shown in Figure 9, the ¯ ν e flux is about an order of magnitude greater in NMO than in IMO. Ourinvestigations (detailed below) reveal that it is the IBD channel that is responsible for why the NMO event1 Mass Detector NMO IMO Unoscillated Helm SFHo Helm SFHo Helm SFHoP150 Hyper-Kamiokande 1 . 77 1 . 78 1 . 74 1 . 75 3 . 02 3 . . 24 0 . 24 0 . 23 0 . 23 0 . 40 0 . . 14 0 . 14 0 . 15 0 . 15 0 . 25 0 . . 10 0 . 10 0 . 10 0 . 10 0 . 17 0 . . 23 50 . 08 43 . 32 41 . 98 85 . 70 84 . . 98 6 . 69 5 . 79 5 . 61 11 . 46 11 . . 95 2 . 78 3 . 17 3 . 06 5 . 30 5 . . 13 3 . 00 2 . 48 2 . 40 5 . 06 4 . count is greater than the IMO event count in the detectors that are sensitive to it. Next, Table II shows thatSK, DUNE and JUNO would detect several and HK would detect several tens of interaction events over thespan of the neutrino signal for a PISN at the upper end of the progenitor mass range. Such a detection wouldindicate high-temperature nuclear burning for an extended period of time. This allows us to easily distinguisha PISN from other types of supernovae: the neutrino signal of a CCSN reaches a peak luminosity during theneutronization burst which occurs ∼ 100 ms after core bounce [53], the entire duration of the neutrino emissionfrom a Type Ia is ∼ ∼ 30 seconds and takes ∼ d = 10 kpc and that the P150 simulationproduces 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 caninvestigate the observability of the various temporal and spectral features in the neutrino signal. Figure 10displays the interaction event rates expected to be observed in the detectors under consideration from a PISNat 10 kpc. The top four panels show that, in log scale, NMO and IMO give roughly the same event time-profileat all detectors however, the bottom four panels, using a linear scale, reveal the aforementioned differencesbetween the rates from NMO and IMO. For both orderings our general expectation is for a Gaussian-like burstof neutrinos over a period of ∼ 30 s. Solid and dashed lines represent results from the Helmholtz and SFHoEOS respectively and it is clear that the choice of EOS has little impact. If we look very carefully in the P250simulation we find in the IMO case a slight increase between t = − t = 0 s due to shock induced diabaticevolution. This burst structure shows that at the peak, HK would be seeing a little more than 5 interactionevents per second for the NMO. However, for SK, DUNE and JUNO, bigger time bins (or a nearer SN) would2be required to confidently see activity per bin. As expected, the P150 simulation has a much wider and dimmersignal.Figure 11 shows the energy structure of the neutrino burst released from a 10 kpc PISN as seen in thedetectors under consideration. From all 8 plots, it is clear that the majority of the signal is below 5 MeV in allcases. The top plots reveals that even though the 10 MeV spectral feature from electron capture on copper is FIG. 10. Detector interaction event rate from a PISN at 10 kpc. The top four (bottom four) plots are on a log (linear)vertical axis. The left (right) plots are for NMO (IMO). For the log plots, the dashed horizontal lines show how the eventrates would shift for a closer PISN. The purple line representing the event rate in HK in the bottom four plots has beenrescaled to a tenth of its proper value for plotting convenience. FIG. 11. Detector interaction event differential spectrum. Event count is the total for the whole neutrino burst. The topfour (bottom four) plots are on log (linear) scales. The left (right) plots are for NMO (IMO). The horizontal lines in thetop four plots show how the event rates would shift for a closer PISN. The purple line representing the event rate in HKin the bottom plot has been rescaled to a tenth of its proper value for plotting convenience. FIG. 12. Detector interaction event differential spectrum by channel. Event count is the total for the whole 30s neutrinoburst. The top (bottom) plots are for HK (DUNE) and the left (right) plots are for NMO (IMO). Results are for the P250simulation using the SFHo EOS. We now examine the spectral event rate per channel in order to quantify and illustrate the previous pointsabout the strength of the IBD in NMO over IMO. Figure 12 shows the interaction event count per 0.5 MeVenergy bins for the full ∼ 30 s duration of a P250 PISN at 10 kpc (using the SFHo EOS). The reason that theIBD channel is present for HK but not for DUNE is because it is dominant for detectors built from materialsthat are composed of hydrogen e.g. water and scintillator. As previously noted, Fig. 9 reveals that the ¯ ν e flux isabout an order of magnitude greater in NMO than IMO. This is the reason why the IBD contribution is shownby Fig. 12 to be much larger in NMO than in IMO. This gives rise to the strong double peak seen in Fig. 11 forthe water and scintillator detectors in NMO. Figure 12 also shows that the neutrino elastic scattering on freeelectrons form a significant contribution to the total event count. For HK, contributions from interactions withoxygen are quite minimal but for DUNE, contributions from interactions with argon are important, especiallyat energies above 4 MeV. Below ∼ ∼ ∼ FIG. 13. Detector threshold analysis. The x-axis represents the simulated detector threshold level and the y-axis is thepercentage of events that would be detected above this threshold. All results are for the full ∼ 30 s neutrino signal.The rows differentiate the two mass ordering cases. The left column is for interaction events and the right column is forsmeared events according to the specifications employed by SNOwGLoBES. Results are presented for the Helmholtz EOSonly. In our analysis, Super-K and Hyper-K have identical threshold structure and are here labeled as WC (for WaterCherenkov). Finally, we examine the effects of detector thresholds and energy smearing which, considering that the bulkof the signal is below 5 MeV, are expected to be an important consideration. Figure 13 shows what percentageof events would survive given a particular detector event threshold. This analysis was for interaction events (leftcolumn), for energy-smeared events (right column), and for both mass orderings. The solid and dashed linesare for the P250 and P150 simulation respectively. The figure reveals that a threshold of 5 MeV would reducethe detectable signal in all cases essentially to zero. However, a threshold of 2 MeV means that, for the P250simulation, Hyper-K, Super-K and JUNO would retain more than 50% of their interaction events while DUNEwould retain closer to 40%. The case is worsened when energy-smearing is added but for a threshold of 2 MeVwe see Hyper-K and Super-K still retain more than 50% of their events, JUNO drops to 30-40% (dependenton mass ordering) and DUNE drops to a bit more than 10% retention. In short, energy smearing and detectorthresholds 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 the6spectral 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, showsthat the Hyper-K and Super-K would still see 50% of their smeared events, while for DUNE and JUNO, thesignal is significantly attenuated. A. IceCube Event Rate The PISN neutrino signal in IceCube must be treated differently than the other detectors listed in TableI. The IceCube neutrino detector is located inside the ice that covers Antarctica and, while it is not primarilydesigned to investigate low energy neutrinos, it nonetheless is sensitive to them. These ‘detections’ will consistof an increase of the background hum of low-energy events in the individual Photo-Multiplier Tubes. The eventsIceCube measures from a PISN are mostly from IBD but there is some contribution from elastic scattering offof electrons in the ice. However no direction or energy information of individual neutrino interactions will beextracted so the two event types cannot be distinguished. What can be determined is the overall flux of theneutrinos and the time structure. In order to be detected the PISN neutrino signal needs to be sufficientlyabove statistical fluctuations of the usual background event rate. Thus, even though the PISN will producemany events in IceCube, the metric of detectability is different from the other detector types considered. Figure FIG. 14. The detection confidence of a PISN in the IceCube neutrino detector. The x-axis is the PISN distance and they-axis is the confidence level of detection above background. Both neutrino mass orderings are represented. Results arefor the P250 simulation. The choice of EOS has little impact in this analysis. 14 shows the detection confidence of a P250 PISN at IceCube for a given PISN distance. As expected, giventhe abundance of ¯ ν e and hence IBD events, detection is easier in NMO than in IMO. However, in either massordering, it is clear that a PISN would have to be within one kiloparsec in order to be confidently detected abovebackground at IceCube. The existence of stars with mass high enough to explode as a PISN at these relatively7nearby distances seems quite unlikely. The P150 PISN case is even less likely to be detected by IceCube (3 σ requires the PISN to be within ∼ 70 pc). VI. ASTRONOMICAL CONTEXT Thus far, we assumed that the PISN was located at 10 kpc. This was chosen because it is roughly thedistance to the Galactic center and the distance used in the SN community as the most probable distance forthe next Galactic SN. If we adopt a conservative lower limit for the ZAMS mass of a PISN of 100 M (cid:12) then thisproximal distance is not unrealistically close: within the Milky Way there are several stars with masses around100 M (cid:12) or above. These include:1. several stars in the Arches cluster [67, 68], close to the Galactic center at a distance of d ≈ (cid:12) [72],2. the OB association Cygnus OB2 [73] at a distance of d = 1 . (cid:12) ,3. the mass of the primary in the HD 15558 system, at a distance of d = 2 . ± 50 M (cid:12) [75] (see also [76]),4. the primary star in the binary η Carina which is, coincidentally, also at a distance of d = 2 . (cid:12) [77, 78].Whether these Galactic very high mass stars are too metal rich to explode as a PISN will have to be determinedwhen the event occurs. The neutrino signal will easily distinguish the explosion type. A number of very highmass stars are also known within the Large Magellanic Cloud (LMC) - which has a lower metallicity of 0 . 43 Z (cid:12) - and which is at a distance of d = 49 . 97 kpc [80]:1. The R136 open cluster contains nine stars with masses greater than 100 M (cid:12) at 1-sigma [15].2. NGC 3603-A1 with the more massive component of the binary having a mass 116 ± 31 M (cid:12) [16].3. WR21a whose more massive component has a mass of 103 . ± . (cid:12) [81].For stars at this distance the event rates given in table II need to be scaled downward by a factor of 25. Thuswe see Hyper-Kamiokande is the only detector at the present time or in the near future capable of detectingevents from a large PISN at LMC distances. Again, the neutrinos from the explosion of any of these stars canbe 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 caseof very massive stars. They could also potentially be the source of some observed super-luminous supernovae Computed from log( Z/Z (cid:12) ) = [Fe / H] and the values of [Fe/H] found in [79] NuLib code. We determined that the dominant emission process was the thermal process of electron positronannihilation into neutrino pairs of all flavors. Also of great importance was the emission via weak processes whichproduced a 10 MeV peak in the spectrum that is familiar from SNe Ia investigations as coming predominantlyfrom electron capture on copper. The neutrino emission is significant for a duration of ∼ 30 s, peaks at thetime 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 decoherencebetween the SN and Earth. The overall effect of flavor oscillations is to convert most of the electron neutrinosinto 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, whichat 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 interactionevent rate at various detectors. For a P250 PISN at 10 kpc, we find that HK could measure several tens of eventsand SK, DUNE and JUNO would measure several events. For the P150 PISN at 10 kpc only HK would detectenough events to observe the explosion. These predictions are not sensitive to the choice of EOS. Thus presentand near-future neutrinos detectors can identify a PISN at the Galactic center providing a useful discriminatorof 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 peakstructure 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 thesignal contains distinct signatures which are potentially detectable with present neutrino detectors should sucha supernova occur in the Milky Way, and with Hyper-Kamiokande if a high mass PISN is located in the LargeMagellanic Cloud. Additionally, the signal has spectral features that could help determine the neutrino massordering if the supernova is sufficiently close. Further refinement of the model and consideration of the smallneglected effects will reduce the uncertainty of the predictions and allow for better extraction of quantitativeinformation should a nearby PISN occur.9 VIII. ACKNOWLEDGMENTS We are grateful to Evan O’Connor for his help with the implementation of NuLib . This work was supportedat NC State by DOE grants DE-FG02-02ER41216 and SC0010263. [1] Z. Barkat, G. Rakavy, and N. Sack, Physical Review Letters , 379 (1967).[2] G. Rakavy and G. Shaviv, ApJ , 803 (1967).[3] G. S. Fraley, Astrophysics and Space Science , 96 (1968).[4] A. Heger and S. E. Woosley, ApJ , 532 (2002), astro-ph/0107037.[5] E. Chatzopoulos and J. C. Wheeler, ApJ , 42 (2012), arXiv:1201.1328 [astro-ph.HE].[6] T. Karlsson, J. L. Johnson, and V. Bromm, ApJ , 6-16 (2008), arXiv:0709.4025.[7] J. H. Wise, M. J. Turk, M. L. Norman, and T. Abel, ApJ , 50 (2012), arXiv:1011.2632.[8] D. Whalen, B. van Veelen, B. W. O’Shea, and M. L. Norman, ApJ , 49-67 (2008), arXiv:0801.3698.[9] A. P. Ji, A. Frebel, and V. Bromm, MNRAS , 659 (2015), arXiv:1508.06137.[10] W. Aoki, N. Tominaga, T. C. Beers, S. Honda, and Y. S. Lee, Science , 912 (2014).[11] J. S. Vink, “Mass-loss rates of very massive stars,” in Very Massive Stars in the Local Universe , edited by J. S. Vink(Springer International Publishing, Cham, 2015) pp. 77–111.[12] N. Langer, C. A. Norman, A. de Koter, J. S. Vink, M. Cantiello, and S.-C. Yoon, A&A , L19 (2007),arXiv:0708.1970.[13] C. Georgy, G. Meynet, S. Ekström, G. A. Wade, V. Petit, Z. Keszthelyi, and R. Hirschi, A&A , L5 (2017),arXiv:1702.02340 [astro-ph.SR].[14] V. K. Kudritzki R.-P., Lennon D.J., Haser S.M., Puls J., Pauldrach A.W.A., in Science with the Hubble SpaceTelescope , edited by P. Benvenuti, F.D. Macchetto, and E.J. Schreier.[15] P. A. Crowther, S. M. Caballero-Nieves, K. A. Bostroem, J. Maíz Apellániz, F. R. N. Schneider, N. R. Walborn,C. R. Angus, I. Brott, A. Bonanos, A. de Koter, S. E. de Mink, C. J. Evans, G. Gräfener, A. Herrero, I. D. Howarth,N. Langer, D. J. Lennon, J. Puls, H. Sana, and J. S. Vink, MNRAS , 624 (2016), arXiv:1603.04994 [astro-ph.SR].[16] O. Schnurr, J. Casoli, A.-N. Chené, A. F. J. Moffat, and N. St-Louis, MNRAS , L38 (2008), arXiv:0806.2815.[17] F. Martins, D. J. Hillier, T. Paumard, F. Eisenhauer, T. Ott, and R. Genzel, A&A , 219 (2008), arXiv:0711.0657.[18] D. F. Figer, Nature , 192 (2005), astro-ph/0503193.[19] C. Koen, MNRAS , 590 (2006).[20] D. Kasen, S. E. Woosley, and A. Heger, ApJ , 102 (2011), arXiv:1101.3336 [astro-ph.HE].[21] D. J. Whalen, W. Even, L. H. Frey, J. Smidt, J. L. Johnson, C. C. Lovekin, C. L. Fryer, M. Stiavelli, D. E. Holz,A. Heger, S. E. Woosley, and A. L. Hungerford, ApJ , 110 (2013), arXiv:1211.4979.[22] D. J. Whalen, J. Smidt, W. Even, S. E. Woosley, A. Heger, M. Stiavelli, and C. L. Fryer, ApJ , 106 (2014),arXiv:1311.1221.[23] A. Kozyreva, M. Gilmer, R. Hirschi, C. Fröhlich, S. Blinnikov, R. T. Wollaeger, U. M. Noebauer, D. R. van Rossum,A. Heger, W. P. Even, R. Waldman, A. Tolstov, E. Chatzopoulos, and E. Sorokina, MNRAS , 2854 (2017), arXiv:1610.01086 [astro-ph.HE].[24] N. Smith, W. Li, R. J. Foley, J. C. Wheeler, D. Pooley, R. Chornock, A. V. Filippenko, J. M. Silverman, R. Quimby,J. S. Bloom, and C. Hansen, The Astrophysical Journal , 1116 (2007).[25] A. Gal-Yam, P. Mazzali, E. O. Ofek, P. E. Nugent, S. R. Kulkarni, M. M. Kasliwal, R. M. Quimby, A. V. Filippenko,S. B. Cenko, R. Chornock, R. Waldman, D. Kasen, M. Sullivan, E. C. Beshore, A. J. Drake, R. C. Thomas, J. S.Bloom, D. Poznanski, A. A. Miller, R. J. Foley, J. M. Silverman, I. Arcavi, R. S. Ellis, and J. Deng, Nature ,624 (2009), arXiv:1001.1156.[26] J. Cooke, M. Sullivan, A. Gal-Yam, E. J. Barton, R. G. Carlberg, E. V. Ryan-Weber, C. Horst, Y. Omori, and C. G.Díaz, Nature , 228 (2012), arXiv:1211.2003 [astro-ph.CO].[27] R. Lunnan, R. Chornock, E. Berger, D. Milisavljevic, D. O. Jones, A. Rest, W. Fong, C. Fransson, R. Margutti, M. R.Drout, P. K. Blanchard, P. Challis, P. S. Cowperthwaite, R. J. Foley, R. P. Kirshner, N. Morrell, A. G. Riess, K. C.Roth, D. Scolnic, S. J. Smartt, K. W. Smith, V. A. Villar, K. C. Chambers, P. W. Draper, M. E. Huber, N. Kaiser,R.-P. Kudritzki, E. A. Magnier, N. Metcalfe, and C. Waters, ApJ , 144 (2016), arXiv:1605.05235 [astro-ph.HE].[28] D. J. Whalen, J. Smidt, A. Heger, R. Hirschi, N. Yusof, W. Even, C. L. Fryer, M. Stiavelli, K.-J. Chen, and C. C.Joggerst, ApJ , 9 (2014), arXiv:1312.5360 [astro-ph.HE].[29] A. Kozyreva, S.-C. Yoon, and N. Langer, A&A , A146 (2014), arXiv:1405.6340 [astro-ph.HE].[30] T. Kunugise and K. Iwamoto, PASJ , L57 (2007).[31] A. Odrzywolek and T. Plewa, A&A , A156 (2011), arXiv:1006.0490 [astro-ph.SR].[32] I. R. Seitenzahl, M. Herzog, A. J. Ruiter, K. Marquardt, S. T. Ohlmann, and F. K. Röpke, Phys. Rev. D , 124013(2015), arXiv:1511.02542 [astro-ph.SR].[33] W. P. Wright, G. Nagaraj, J. P. Kneller, K. Scholberg, and I. R. Seitenzahl, Phys. Rev. D , 025026 (2016),arXiv:1605.01408 [astro-ph.HE].[34] W. P. Wright, J. P. Kneller, S. T. Ohlmann, F. K. Röpke, K. Scholberg, and I. R. Seitenzahl, Phys. Rev. D ,043006 (2017), arXiv:1609.07403 [astro-ph.HE].[35] S. Ekström, C. Georgy, P. Eggenberger, G. Meynet, N. Mowlavi, A. Wyttenbach, A. Granada, T. Decressin,R. Hirschi, U. Frischknecht, C. Charbonnel, and A. Maeder, A&A , A146 (2012), arXiv:1110.5049 [astro-ph.SR].[36] N. Yusof, R. Hirschi, G. Meynet, P. A. Crowther, S. Ekström, U. Frischknecht, C. Georgy, H. Abu Kassim, andO. Schnurr, MNRAS , 1114 (2013), arXiv:1305.2099 [astro-ph.SR].[37] M. S. Gilmer, A. Kozyreva, R. Hirschi, C. Fröhlich, and N. Yusof, ApJ , 100 (2017), arXiv:1706.07454 [astro-ph.SR].[38] B. Fryxell, K. Olson, P. Ricker, F. X. Timmes, M. Zingale, D. Q. Lamb, P. MacNeice, R. Rosner, J. W. Truran, andH. Tufo, ApJS , 273 (2000).[39] A. Dubey, L. B. Reid, K. Weide, K. Antypas, M. K. Ganapathy, K. Riley, D. Sheeler, and A. Siegal, ArXiv e-prints(2009), arXiv:0903.4875 [cs.SE].[40] F. X. Timmes and F. D. Swesty, ApJS , 501 (2000).[41] D. Lee and A. E. Deane, Journal of Computational Physics , 952 (2009).[42] S. M. Couch, C. Graziani, and N. Flocke, ApJ , 181 (2013), arXiv:1307.3135 [astro-ph.HE].[43] G. M. Fuller, W. A. Fowler, and M. J. Newman, ApJ , 715 (1982).[44] K. Langanke and G. Martínez-Pinedo, Nuclear Physics A , 481 (2000), nucl-th/0001018. [45] K. Langanke, G. Martínez-Pinedo, J. M. Sampaio, D. J. Dean, W. R. Hix, O. E. Messer, A. Mezzacappa,M. Liebendörfer, H.-T. Janka, and M. Rampp, Physical Review Letters , 241102 (2003), astro-ph/0302459.[46] A. Burrows, S. Reddy, and T. A. Thompson, Nuclear Physics A , 356 (2006), astro-ph/0404432.[47] T. Oda, M. Hino, K. Muto, M. Takahara, and K. Sato, Atomic Data and Nuclear Data Tables , 231 (1994).[48] A. W. Steiner, M. Hempel, and T. Fischer, ApJ , 17 (2013), arXiv:1207.2184 [astro-ph.SR].[49] E. O’Connor and C. D. Ott, Classical and Quantum Gravity , 114103 (2010), arXiv:0912.2393 [astro-ph.HE].[50] C. Sullivan, E. O’Connor, R. G. T. Zegers, T. Grubb, and S. M. Austin, ApJ , 44 (2016), arXiv:1508.07348[astro-ph.HE].[51] N. Itoh, H. Hayashi, A. Nishikawa, and Y. Kohyama, ApJS , 411 (1996).[52] T. Totani, K. Sato, H. E. Dalhed, and J. R. Wilson, ApJ , 216 (1998), astro-ph/9710203.[53] T. Fischer, S. C. Whitehouse, A. Mezzacappa, F.-K. Thielemann, and M. Liebendörfer, A&A , A80 (2010),arXiv:0908.1871 [astro-ph.HE].[54] H. A. Bethe, Physical Review Letters , 1305 (1986).[55] H. Duan and J. P. Kneller, Journal of Physics G Nuclear Physics , 113201 (2009), arXiv:0904.0974 [astro-ph.HE].[56] H. Duan, G. M. Fuller, and Y.-Z. Qian, Annual Review of Nuclear and Particle Science , 569 (2010),arXiv:1001.2799 [hep-ph].[57] A. Mirizzi, I. Tamborra, H.-T. Janka, N. Saviano, K. Scholberg, R. Bollig, L. Hudepohl, and S. Chakraborty, ArXive-prints (2015), arXiv:1508.00785 [astro-ph.HE].[58] H. Duan, G. M. Fuller, J. Carlson, and Y.-Z. Qian, Phys. Rev. D , 105014 (2006), astro-ph/0606616.[59] A. S. Dighe and A. Y. Smirnov, Phys. Rev. D , 033007 (2000), hep-ph/9907423.[60] S. J. Parke, Physical Review Letters , 1275 (1986).[61] C. Giunti, C. W. Kim, and U. W. Lee, Physics Letters B , 237 (1998), hep-ph/9709494.[62] S. Fukuda, Y. Fukuda, T. Hayakawa, E. Ichihara, M. Ishitsuka, Y. Itow, T. Kajita, J. Kameda, K. Kaneyuki,S. Kasuga, K. Kobayashi, Y. Kobayashi, Y. Koshio, M. Miura, S. Moriyama, M. Nakahata, S. Nakayama, T. Namba,Y. Obayashi, A. Okada, M. Oketa, K. Okumura, T. Oyabu, N. Sakurai, M. Shiozawa, Y. Suzuki, Y. Takeuchi,T. Toshito, Y. Totsuka, S. Yamada, S. Desai, M. Earl, J. T. Hong, E. Kearns, M. Masuzawa, M. D. Messier, J. L.Stone, L. R. Sulak, C. W. Walter, W. Wang, K. Scholberg, T. Barszczak, D. Casper, D. W. Liu, W. Gajewski,P. G. Halverson, J. Hsu, W. R. Kropp, S. Mine, L. R. Price, F. Reines, M. Smy, H. W. Sobel, M. R. Vagins,K. S. Ganezer, W. E. Keig, R. W. Ellsworth, S. Tasaka, J. W. Flanagan, A. Kibayashi, J. G. Learned, S. Matsuno,V. J. Stenger, Y. Hayato, T. Ishii, A. Ichikawa, J. Kanzaki, T. Kobayashi, T. Maruyama, K. Nakamura, Y. Oyama,A. Sakai, M. Sakuda, O. Sasaki, S. Echigo, T. Iwashita, M. Kohama, A. T. Suzuki, M. Hasegawa, T. Inagaki, I. Kato,H. Maesaka, T. Nakaya, K. Nishikawa, S. Yamamoto, T. J. Haines, B. K. Kim, R. Sanford, R. Svoboda, E. Blaufuss,M. L. Chen, Z. Conner, J. A. Goodman, E. Guillian, G. W. Sullivan, D. Turcan, A. Habig, M. Ackerman, F. Goebel,J. Hill, C. K. Jung, T. Kato, D. Kerr, M. Malek, K. Martens, C. Mauger, C. McGrew, E. Sharkey, B. Viren,C. Yanagisawa, W. Doki, S. Inaba, K. Ito, M. Kirisawa, M. Kitaguchi, C. Mitsuda, K. Miyano, C. Saji, M. Takahata,M. Takahashi, K. Higuchi, Y. Kajiyama, A. Kusano, Y. Nagashima, K. Nitta, M. Takita, T. Yamaguchi, M. Yoshida,H. I. Kim, S. B. Kim, J. Yoo, H. Okazawa, M. Etoh, K. Fujita, Y. Gando, A. Hasegawa, T. Hasegawa, S. Hatakeyama,K. Inoue, K. Ishihara, T. Iwamoto, M. Koga, I. Nishiyama, H. Ogawa, J. Shirai, A. Suzuki, T. Takayama, F. Tsushima,M. Koshiba, Y. Ichikawa, T. Hashimoto, Y. Hatakeyama, M. Koike, T. Horiuchi, M. Nemoto, K. Nishijima, H. Takeda, H. Fujiyasu, T. Futagami, H. Ishino, Y. Kanaya, M. Morii, H. Nishihama, H. Nishimura, T. Suzuki, Y. Watanabe,D. Kielczewska, U. Golebiewska, H. G. Berns, S. B. Boyd, R. A. Doyle, J. S. George, A. L. Stachyra, L. L. Wai, R. J.Wilkes, K. K. Young, H. Kobayashi, and Super-Kamiokande Collaboration, Nuclear Instruments and Methods inPhysics Research A , 418 (2003).[63] (2016).[64] R. Acciarri, M. A. Acero, M. Adamowski, C. Adams, P. Adamson, S. Adhikari, Z. Ahmad, C. H. Albright, T. Alion,E. Amador, and et al., ArXiv e-prints (2016), arXiv:1601.02984 [physics.ins-det].[65] F. An, G. An, Q. An, V. Antonelli, E. Baussan, J. Beacom, L. Bezrukov, S. Blyth, R. Brugnera, M. Buizza Avanzini,J. Busto, A. Cabrera, H. Cai, X. Cai, A. Cammi, G. Cao, J. Cao, Y. Chang, S. Chen, S. Chen, Y. Chen, D. Chiesa,M. Clemenza, B. Clerbaux, J. Conrad, D. D’Angelo, H. De Kerret, Z. Deng, Z. Deng, Y. Ding, Z. Djurcic, D. Dornic,M. Dracos, O. Drapier, S. Dusini, S. Dye, T. Enqvist, D. Fan, J. Fang, L. Favart, R. Ford, M. Göger-Neff, H. Gan,A. Garfagnini, M. Giammarchi, M. Gonchar, G. Gong, H. Gong, M. Gonin, M. Grassi, C. Grewing, M. Guan,V. Guarino, G. Guo, W. Guo, X.-H. Guo, C. Hagner, R. Han, M. He, Y. Heng, Y. Hsiung, J. Hu, S. Hu, T. Hu,H. Huang, X. Huang, L. Huo, A. Ioannisian, M. Jeitler, X. Ji, X. Jiang, C. Jollet, L. Kang, M. Karagounis, N. Kazar-ian, Z. Krumshteyn, A. Kruth, P. Kuusiniemi, T. Lachenmaier, R. Leitner, C. Li, J. Li, W. Li, W. Li, X. Li, X. Li,Y. Li, Y. Li, Z.-B. Li, H. Liang, G.-L. Lin, T. Lin, Y.-H. Lin, J. Ling, I. Lippi, D. Liu, H. Liu, H. Liu, J. Liu,J. Liu, J. Liu, Q. Liu, S. Liu, S. Liu, P. Lombardi, Y. Long, H. Lu, J. Lu, J. Lu, J. Lu, B. Lubsandorzhiev, L. Lud-hova, S. Luo, V. Lyashuk, R. Möllenberg, X. Ma, F. Mantovani, Y. Mao, S. M. Mari, W. F. McDonough, G. Meng,A. Meregaglia, E. Meroni, M. Mezzetto, L. Miramonti, T. Mueller, D. Naumov, L. Oberauer, J. P. Ochoa-Ricoux,A. Olshevskiy, F. Ortica, A. Paoloni, H. Peng, J.-C. Peng, E. Previtali, M. Qi, S. Qian, X. Qian, Y. Qian, Z. Qin,G. Raffelt, G. Ranucci, B. Ricci, M. Robens, A. Romani, X. Ruan, X. Ruan, G. Salamanna, M. Shaevitz, V. Sinev,C. Sirignano, M. Sisti, O. Smirnov, M. Soiron, A. Stahl, L. Stanco, J. Steinmann, X. Sun, Y. Sun, D. Taichenachev,J. Tang, I. Tkachev, W. Trzaska, S. van Waasen, C. Volpe, V. Vorobel, L. Votano, C.-H. Wang, G. Wang, H. Wang,M. Wang, R. Wang, S. Wang, W. Wang, Y. Wang, Y. Wang, Y. Wang, Z. Wang, Z. Wang, Z. Wang, Z. Wang,W. Wei, L. Wen, C. Wiebusch, B. Wonsak, Q. Wu, C.-E. Wulz, M. Wurm, Y. Xi, D. Xia, Y. Xie, Z.-z. Xing, J. Xu,B. Yan, C. Yang, C. Yang, G. Yang, L. Yang, Y. Yang, Y. Yao, U. Yegin, F. Yermia, Z. You, B. Yu, C. Yu, Z. Yu,S. Zavatarelli, L. Zhan, C. Zhang, H.-H. Zhang, J. Zhang, J. Zhang, Q. Zhang, Y.-M. Zhang, Z. Zhang, Z. Zhao,Y. Zheng, W. Zhong, G. Zhou, J. Zhou, L. Zhou, R. Zhou, S. Zhou, W. Zhou, X. Zhou, Y. Zhou, Y. Zhou, andJ. Zou, Journal of Physics G Nuclear Physics , 030401 (2016), arXiv:1507.05613 [physics.ins-det].[66] R. Abbasi, Y. Abdou, T. Abu-Zayyad, M. Ackermann, J. Adams, J. A. Aguilar, M. Ahlers, M. M. Allen, D. Altmann,K. Andeen, and et al., A&A , A109 (2011), arXiv:1108.0171 [astro-ph.HE].[67] T. Nagata, C. E. Woodward, M. Shure, and N. Kobayashi, AJ , 1676 (1995).[68] A. S. Cotera, E. F. Erickson, S. W. J. Colgan, J. P. Simpson, D. A. Allen, and M. G. Burton, ApJ , 750 (1996).[69] M. J. Reid, Ann. Rev. Astron. Astrophys. , 345 (1993).[70] Z. Malkin, ArXiv e-prints (2012), arXiv:1202.6128 [astro-ph.GA].[71] T. Camarillo, V. Mathur, T. Mitchell, and B. Ratra, ArXiv e-prints (2017), arXiv:1708.01310.[72] G. Gräfener, J. S. Vink, A. de Koter, and N. Langer, A&A , A56 (2011), arXiv:1106.5361 [astro-ph.SR].[73] P. Massey and A. B. Thompson, AJ , 1408 (1991). [74] M. T. Schuster, M. Marengo, J. L. Hora, G. G. Fazio, R. M. Humphreys, R. D. Gehrz, P. M. Hinz, M. A. Kenworthy,and W. F. Hoffmann, The Astrophysical Journal , 1423 (2009).[75] M. De Becker, G. Rauw, J. Manfroid, and P. Eenens, A&A , 1121 (2006), astro-ph/0606379.[76] C. D. Garmany and P. Massey, Publications of the Astronomical Society of the Pacific , 500 (1981).[77] A. Kashi and N. Soker, ApJ , 602 (2010), arXiv:0912.1439 [astro-ph.SR].[78] N. Clementel, T. I. Madura, C. J. H. Kruip, J.-P. Paardekooper, and T. R. Gull, MNRAS , 2445 (2015),arXiv:1412.7569 [astro-ph.SR].[79] S. Choudhury, A. Subramaniam, and A. A. Cole, MNRAS , 1855 (2016), arXiv:1510.05769.[80] G. Pietrzyński, D. Graczyk, W. Gieren, I. B. Thompson, B. Pilecki, A. Udalski, I. Soszyński, S. Kozłowski,P. Konorski, K. Suchomska, G. Bono, P. G. P. Moroni, S. Villanova, N. Nardetto, F. Bresolin, R. P. Kudritzki,J. Storm, A. Gallenne, R. Smolec, D. Minniti, M. Kubiak, M. K. Szymański, R. Poleski, Ł. Wyrzykowski, K. Ulaczyk,P. Pietrukowicz, M. Górski, and P. Karczmarek, Nature , 76 (2013), arXiv:1303.2063 [astro-ph.GA].[81] F. Tramper, H. Sana, N. E. Fitzsimons, A. de Koter, L. Kaper, L. Mahy, and A. Moffat, MNRAS455