Strong Supernova 1987A Constraints on Bosons Decaying to Neutrinos

Majoron-like bosons would emerge from a supernova (SN) core by neutrino coalescence of the form $\nu\nu\to\phi$ and $\bar\nu\bar\nu\to\phi$ with 100 MeV-range energies. Subsequent decays to (anti)neutrinos of all flavors provide a flux component with energies much larger than the usual flux from the"neutrino sphere."The absence of 100 MeV-range events in the Kamiokande-II and Irvine-Michigan-Brookhaven signal of SN 1987A implies that less than 1% of the total energy was thus emitted and provides the strongest constraint on the Majoron-neutrino coupling of $g\lesssim 10^{-9}\,{\rm MeV}/m_\phi$ for $100~{\rm eV}\lesssim m_\phi\lesssim 100~{\rm MeV}$. It is straightforward to extend our new argument to other hypothetical feebly interacting particles.

If the FIPs interact so strongly that they are trapped themselves or decay before leaving the SN, they contribute to energy transfer [19] and may strongly affect overall SN physics and the explosion mechanism.A class of low-explosion-energy SNe provides particularly strong constraints on such scenarios [20].FIPs on the trapping side of the SN-excluded regime are often constrained by other arguments, although allowed gaps may remain, such as the historical "hadronic axion window" or more recently the "cosmic triangle" for axion-like particles, both meanwhile closed.
In other cases, FIP decays include active neutrinos.In the free-streaming limit, FIPs escape from the inner SN core and so their decays provide 100-MeV-range events, much larger than the usual neutrino burst of few 10 MeV that emerges from the "neutrino sphere" at the edge of the SN core.The background of atmospheric muons has yet larger energies and so the new signal would stick out in a future SN neutrino observation.This argument was first advanced in Ref. [7], and offers an intriguing future detection opportunity.
Our main point is that, by the same token, SN 1987A already provides restrictive limits because the legacy data do not sport any events with such intermediate energies.This constraint, which is available today without the need to wait for the next galactic SN, is far more restrictive than the traditional energy-loss argument.
We illustrate our new argument with the simple case of nonstandard or "secret" neutrino-neutrino interactions [4][5][6][7][8], mediated by a (pseudo)scalar φ (mass m φ ) that we FIG. 1. Constraints on the Majoron coupling in the m φ -g φ m φ plane from SN 1987A energy loss (green) and the absence of 100 MeV-range ("high-E") events (blue).The shaded range brackets the cold (upper curves) vs. hot (lower curves) SN models, i.e., the Garching muonic models SFHo-18.8 and LS220-s20.0[29].Above the dashed line, Majorons with a reference kinetic energy of 100 MeV decay before leaving the SN core.The "ceiling" of the energy-loss bound is probably outside this figure, but we are not confident about its exact location.The schematic big bang nucleosynthesis (BBN) bounds are taken from Fig. 1 of Ref. [30], based on the cosmic radiation density.Somewhat more restrictive limits may follow from the cosmic microwave background (CMB) (see text).
call Majoron and take to interact with all flavors with the same strength g.We consider m φ > ∼ 100 eV so that neutrino masses and refractive matter potentials can be ignored.The lepton-number violating production channels ν ν → φ and νν → φ and corresponding decays yield the constraints previewed in Fig. 1.
Majoron decay and production.-Auniversal ν-ν interaction by Majoron exchange is given by [39] where ψ ν is a two-component Majorana field and g a real number.In the relativistic limit we refer to the Majorana helicity states as ν and ν in the usual sense.
The decay into pairs of relativistic neutrinos requires equal helicities, implying the lepton-number violating channels φ → νν or ν ν.Each individual rate is which includes a symmetry factor 1/2 for identical finalstate particles.(We always use natural units with h = c = k B = 1.)The total rate requires a factor of 6 for six species [40].For a relativistic Majoron, this rate is slower by the Lorentz factor m φ /E φ , implying that the laboratory decay rate depends only on the combination gm φ .
The requirement that Majorons with E φ = 100 MeV decay beyond the neutrino-sphere radius of 20 km thus implies gm φ < ∼ 10 −7 MeV, shown as a dashed line in Fig. 1.On the other hand, the decay neutrinos should not be delayed by more than a few seconds.The requirement MeV for E φ = 100 MeV.The time-of-flight difference is much smaller for relativistic Majorons, so for the constraints shown in Fig. 1 the signals are indeed contemporaneous, although somewhat marginally for m φ around 100 MeV.
The neutrino decay spectrum is flat between In a neutrino gas of one species α, occupation number f α (E ν ), the spectral Majoron emission rate from ν α ν α coalescence then is For local thermal equilibrium with temperature T and neutrino chemical potential µ α , the corresponding Fermi-Dirac distribution is f α (E ν ) = e (Eν −µα)/T + 1 −1 .The chemical potential for a flavor ν enters with opposite sign, depending on α denoting a ν or ν.Notice that the lepton-number violation caused by the φ interaction implies µ ν = 0 in true equilibrium.All Majorons decay close to the SN equally into all six neutrino species with a flat spectrum.Therefore, the effective single-species spectral neutrino emission rate is The minimal E φ to produce a neutrino of energy E ν is The first factor of 2 is for two neutrinos per decay, whereas 1/6 appears because this is the rate into one of six species.
One-zone SN model.-For a first estimate we use a one-zone model of the collapsed SN core with a chemical potential µ ν = 100 MeV for ν e and vanishing for the other flavors, volume (4π/3)R 3 with R = 10 km for the emitting region, and duration for substantial deleptonization of τ = 1 s [41].After collapse, the SN core is cold (T 10 MeV) and heats up from outside in as the material deleptonizes.Majoron emission is thus from the coalescence of ν e ν e alone which we take as perfectly degenerate.(In contrast, novel particle emission usually becomes large only after the SN core has heated up at around 1 s after collapse [24].) For m φ = 0 the integral in Eq. ( 3) is a "triangle function" that rises linearly to the value µ ν at E φ = µ ν and then decreases linearly to zero at E φ = 2µ ν .The energy-loss rate per unit volume is 3 with L ν 2 × 10 52 erg/s as recommended by a simple recipe [2] implies gm φ < ∼ 4π 3L ν /R 3 µ 3 ν = 5.5 × 10 −9 MeV.Likewise, the effective ν α production rate per unit volume is Ṅα = (g 2 m 2 φ /64π 3 ) µ 2 ν /3 and therefore the total emitted number is N α = Ṅα (4π/3)R 3 τ .The fluence at Earth is N α /(4πd 2 SN ) where d SN = 49.6 kpc is the distance to SN 1987A [66].The largest detector was IMB with a fiducial mass of 6.8 kton [15] and thus N p = 4.5 × 10 32 fiducial protons.The detection cross section is very roughly σ σE 2 ν with σ 10 −43 cm 2 /MeV 2 and E 2 ν = 7µ 2 ν /18.The total number of 100-MeV-range events therefore is MeV. Numerical SN models.-Thisconstraint is much more restrictive than from energy loss, motivating a detailed study.To this end we use the Garching 1D models SFHo-18.8 and LS220-s20.0that were evolved with the Prometheus Vertex code with six-species neutrino transport [67].These muonic models were recently also used for other particle constraints [24,29].With different final neutron-star masses and different equations of state, these models were taken to span the extremes of a cold and a hot case, reaching internal T of around 40 vs. 60 MeV.On the other hand, the initial µ νe profiles are much more similar, in both cases around 150 MeV in the center and a "lepton core" reaching up to around 10 km.The lepton number of the outer core layers is released within a few ms after core bounce in the form of the prompt ν e burst.More details about these models are provided in the Supplemental Material [42].
SN neutrinos follow a quasi-thermal spectrum that can be represented by a Gamma distribution [68][69][70].We thus write the time-integrated spectrum in the form where E tot is the total SN energy release, E 0 the average νe energy, α a parameter that would be 2 for a Maxwell-Boltzmann distribution, and Γ the Gamma function, not to be confused with a Gamma distribution.3) that we correct for gravitational redshift through the tabulated lapse factors as described in Ref. [24].In the cold model, we find a Majoron luminosity at 1 s post bounce of L φ (1 s) = (gm MeV ) 2 6.46 × 10 68 erg/s, where m MeV = m φ /MeV.According to the traditional SN 1987A cooling argument [2,24,71] we compare it with L ν (1 s) = 4.40×10 52 erg/s, leading to gm φ < 0.83 × 10 −8 MeV shown in Fig. 1.For larger masses, we include a cutoff for those Majorons that are produced with insufficient energy to escape the gravitational potential as explained in the Supplemental Material of Ref. [20].The total emission is E tot φ = (gm MeV ) 2 1.94 × 10 69 erg and nominally E tot ν = E tot φ for gm φ = 0.99 × 10 −8 MeV, practically identical to the luminosity comparison at 1 s.
For the hot model we find L φ (1 s) = (gm MeV )  [9][10][11] and IMB (6.8 kton) [14][15][16].They observed events with energies up to 40 MeV via inverse beta decay νe + p → e + + n, whereas elastic scattering on electrons is small (but dominates for solar ν e detection).For our 100 MeV-range energies, charged current (CC) reactions on oxygen of the form MeV. For energies above the muon production threshold (m µ = 105.7 MeV), the corresponding muonic CC processes also happen, especially of course for atmospheric neutrinos at yet larger energies.Muons quickly come to rest by ionization and produce "Michel e ± " with a characteristic spectrum ending at 53 MeV, half the muon mass.Below the muon Cherenkov threshold of about 160 MeV, they are termed "invisible muons."(For more details about these processes see the Supplemental Material [42].) Figure 2 shows the spectral fluence (time-integrated flux) for the standard SN neutrinos from the cold model, averaged over νe , νµ and ντ .The energy-integrated fluence is 5.10 × 10 9 cm −2 for one species.We also show the corresponding e ± spectrum in the detector; the total event number is 5.07 per kton (for 100% detection efficiency).Next we show the ν spectrum from φ decay which is the same in every species; the total fluence in one species is (gm MeV ) 2 1.90 × 10 25 cm −2 .The e ± event number times (gm MeV ) 2 /kton is 3.62 × 10 17 produced by νe and ν e in CC reactions and 0.37 × 10 17 from Michel e ± (E < ∼ 53 MeV) caused by invisible muons, and a total of 3.99 × 10 17 .
Above the muon Cherenkov threshold of 160 MeV, and assuming the same detection efficiency as for e ± , visible µ ± contribute another 11% to the total events.After each such event, the IMB detector would be blind by trigger dead time, so we should not include the subsequent Michel events.However, even for µ ± themselves, the Cherenkov threshold behavior and the detection efficiency are not available.Therefore, we do not include visible muons, making our Majoron bounds more conservative by some 5%.
A single event with 100% detection efficiency in IMB thus requires gm φ = 6.06 × 10 −10 MeV.For the hot model, the corresponding result is gm φ = 3.71 × 10 −10 MeV, both smaller than the estimate from the onezone model, where we underestimated the cross section.Once more, the exact SN model is not crucial and we essentially find the limits shown in Fig. 1.
Analysis of SN 1987A data.-Wenow turn to a detailed analysis of the Kamiokande II and IMB data.We summarize several details in the Supplemental Material [42] and here only remark that event information was recorded depending on a hardware trigger.In an off-line analysis, one searched for low-energy few-seconds event clusters."Low energy" was defined in Kamiokande-II as less than 170 photo electrons in the inner detector or E e < ∼ 50 MeV [9][10][11], whereas IMB used maximally 100 PMTs firing or E e < ∼ 75 MeV [14][15][16].However, as discussed in Supplemental Material [42], we can conclude that no high-energy events were actually observed even above these thresholds during the SN 1987A burst.
The events from φ decay overlap with the standard SN signal, so one should perform a maximum likelihood analysis with g and m φ as fit parameters.However, the standard SN signal depends on the chosen SN model.For example, our cold (hot) model (using the average νe -ν µντ spectrum) would have produced 9.12 (21.3) events in Kamiokande II with average detected electron energy of 20.1 (22.6)MeV, to be compared with the actually observed 12 events with 14.7 MeV average energy.In IMB they would have produced 3.49 (12.5) events on average with 31.3 (34.4) MeV, to be compared with 8 events with 31.9MeV average.Neither of these models fits the data well and the Kamiokande II and IMB data are themselves in tension with each other, although in terms of the E tot -E 0 -α parameters one finds credible overlapping values [72,73].
We do not have a suite of SN models that would allow us to find the one that best fits the SN 1987A data.Instead we represent the signal in the form of Eq. ( 5) and use an unbinned likelihood for the energies of the events in each detector, as defined in the Supplemental Material [42].First we verify that the maximum of the likelihood for both experiments is at g = 0, i.e., neither of them prefers the new signal.Next we marginalize the combined likelihood by maximizing it for each value of g and m φ over E 0 and E tot .This guarantees our constraints to be conservative, because for each choice of the Majoron parameters we choose the SN neutrino spectral shape as the one that maximizes the agreement with the data.We then follow the procedure outlined in Ref. [74] to set upper bounds on the Majoron coupling for each value of the Majoron mass; more details on our statistical procedure are given in the Supplemental Material [42].We show the corresponding constraints, dominated by the IMB data, in Fig. 1.
Discussion and outlook.-Wehave considered FIPs that escape from the inner SN core and later decay into active neutrinos.Our main result is that the lack of 100-MeV-range events in the SN 1987A data provides surprisingly restrictive constraints.Specifically, the energy loss by νν → φ Majoron emission must be less than 1% of the total binding energy, much more restrictive than the usual SN 1987A cooling limit.
Moreover, our new bound depends mainly on emission during the first second and not on the sparse latetime events or the predicted cooling speed that depends, e.g., on PNS convection.Our result is also insensitive to a concern that the SN 1987A neutron star has not yet been found (see however [75,76]) and that the late events could have been caused by black-hole accretion [77].(See however [29] for a rebuttal of this scenario.) Our limit implies that the impact on SN physics and the explosion mechanism is small.However, our discussion leaves open what happens for much stronger couplings when Majorons do not freely escape.The SN core could deleptonize already during infall, perhaps preventing a successful explosion.On the other hand, a thermal bounce may still occur [35,78].If the interactions are yet stronger, neutrinos and Majorons form a viscous fluid that is more strongly coupled to itself than to the nuclear medium.This peculiar case was recently examined [8]; the SN 1987A signal may exclude a certain range of parameters beyond the upper edge of Fig. 1.
For m φ < ∼ 1 MeV, the cosmic radiation density measured by BBN provides comparable bounds (Fig. 1 of Ref. [30], see also Refs.[79][80][81]), and those from the CMB may be more restrictive, but the exact reach in mass and coupling strength was not directly provided.Having different systematic issues, the cosmological and SN 1987A arguments are nicely complementary for m φ < ∼ 1 MeV, whereas the SN 1987A sensitivity is unique for larger m φ .
Our method can be applied to any class of FIPs decaying to neutrinos.Examples include heavy neutral leptons [82,83] and gauge bosons arising from new symmetries like U (1) Lµ−Lτ [84,85], which can be further constrained relative to the existing bounds from energy loss [86,87].Notice also that bosons coupling exclusively to neutrinos have different production rates if the coalescence process is lepton-number conserving (ν ν → φ) or violating (νν → φ) because in the PNS core, the neutrino and antineutrino distributions differ.
At present it remains open if there exist allowed Majoron parameters somewhere in the trapping regime, a question left for future study.Couplings below our limit leave open the exciting possibility of a detection in the neutrino signal of a future galactic SN [7] that would reveal FIP emission from the inner SN core.
Note Added.-Sinceour paper had appeared on arXiv, our new argument was used to constrain the heavy-lepton model of Ref. [88].
[71] Scaling the original axion bounds to the Majoron case with that simple recipe ignores that here the particle emission is largest directly after core bounce, whereas in the axion or similar cases, the core first has to heat up and the emission is largest perhaps around 1 s post bounce.Moreover, Majorons remove both energy and lepton number.We suspect that the impact on the SN 1987A neutrino signal would be larger than implied by scaling the axion case.A detailed analysis would require including Majoron losses in self-consistent SN models.In view of the much more restrictive counting-rate argument, this exercise is not needed and we can post-process existing models.

Supplemental Material for the Paper Strong Supernova 1987A Constraints on Bosons Decaying to Neutrinos
We summarize some details about the detection cross sections for SN neutrinos in a water Cherenkov detector used in our analysis, the historical SN 1987A observations, our statistical analysis, and the Garching SN models.

A. Detection cross sections
The primary channel for neutrino detection from SN 1987A was inverse beta decay (IBD) νe + p → e + + n on the hydrogen nuclei of the water molecules.Neglecting the recoil of the nucleus, the final positron has an energy E e = E ν − Q νep , with Q νep = 1.29 MeV, and it emits Cherenkov radiation visible in the detector.At typical SN energies, νµ and ντ are kinematically unable to interact via charged current (CC).
Above about 70 MeV, neutrino interactions in a water Cherenkov detector start to be dominated by CC reactions on oxygen of the form ν e + 16 O → e − + X, where X is a final excited nuclear state dominated by 16 F * [43-45] and a similar reaction for antineutrinos, where the dominant final state is 16 N * .The final state e ± retains memory of the initial neutrino energy.Specifically we use E e − = E ν − Q νeO , with Q νeO = 15.4MeV, and the positron energy is The cross sections are shown in Fig. S1, where the one for IBD is taken from Ref. [46], the one for νµ p scattering from Ref. [47], the ones for νe O and ν e O from Ref. [44], and the ones for νµ O and ν µ O from Ref. [48].
In this low-energy range, muon and tau neutrinos can only interact with nucleons via neutral-current interactions.In the interaction, nuclei can be excited and promptly decay to photons, leading to a potentially observable signature [49].For a future Galactic SN, this signature is likely to be observed.However, due to the lower cross sections of the neutral-current scattering, this process played no role for SN 1987A and we will not consider it even for our 100-MeV-range neutrinos.
At energies above the muon production threshold (m µ = 105.6MeV), the muon-flavored neutrinos from φ decay also contribute to the analogous CC rates.Due to large energy losses by ionization, these µ ± are stopped within a short length of the order of 1 m from their interaction vertex, and they finally decay at rest and produce a visible e ± .They follow the well-known Michel spectrum, behavior and detection efficiency.Leaving out this signal causes only a small and conservative error in the Majoron bounds (see main text).

B. SN 1987A Neutrino Observations
Supernova 1987A, in the Large Magellanic Cloud at a distance of 49.59 ± 0.09 stat ± 0.54 syst kpc from Earth [66], was discovered independently by Ian Shelton, Oscar Duhalde, and Albert Jones [50] on February 23, 1987, and later targeted by searches in the entire electromagnetic spectrum.The first evidence for optical brightening was found at 10:38 UT (Universal Time) on plates taken by McNaught.The first naked-eye visible (in the southern hemisphere) SN since the invention of the telescope, its observation is narrated in detail in a review by Koshiba [12].This was also the first SN explosion with a known progenitor star, Sanduleak −69 202, a blue supergiant catalogued by Nicholas Sanduleak in 1970 [51].At the time of the explosion, there were four running experiments that were big enough that they could have detected the gargantuan flux of neutrinos emitted in the collapse of a stellar core.
The largest one was the Irvine-Michigan-Brookhaven (IMB) water Cherenkov detector, an experiment built to look for proton decay [52], that was located in the Morton-Thiokol salt mine (Fairport, Ohio, USA).It was equipped with 2048 8-inch photomultiplier tubes (PMTs) such that 6,800 tons of water (of a total of 8,000 tons) were within the PMT planes, taken as the fiducial volume for the SN 1987A search [15].A failure of a high-voltage power supply shortly before SN 1987A left a contiguous quarter of the PMTs off-line with a geometric effect on the trigger efficiency that was later calibrated.The detector was triggered when at least 20 PMTs fired in 50 ns, corresponding to an energy threshold of 15-25 MeV for showering particles [15].(A trigger of 25 PMTs is mentioned in Ref. [14]).At the relatively shallow depth of 1570 m water equivalent, the flux of atmospheric muons caused a trigger rate of 2.7 Hz. Muons are recognized by tracks entering the detector from the outside and of course coming mostly from above.The detector is dead for 35 ms after each trigger.The SN 1987A signal consisted of 8 events and in addition 15 muons were recorded [16], a total of 23 triggers, amounting to 23 × 35 ms = 0.8 s dead time, or 13% of the SN signal duration of 6 s.In Fig. S3 we show the geometrically averaged detection efficiency, including the 0.87 reduction by dead time.
Atmospheric neutrinos are recognized as contained events and occurred at a rate of around 2/day in the energy range 20-2000 MeV [14].Our new neutrino signal in the 100 MeV range would look like low-energy atmospheric neutrinos.
The SN 1987A burst was found by looking in the recorded data for low-energy few-second event clusters, where "low energy" was defined as fewer than 100 PMTs firing, corresponding roughly to a 75 MeV energy cut.However, other than the 8 SN 1987A events and 15 muons, no other triggers occurred that would have been interpreted as a rare atmospheric neutrino.
The SN 1987A events must be due to IBD with a practically isotropic distribution of final-state e + .However, IMB found a conspicuous directional correlation in the opposite direction of SN 1987A, i.e., the events look "forward peaked."This effect is not explained by the detector's geometrical bias due to the 25% PMT failure.One idea held that the signal was not caused by neutrinos but instead some new X 0 bosons that scatter coherently on oxygen and thus generate the observed angular characteristic [54].However, the required cross section is excluded by stellar cooling bounds from the reverse process [55].No viable explanation other than a rare statistical fluctuation is available.
With a fiducial mass of 6.8 kton, IMB would have seen the largest number of 100 MeV-range events.At lower energies it suffered from a trigger efficiency of only 15% at 20 MeV, but rising to 80% at 70 MeV.During the SN 1987A burst, no events besides the 15 background muons + 8 SN events = 23 triggers were observed. 1We conclude that there were no unreported events above the low-energy criterion of 75 MeV.1987A neutrino data collected at Kamiokande-II, IMB, and Baksan.We show the detected positron energy as a function of time after the first event in each detector.Because of clock uncertainties, the exact temporal offset between the observations is not fixed.We do not show events which are attributed to background.
The second largest detector was the Kamiokande II water Cherenkov detector (Mozumi Mine, Kamioka section of Hida, Gifu Prefecture, Japan) with a fiducial mass of 2,140 metric tons for the SN 1987A search, where again the entire volume up to the PMTs was taken [9,10].This detector was built in 1983 to search for nucleon decay (Kamiokande = Kamioka Nucleon Decay Experiment) [56] and later upgraded to Kamiokande II to search for solar ν e in the 10 MeV range.The photo cathode coverage was increased and radioactive backgrounds decreased to lower the threshold and solar data were taken since the end of 1986.Despite its smaller mass, the low threshold made Kamiokande II competitive for the SN 1987A discovery (see Fig. S3 for the trigger efficiency) although for our 100-MeV-range events, IMB is better suited.
At a greater depth of 2700 m w.e., the atmospheric trigger rate was 0.37 Hz and indeed, 4 muons were found in the 20 s interval preceding the SN 1987A burst and several after the SN burst, but none until just before the 12th event.Atmospheric neutrinos, in the form of fully contained events, show up once every few days.Lowenergy radioactive backgrounds triggered with 0.23 Hz.The trigger dead time is less than 50 ns after an event.
To find the SN 1987A burst, the data recorded on a magnetic tape were searched for low-energy event clusters, where here the definition was less than 170 PMTs firing (E e < ∼ 50 MeV).We show the burst in Fig. S2 as a function of time after the first event.The absolute timing is poorly known, probably to within ±15 s based on comparing the computer clock with a wrist watch, but a conservative uncertainty of ±1 min was officially stated.A power outage in the mine on February 26 prevented a recalibration of the computer clock [57].The signal arrived at 4:35 pm on Monday, 23 February 1987, but this was a substitute holiday.According to working-day schedule, the magnetic tape would have been exchanged at 4:30 pm and the signal might have been missed.
The highest-energy events are also forward peaked, in analogy to IMB, while most of the events are isotropically distributed as expected for IBD.There is a conspicuous gap of 7.3 s between event 9 and 10, filled however with IMB data and probably has nothing to do with SN 1987A.Very recently, one member of the Kamiokande collaboration has speculated that the gap could have been caused by a fault of the magnetic tape drive.He noted that during that gap, there are also no other events (low-energy background or atmospheric muons) and that the probability for such a long gap was very small [13]. 2or our analysis, we are mainly interested in the highenergy events that Kamiokande would have seen during the SN 1987A burst.Contained events with 30 MeV < visible energy < 1.33 GeV would have gone into the atmospheric neutrino analysis, but none were found in the period around SN 1987A.For this analysis, the fiducial volume may have been as low as 780 tons (more than 2 m from the wall). 3We conclude that conservatively no event of our interest was observed in this volume.
The third experiment was the Baksan Scintillator Underground Telescope (BUST) under Mount Andyrchi in the North Caucasus at a depth of 850 m w.e., operated by the Institute of Nuclear Research (Moscow) [17,18].It started operation in June 1980 and is still running today, with SN 1987A the only SN neutrino burst observed in more than four decades.BUST consists of 3156 segments of 70 × 70 × 30 cm.A possible SN 1987A event was selected as one that triggers one and only one segment and with E e < ∼ 50 MeV.The fiducial inner part has a mass of 130 t that was opened for the SN 1987A analysis to 200 t.Its burst was reported at 7:36:06:571 UT and thus 30 s later than IMB.While the clock synchronization with UT is usually ±2 s, the clock was observed to have shifted forward by 54 s between February 17 and March 11 for unknown reasons.So the observed signal is probably contemporaneous with IMB and Kamiokande II.Because of its small size, BUST is least useful for us and so we have not investigated how our 100 MeV-range events would have shown up there.
A fourth instrument was the Liquid Scintillation Detector (LSD), located in the gallery of the Mont Blanc tunnel, between Italy and France [58,59].It was specifically built to search for a galactic SN burst with a typical assumed distance of 10 kpc.LSD used 72 100 × 150 × 100 cm 3 liquid scintillator modules, arranged in three horizontal layers for a total mass of 90 tons.Each module was equipped with three PMTs of 15 cm diameter, and the signal was recorded whenever a threefold coincidence occurred within 150 ns.
The LSD collaboration was the first to declare the (possible) discovery of SN neutrinos due to the detection of 5 events, above the 7 MeV threshold, in an interval of 7 seconds, beginning at UT 2:52:36.79 and compatible with the core-collapse standard model at 50 kpc.This signal is almost five hours earlier than the other detectors which observed nothing special at the LSD time and LSD observed nothing special the time of the others.While high multiplicity events can be caused e.g. by spallation of oxygen induced by primary muons, no similar event was found during the entire LSD operation which ended with the devastating fire in the Mont Blanc tunnel March 24, 1999.
The community has settled for the LSD event as being a rare or unexplained fluctuation.A credible physical origin at SN 1987A is astrophysically hard to construct.Schaeffer, Declais and Jullian computed that, assuming a SN origin for the events seen by LSD, the total energy emitted by SN 1987A would have been 3×10 54 erg, much larger than the value expected by standard core-collapse supernova theory [60].

C. Statistical Analysis
We perform our maximum likelihood analysis along the lines of similar previous studies [72,73].For the standard SN νe signal we assume a quasi-thermal distribution of the form Eq. ( 5) described by the three parameters E tot , E 0 and α.We compute the standard e + signal from the IBD cross section discussed in Sec.A and for the event spectrum in each detector use the efficiencies discussed earlier, including the IMB dead-time effect of 0.87.
The SN 1987A are not informative about α [73], so we do not try to fit it, but rather use a range of values motivated by numerical SN models.In particular, we use α = 2.39 (2.07) for the cold (hot) model.The instantaneous neutrino spectra are pinched, i.e., their variance is smaller than that of a Maxwell-Boltzmann spectrum (α > 2), whereas time-integrated spectra are close to Maxwell Boltzmann.The SN spectra somewhat depend on flavor, but the effect of flavor oscillations is not yet well understood and moreover, because of the LESA effect [61], the spectrum depends on the observer direction relative to the 3D structure of the SN explosion.
For each of the two experiments, we thus define an unbinned likelihood ) where E i are the observed energies, and E det is the energy reconstructed from the number of firing PMTs, drawn from a Poisson distribution as in Ref. [72].An unimportant normalization constant has been removed, because we will only deal with likelihood ratios.
In the event rate, we also include the new signal prediction that depends on the parameters g and m φ ; at small masses these appear in the combination gm φ and thus collapse to essentially a single parameter.For Kamiokande II, we reduce the fiducial volume from 2140 tons to 780 tons, as discussed above.We only consider the final-state e ± from CC reactions as well as from muon decay, but not the Cherenkov signal caused by muons above the Cherenkov threshold as discussed in the main text.We keep α fixed at the predicted value for the cold and hot SN model.We then marginalize over E 0 and E tot as explained in the main text.In this way, we obtain an effective two-dimensiona likelihood We now define a test statistic, The asymptotic distribution of this variable under the assumption that Majorons exist is a half-chi-squared distribution [74], which allows us to set a threshold value for 95% C.L. exclusion at χ 2 = 2.7.With this procedure, we find the limit contours shown in Fig. 1.

D. Garching Supernova Models
In our numerical analysis we use the SN models SFHo-18.8 and LS220-s20.0from the Garching group that were evolved with the Prometheus Vertex code with sixspecies neutrino transport [67] in spherical symmetry.These "muonic models" were recently also used for other particle constraints [24,29], where more details are described and radial profiles of various physical quantities are given for specific snapshots of time.PNS convection was taken into account by a mixing-length treatment.Explosions were triggered by hand a few 100 ms after bounce at the Fe/Si or Si/O composition interface of the progenitor star.
Following Ref. [29], we note that the SFHo equation of state is fully compatible with all current constraints from nuclear theory and experiment and astrophysics, including pulsar mass measurements and the radius constraints deduced from gravitational-wave and Neutron Star Interior Composition Explorer measurements.For comparison, some of the Garching muonic models also use the traditional LS220 equation of state.
The model SFHo-18.8 [29] uses a progenitor star with mass 18.8 M that reaches a final neutron-star baryonic mass of 1.351 M and gravitational mass of 1.241 M , hence a gravitational binding energy of (1.351 − 1.241)M = 0.110M = 1.98 × 10 53 erg.It is at the lower end of plausible neutron-star masses and released binding energy.It reaches a maximum core temperature near 40 MeV, the coldest of this suite of models.We thus refer to it as our "cold" model and it is taken to bracket the lower end of neutron-star mass and core temperature.
The "hot" model LS220-s20.0reaches a maximum core temperature of around 60 MeV.It has a progenitor mass of 20.0 M and reaches a neutron-star mass of 1.926 M , near the upper end of observed neutron-star masses.Its final gravitational mass is 1.707 M and thus releases 0.219M = 3.93 × 10 53 erg.This model is taken to bracket the upper end of both energy release and internal temperature.
In Figs.S4 and Figs.S5 we show several internal properties of these models as a function of time and mass coordinate for these two models.The left panels show the temperature and we see that after collapse the models are cold.They heat up at the edge of the inner core as they contract, with the maximum T and largest extent of the hot region achieved at around 1 s.Therefore, the emission rate of new particles would be largest around this time if the emission rate depends on temperature, as it often happens in other extensions of the Standard Model because it is the thermal energy of the medium constituents that is emitted.
However, in our case of Majoron emission by neutrino coalescence, the process ν e ν e → φ dominates by far, and so the chemical potential µ νe rather than T is the key quantity.It is shown in the middle panels, and we see that it is some 100 MeV up to roughly the inner 0.5 M , corresponding roughly to a radius of 10 km.At 1-2 s it drops quickly as the core deleptonizes.Beta equilibrium implies that ∆µ = µ e − µ νe = µ µ − µ νµ = µ n − µ p , whereas the number densities of ν e and e − must add up to the trapped lepton number of around 0.30 per baryon.However, the exact value of ∆µ depends on the nucleon properties in the medium and thus on the equation of state.Using free protons and neutrons provides the right order of magnitude, but is not a good approximation to estimate the emission rate, because in our case the latter scales rapidly as µ 3 νe (see main text).In these models with six-species neutrino transport, a chemical potential also builds up for ν µ in the sense that a significant population of νµ builds up, but the maximum of |µ νµ | remains a factor of 2-3 smaller than µ νe .As the emission rate scales with µ 3 ν , the muonic contribution remains only an order 10% correction.In Fig. S6, we finally show contours of the Majoron emission rate per unit mass.While the Majoron emission rate per unit volume scales as µ 3 ν and thus peaks at the center of the star, the emission rate per unit mass peaks at the edge of the inner core shown by the "yellow peak".This is because of the larger volume associated with the outer shells of the core.The chosen coupling strength for both "hot" and "cold" model coincides with the corresponding energy loss criterion detailed in the text, so that the Majoron luminosity at 1 s coincides with the neutrino luminosity.For the chosen coupling strength of g φ m φ = 8.3 × 10 −9 MeV (cold) and g φ m φ = 7.7 × 10 −9 MeV (hot), the emission rate is around 1 × 10 20 erg/gs throughout the inner core up to 0.50 M for the first second and then drops quickly.In the hot model, there is significant emission at larger mass coordinate around 0.5 s, deriving from the relatively large νµ population.

E. Neutrino chemical potentials and older models
Previous authors have derived SN 1987A energy-loss bounds, based on the same coalescence process, or have provided sensitivity forecasts for 100-MeV-range events from a future galactic SN [4,7].The emitting SN core was approximated as a one-zone model with µ νe = 200 MeV over a volume with R = 10 km and, in the case of Ref. [7], for a time scale of 10 s.These assumptions yield far more restrictive limits or far more ambitious signal predictions than our one-zone model or the numerical Garching models.
In Ref. [7], the chemical potential was taken from the pioneering paper [62] (see their Fig.11).In this protoneutron star (PNS) cooling simulation, the nuclear equation of state was still relatively rough.Moreover, the starting value of trapped lepton number per baryon of Y L = 0.35 was chosen as an initial condition and did not follow from a self-consistent SN simulation.More recent systematic PNS cooling simulations [63] used more sophisticated nuclear and microphysics, chose a similar initial Y L = 0.35, and found an initial value at the center of µ νe ∼ 170 MeV (see Fig. 9 for their baseline model).
Modern self-consistent simulations that include the infall phase find much smaller values of the trapped lepton number, 0.30 being a more typical number, depending on the progenitor model, also leading to smaller µ νe .In the muonic Garching models used here, the trapped lepton number in the center at core bounce is around 0.28 for the hot and 0.29 for the cold model and the initial µ νe ∼ 150 MeV at the center.
For Majoron emission, the geometrically largest region, very roughly around a mass coordinate of 0.5 M , is more relevant than the values at the center and so this region is indicative of the parameters that one could use for a one-zone description.This point is especially relevant for the time evolution because deleptonization occurs earlier at larger radii.Figure 11 of Ref. [62] reveals that after only a few seconds, µ νe strongly drops, and considering that the emission rate varies as µ 3 νe , the signal would strongly quench at 2-3 s and a similar conclusion follows from Fig. 9 of Ref. [63].
However, the deleptonization time scale can be much faster if the effect of PNS convection is included, in contrast to Refs.[62,63] or recently Ref. [64] who studied the late neutrino signal.We refer to a recent study of PNS evolution [65] (see this paper for references to the earlier literature) who found that convection, implemented with a mixing-length approximation, speeds up deleptonization by about a factor of 4 (see especially their Sec.4.1) and as such is crucial for determining the overall time scale.Of course, the exact quantitative impact on Majoron emission or on SN neutrino signal properties may not be captured by this single number which refers to deleptonization at the center of the star.
In our study we have used the numerical Garching models described earlier that include a mixing-length treatment of PNS convection, use nuclear equations of state that agree with modern information (notably on neutron-star masses and radii), and find trapped lepton abundances and chemical potentials commensurate with other modern simulations.
For the case of Majoron emission one can actually characterize the different models with a single figure, the trapped number of ν e in the core.In the degenerate limit, the Majoron luminosity of the SN core happens to be proportional to N νe , the total number of ν e present in the core as explained around Eq. (S15) below.For our cold model at core bounce, we find Nν = 3.5 in units (100 MeV) 3 (10 km) 3 , whereas at 1 s postbounce it is 0.74.These numbers justify the one-zone parameters adopted in the main text.
If instead one uses µ ν = 200 MeV with the same onezone radius 10 km, at 1 s one finds Nν = 8, about a factor of 11 larger and thus leading to much more restrictive energy-loss bounds as reported, for example, in Refs.[4,7].
For our argument about missing 100-MeV-range neutrinos in the 1987A data or the earlier forecasts for a future galactic SN [7], what matters is a somewhat different quantity.In the degenerate limit, the number emission rate scales with µ 2 ν .If we assume very roughly that the detection cross section scales with energy-squared, the count rate arising from a one-zone model scales with µ 4 ν R 3 τ as discussed in the main text.Therefore, one simple figure of merit for the source model is Ĉν = 3 dt dr r 2 µ 4 ν (r, t).Our cold model yields Ĉν = 2.30 in units of (100 MeV) 4 (10 km) 3 s.If one were to use a onezone model with µ ν = 200 MeV, R = 10 km and τ = 10 s, one instead finds Ĉν = 160, a factor of 70 larger than our value.
The sensitivity of the ν e abundance in the SN core, and its time-integrated value, to the microphysics input as well as deleptonization speed of the SN model mandates a somewhat careful gauging of one-zone parameters.

F. Majoron decay rate and emissivity
The matrix element for the decay of a single Majoron into a pair of neutrinos is This is also the matrix element for coalescence of a pair of neutrinos into a Majoron; notice that there are no additional factors coming from averages over spin states since we consider Majorana neutrinos.
The decay rate of a Majoron into a pair of neutrinos is where we denote by p 1 , p 2 , and p φ the four-momenta of the two neutrinos and the Majoron respectively, and in bold we denote their three-momenta.The factor 1/2 accounts for the presence of two identical particles in the final state.Performing the phase-space integral, we obtain In the case of neutrino coalescence, the rate of Majoron production from a pair of neutrinos, restricting to a single flavor, is (see, e.g., Ref. [1]) where f ν (E) is the neutrino phase-space distribution function.Performing the integral we recover as reported in the main text.Actually it is instructive to compare the rate of absorption Γ A (E φ ) of a Majoron in the neutrino background with the spontaneous rate of emission Γ E (E φ ).The rate of absorption is given by the vacuum decay rate Eq.(S7) times a Lorentz factor m φ /E φ .Moreover, the final-state neutrinos are Pauli-blocked so that overall we find (S10) where E ± = 1 2 (E φ ± p φ ) as defined in the main text.The integral expression is equal to 1 in the absence of Pauli blocking because the interval of integration has length p φ .
On the other hand, the emission rate per unit volume given in Eq. (S9) is equal to FIG. S1.Charged-current neutrino cross sections in a waterCherenkov detector.
FIG. S2.1987A neutrino data collected at Kamiokande-II, IMB, and Baksan.We show the detected positron energy as a function of time after the first event in each detector.Because of clock uncertainties, the exact temporal offset between the observations is not fixed.We do not show events which are attributed to background.
FIG. S3.Detection efficiencies for electrons and positrons at Kamiokande and IMB, taken from Ref.[72], including the dead-time effect in IMB.We continue them at energies higher than 60 MeV by extrapolation.
FIG. S4.Temperature (left), chemical potential of electron neutrinos (center), and chemical potential of muon neutrinos (right) as a function of post-bounce time and mass coordinate for the Garching "cold" model.The red line identifies the density 3 × 10 12 g cm −3 and thus essentially the edge of the PNS.The final neutron-star mass is 1.351 M .
The factor 1/6 represents assumed flavor equipartition.The parameters are chosen such that E tot , E 0 = E ν , and E 2 ν agree with the numerical spectrum.The cold model releases E tot = 1.98 × 10 53 erg.The exact impact of flavor oscillations on SN neutrinos is not yet fully understood.Averaging over all three ν flavors, we find E 0 = 12.7 MeV and α = 2.39.For the hot model, these parameters are E tot = 3.93 × 10 53 erg, E 0 = 14.3 MeV and α = 2.07.SN 1987A cooling limit.-Thelocal Majoron energy loss follows from Eq. ( 21.39 × 10 69 erg/s, to be compared with L ν (1 s) = 8.29 × 10 52 erg/s, leading to gm φ < 0.77 × 10 −8 MeV.Moreover, E tot φ = (gm MeV ) 2 4.39 × 10 69 erg and E tot ± from φ decay FIG.2.Normalized particle spectra from the time-integrated emission of the cold model SFHo-18.8."Standard ν" is the flavor average of the usual SN ν and "Standard e ± " the corresponding e ± spectrum in the detector (ignoring detection efficiencies), whereas the new contributions are marked "from φ decay."Theyinclude Michel e ± (endpoint 53 MeV) from µ ± decays at rest, which themselves emerge from CC interactions of νµ and νµ that come from φ decay.νe+ O → e + + X and ν e + O → e − + Y with X and Y excited final-state nuclei, dominate for E ν > ∼ 70 The absolute time of an event was recorded to an uncertainty ±50 ms thanks to the WWVB clock, a time signal radio station operated by the National Institute of Standards and Technology [53].The first IMB event occurred at 7:35:41.374Universal Time on 23 February 1987, corresponding to 2:35 am local time on a Monday very early morning.