#  CTPU-PTC-23-55 Limits on heavy neutral leptons, $Z^{\prime}$ bosons and
majorons
from high-energy supernova neutrinos

Kensuke Akita1***<EMAIL_ADDRESS>, Sang Hui<EMAIL_ADDRESS>,
Mehedi<EMAIL_ADDRESS>and Seokhoon<EMAIL_ADDRESS>
1Particle Theory and Cosmology Group, Center for Theoretical Physics of the
Universe,
Institute for Basic Science (IBS), Daejeon, 34126, Korea
2High Energy Physics Center, Chung-Ang University,
Dongjak-gu, Seoul 06974, Korea

( Abstract
Light hypothetical particles with masses up to $\mathcal{O}(100)\ {\rm MeV}$
can be produced in the core of supernovae. Their subsequent decays to
neutrinos can produce a flux component with higher energies than the standard
flux. We study the impact of heavy neutral leptons, $Z^{\prime}$ bosons, in
particular ${\rm U(1)}_{L_{\mu}-L_{\tau}}$ and ${\rm U(1)}_{B-L}$ gauge
bosons, and majorons coupled to neutrinos flavor-dependently. We obtain new
strong limits on these particles from no events of high-energy SN 1987A
neutrinos and their future sensitivities from observations of galactic
supernova neutrinos. )

###### Contents

  1. 1 Introduction
  2. 2 Reference SN core model
  3. 3 Neutrino oscillations in supernovae and secondary neutrino fluxes on Earth
  4. 4 Secondary neutrino fluxes by decays of light particles
    1. 4.1 General description
    2. 4.2 Heavy neutral lepton
    3. 4.3 ${\rm U(1)}_{L_{\mu}-L_{\tau}}$ gauge boson
    4. 4.4 ${\rm U(1)}_{B-L}$ gauge boson
    5. 4.5 Majoron
  5. 5 Limits from SN 1987A and a future galactic supernova
    1. 5.1 Event rates and statistical analysis
    2. 5.2 Heavy neutral lepton
    3. 5.3 ${\rm U(1)}_{L_{\mu}-L_{\tau}}$ gauge boson
    4. 5.4 ${\rm U(1)}_{B-L}$ gauge boson
    5. 5.5 Majoron
  6. 6 Summary

##  1 Introduction

The hot $(T\sim 30\ {\rm MeV})$ and dense $(\rho\sim 3\times 10^{14}\ {\rm
g/cm^{3}})$ cores in supernovae (SNe) provide an environment to test feebly
interacting particles (FIPs) [1]. In the core, particles with masses up to
$m\sim 100\ {\rm MeV}$ can be produced. If their interactions are considerably
weak, they constitute a channel of energy loss, shortening the duration of
neutrino burst from supernovae [2]. Using this argument, the observations of
SN 1987A neutrino burst constrain various FIP scenarios (e.g., axions and
axion-like particle [3, 4, 5, 6, 7, 8, 9, 10, 11], heavy neutral leptons
(HNLs) [12, 13], dark photons and gauge bosons [14, 5, 15, 9, 16, 17, 18],
majorons [19, 20, 21, 22, 23, 24, 25, 26, 27, 28]).

Besides the search by the energy loss argument, one can search for FIPs by
studying secondary peculiar photon and neutrino fluxes due to their secondary
interactions. If FIPs with a sufficiently long lifetime escape from the SN
core and decay into photons or neutrinos, the produced particles are not
thermalized in the SN medium, giving a secondary flux. The typical energy of
the secondary neutrino flux is $\gtrsim 100$ MeV, corresponding to the FIP
energy related to the core temperature (or their mass, etc.), $E\gtrsim 3T$.
In addition, if FIPs decay into electron-positron pairs, the annihilation of
the produced positron with electrons in the Galaxy can contribute a galactic
$511\ {\rm keV}$ $\gamma$-ray signal. The search for the high-energy secondary
fluxes from SNe in $\gamma$-ray and neutrino telescopes gives improved
constraints on axion-like particles [29, 30, 31, 32, 33, 34, 35, 36], dark
photons [14, 37, 38, 17], majorons [39, 40], and HNLs [38, 12, 41]. One can
also search for light particles by using the excessive energy deposition in
the stellar envelope by FIP decays over the observed energy in Type IIP SNe
[42, 43, 38, 10, 44, 17, 45] and the change of duration of SN 1987A neutrino
burst [46, 47, 48].

In this work, building upon the previous research [12, 39, 40, 41], we derive
new supernova constraints on light particles interacting with neutrinos by
estimating observable secondary neutrino fluxes due to their decay. As an
illustrative example, we explore the implications for heavy neutral leptons
(see e.g., ref. [49] for review and motivation of HNLs), $Z^{\prime}$ boson,
particularly ${\rm U(1)}_{L_{\mu}-L_{\tau}}$ [50, 51, 52] and ${\rm
U(1)}_{B-L}$ [53, 54, 55, 56] gauge bosons. Additionally, we examine
(pseudo-)scalar bosons coupled to neutrinos in a flavor-dependent way,
referred to as flavored majorons [57, 58, 59]. These particles can potentially
address some experimental anomalies [60, 61, 62, 63, 64, 65, 66, 67, 68, 69,
70, 71] and/or cosmological tensions [49, 72, 73, 74, 75, 76, 77, 78, 79, 80].
Secondary neutrino fluxes are produced by their decay in the energy region of
$\mathcal{O}(100)$ MeV, higher than the standard neutrino burst with energy
$E_{\nu}\simeq 15\ {\rm MeV}$ which emerges in the neutrino decoupling region
outside the core (called the neutrino sphere). Furthermore, the detectability
of SN neutrinos in neutrino telescopes is more sensitive to higher energies
due to the scattering cross section, which scales roughly as
$\sigma_{\nu}\propto E_{\nu}^{2}$ [81, 82, 83, 84]. The search for such high-
energy SN neutrinos would significantly improve the supernova limits on these
particles. We derive limits on these light particles from the SN 1987A
observations in Kamiokande-II [85] and Irvine-Michigan-Brookhaven (IMB) [86]
water Cherenkov detectors. The SN 1987A neutrinos are also observed in Baksan
Underground Scintillation Telescope (BUST) [87] but we do not study it because
of its small volume. We also discuss the projected sensitivities from future
observations of a galactic supernova, which occurs a few times per century
[88], using Hyper-Kamiokande (HK) [89].

We comment on the differences between the previous works [12, 39, 40, 41] and
this work. While refs. [39, 40] took into account secondary neutrino fluxes
through flavor-universal decays of light bosons, we investigate a broader
spectrum of particle physics models and flavor-dependent decays of light
bosons to neutrinos. Our analysis extends to estimating the secondary fluxes
at neutrino detectors, accounting for the effects of neutrino oscillations
during their journey to Earth. Refs. [12, 41] focused on secondary fluxes
resulting from decays of HNLs and the future sensitivity from observations of
galactic supernovae. In this work, we improve various calculations, including
the production rate of HNLs in the SN core, neutrino oscillations, and the
event rates at the detector. Furthermore, we obtain the current limits on HNLs
from SN 1987A observations, building upon the above argument.

The paper is organized as follows: In section 2, we present our benchmark
model of the SN core. In section 3, we discuss the effects of neutrino
oscillations on the secondary neutrino fluxes by light particle decays. In
section 4, we briefly characterize the secondary neutrino flux and fluence
(time-integrated flux) by the light particle decays. In section 5, we show the
new limits on HNLs, $Z^{\prime}$ bosons, and flavored majorons, including the
event rates in neutrino telescopes and the methods of our statistical
analysis. Section 6 summarizes the paper.

##  2 Reference SN core model

For our estimations of light particle emissions from the collapsed core, we
use the SN core model SFHo-18.8 [90, 8, 10, 40] in spherical symmetry from the
Garching group. For the actual values of physical quantities (temperature, the
number densities, and the chemical potentials of the Standard Model (SM)
particles) in the SN core, we take them from figures 1 and 5 in ref. [10] and
figure S5 in ref. [40]. For the proton and electron number density $n_{p}$, we
always take $Y_{p}\simeq Y_{e}\equiv n_{p,e}/n_{b}\simeq 0.2$ with the charge
neutrality condition [10, 91], where $n_{b}$ is the baryon number density. The
number density for pions is computed using the relativistic virial expansion
and other physical quantities in SFHo-18.8, based on ref. [92]. The SN model
SFHo-18.8 is a conservative model, where the final neutron star mass is at the
lower edge of the allowed range for SN 1987A [93]. We neglect feedback onto
the dynamics of the collapsed core due to the emission of hypothesis light
particles. This assumption would be correct, at least when the luminosity for
light particles is smaller than that for neutrinos. However, even when the
luminosity for light particles is comparable with or larger than that for
neutrinos, we neglect the feedback because of lacking SN simulations including
the light particles, as the rest of the literature does.

In a dense medium, nucleon masses are reduced due to nucleon forces [94]. We
take the ratio of effective nucleon mass with respect to the vacuum one as a
function of the mass density, based on the relativistic mean-field treatment,
from ref. [6]. In SN physics, the mass density is usually defined as
$\rho\equiv m_{u}n_{b}$, where $n_{b}$ is the baryon number density and
$m_{u}=931.494\ {\rm MeV}$ is the atomic mass. From the nucleon number
densities and effective nucleon masses, we calculate the effective nucleon
chemical potentials.

We summarize the typical values of physical quantities in the core required
for our calculations. Inside the core radius of $r\sim 10\ {\rm km}$, where
the SM particles are in thermal equilibrium, the temperature is $T\sim 30\
{\rm MeV}$ for the duration of $t=5\ {\rm s}$ after the core bounce. After
$t=5\ {\rm s}$, the temperature decreases due to neutrino emissions. Inside
the core, the baryon density is $\rho\sim 3\times 10^{14}\ {\rm g/cm^{3}}$,
corresponding to the nuclear saturation density, with the effective nucleon
masses of $m_{n,p}\sim\mathcal{O}(500)\ {\rm MeV}$ and the effective chemical
potential of $\mu_{n}\sim\mathcal{O}(600)\ {\rm MeV}$ and
$\mu_{p}\sim\mathcal{O}(500)\ {\rm MeV}$ in a neutron-rich star. Such a
sizable difference $\mu_{n}-\mu_{p}\sim\mathcal{O}(100)\,{\rm MeV}$ implies
that negatively charged particles such as electrons can be largely populated
in the core. Muons are indeed populated in the core because highly degenerated
electrons with $\mu_{e}\sim\mu_{n}-\mu_{p}\gtrsim m_{\mu}$, thermally
distributed photons, and neutrinos with $T\sim 30\ {\rm MeV}$ can easily be
converted into muons. Radiating a slight excess of $\nu_{\mu}$ compared to
$\bar{\nu}_{\mu}$ create a net muon number [95], corresponding to
$Y_{\mu}\equiv n_{\mu}/n_{b}\sim 0.02$ with the muon number density $n_{\mu}$.
Negatively charged pions $\pi^{-}$ would also be populated in thermal and
chemical equilibrium due to strong interactions with $Y_{\pi^{-}}\equiv
n_{\pi^{-}}/n_{b}\sim 0.01$ [92, 91]. The chemical potentials for $\nu_{e}$
and $\bar{\nu}_{e}$ are $\mu_{\nu_{e}}=-\mu_{\bar{\nu}_{e}}\sim 100\ {\rm
MeV}$ for duration of $2\ {\rm s}$ [40] while those for the other flavors are
approximately zeros. After 2 s, the chemical potential for $\nu_{e}$ becomes
also zero by the depleptonization due to a large emission of $\nu_{e}$.

The standard SN neutrino fluxes follow a quasi-thermal spectrum [96, 97]. We
write the time-integrated SN neutrino flux as

$\displaystyle\frac{dN_{\nu}}{dE_{\nu}}=\frac{E_{\rm
tot}}{6E_{0}^{2}}\frac{(1+\alpha)^{1+\alpha}}{\Gamma(1+\alpha)}\left(\frac{E_{\nu}}{E_{0}}\right)^{\alpha}e^{-(1+\alpha)E_{\nu}/E_{0}},$
(1)

where $E_{\rm tot}$ is the total energy emitted from a supernova, $E_{0}$ is
the average energy for a neutrino species and $\alpha$ is a ”pinching”
parameter that is 2 for the Maxwell-Boltzmann distribution and not 2 when the
spectrum is pinched from the distribution. $\Gamma$ is the Gamma function.
$E_{\rm tot}/6$ is the total emitted energy for one neutrino species. In the
SN model SFHo-18.8, the total energy is $E_{\rm tot}=1.98\times 10^{53}\ {\rm
erg}$ corresponding to the released gravitational binding energy. The exact
effects of neutrino oscillations in the SN core are still under investigation.
After averaging over all three $\bar{\nu}$, whose spectra produced in the core
are almost the same [97], the other parameters are $E_{0}=12.7\ {\rm MeV}$ and
$\alpha=2.39$ [40].

##  3 Neutrino oscillations in supernovae and secondary neutrino fluxes on
Earth

To estimate the SN neutrino spectrum on Earth, we have to consider the effects
of neutrino oscillations en route to Earth. In the medium of the core,
electron neutrinos are produced through both the neutral and charged current
interactions while non-electron neutrinos
($\nu_{x}=\nu_{\mu},\nu_{\tau},\bar{\nu}_{\mu},\bar{\nu}_{\tau}$) are produced
through the approximately same neutral current interactions. Then the standard
initial SN fluxes for $\nu_{x}$ would be equal. However, in general, the
secondary neutrino fluxes produced by the decays of light particles are
different even for each non-electron neutrino flavor. In this section, we
extend the conventional discussion of neutrino oscillations in supernovae to
its application to the secondary neutrino fluxes, following ref. [98].

Effects of neutrino oscillations on the fluxes depend on the position of light
particle decays and the density of the SN medium. In the very high-density
region, neutrino self-interactions are not negligible. Neutrino oscillations
by the self-interactions are still under investigation [99, 100, 101, 102].
However, the secondary neutrino fluxes distinguishable from the standard one
are produced outside the neutrino sphere, where neutrino self-interactions
decouple. Then we will neglect neutrino oscillations by neutrino self-
interaction in the high-density region.

Neutrinos also interact with the electron background in the SN envelope, which
induces neutrino oscillations in a matter (see e.g., refs. [103, 104] for a
review). For a constant electron background, neutrino states in the matter
basis, which diagonalize the Hamiltonian including neutrino masses and the
electron potential, do not oscillate. For varying electron density, the matter
eigenstates change while neutrinos propagate in matter. Then the eigenstates
in the matter basis with a fixed density oscillate. Even for varying density,
if the change is slow (i.e., adiabatic), the component of the prompt
Hamiltonian eigenstates will retain that of the respective initial states
during their propagation [103]. This approximation is called the adiabatic
approximation, which is valid for [105]

$\displaystyle\gamma\equiv\frac{\Delta
m^{2}}{2E_{\nu}}\frac{\sin^{2}2\theta}{\cos
2\theta}\frac{1}{(1/n_{e})(dn_{e}/dr)}\gg 1,$ (2)

where $\Delta m^{2}$ and $\theta\in[0,\pi/2]$ are the neutrino mass squared
differences and the mixing angle for flavor neutrinos in vacuum, respectively,
and $n_{e}$ is the electron number density. Related with the electron
potential, there are two parameter combinations, $(\Delta
m^{2},\theta)=(\Delta m_{31}^{2},\theta_{13})$ and $(\Delta
m_{21}^{2},\theta_{12})$ with $\Delta m_{ij}^{2}=m_{i}^{2}-m_{j}^{2}$. In the
region of SN with densities of $\gtrsim 1\ {\rm g/cm^{3}}$, the electron
fraction in the SN envelope is $Y_{e}\equiv n_{e}/n_{b}=n_{e}m_{u}/\rho\simeq
0.5$ [98, 100]. We assume the density profile in the SN envelope approximately
as [106, 107]

$\displaystyle\rho Y_{e}\sim 10^{5}\ {\rm g/cm^{3}}\left(\frac{r}{10^{9}\ {\rm
cm}}\right)^{-3}.$ (3)

In the SN envelope, the adiabatic approximation would be valid well for
neutrinos with energy of $E_{\nu}\lesssim 500\ {\rm MeV}$ [98]. For light
particles decaying to two neutrinos, this approximation would be valid for
their mass of $m\lesssim 1\ {\rm GeV}$, which is enough larger than masses
that can be produced in the core. In the following, we adopt the adiabatic
approximation to the effects of neutrino oscillations on the secondary fluxes.

The oscillation probability of flavor neutrinos from a production point of
neutrinos with a given electron density $n_{e}$ to Earth is given by,
averaging over time due to their long propagation,

$\displaystyle P_{\nu_{\alpha}\rightarrow\nu_{\beta}}=\sum_{i}\left|U_{\alpha
i}^{M}\right|^{2}\left|U_{\beta i}\right|^{2},$ (4)

and the same expression holds for antineutrinos. $U_{\beta i}$ is the
Pontecorvo–Maki–Nakagawa–Sakata (PMNS) matrix and $U^{M}_{\alpha i}$ is the
effective PMNS matrix at the production point approximately given by [104]

$\displaystyle U^{M}$ $\displaystyle=\begin{pmatrix}1&0&0\\\
0&c_{23}&s_{23}\\\ 0&-s_{23}&c_{23}\end{pmatrix}\begin{pmatrix}1&0&0\\\
0&1&0\\\ 0&0&e^{i\delta_{\rm
CP}}\end{pmatrix}\begin{pmatrix}c^{M}_{13}&0&s^{M}_{13}\\\ 0&1&0\\\
-s^{M}_{13}&0&c^{M}_{13}\end{pmatrix}\begin{pmatrix}c_{12}^{M}&s_{12}^{M}&0\\\
-s_{12}^{M}&c_{12}^{M}&0\\\ 0&0&1\end{pmatrix},$ (5)

where $c_{ij}=\cos\theta_{ij}$, $s_{ij}=\sin\theta_{ij}$, $\theta_{ij}$ and
$\delta_{\rm CP}$ are the mixing angle and CP-violating phase in vacuum,
respectively. $c^{M}_{ij}=\cos\theta^{M}_{ij}$ and
$s^{M}_{ij}=\sin\theta_{ij}^{M}$ satisfy approximately [104],

$\displaystyle\tan 2\theta^{M}_{13}$ $\displaystyle=\frac{\tan
2\theta_{13}}{1-\frac{2E_{\nu}V_{e}}{\Delta m_{31}^{2}\cos 2\theta_{13}}},$
(6) $\displaystyle\tan 2\theta^{M}_{12}$ $\displaystyle=\frac{\tan
2\theta_{12}\cos\theta^{M}_{13}}{1-h},$ (7)

with

$\displaystyle h=\left\\{\begin{array}[]{ll}\frac{\Delta m_{31}^{2}}{\Delta
m_{21}^{2}}-\frac{\sin^{2}\theta_{12}}{\cos 2\theta_{12}}&{\rm for}\ \
2E_{\nu}V_{e}/\Delta m_{31}^{2}\cos 2\theta_{13}>1\\\
\frac{2E_{\nu}V_{e}}{\Delta_{21}^{2}\cos 2\theta_{12}}&{\rm for}\ \
2E_{\nu}V_{e}/\Delta m_{31}^{2}\cos 2\theta_{13}<1\end{array}\right.$ (10)

Here $\Delta m_{ij}\equiv m_{i}^{2}-m_{j}^{2}$, $V_{e}=\sqrt{2}G_{F}n_{e}$ for
neutrinos, $V_{e}=-\sqrt{2}G_{F}n_{e}$ for antineutrinos and $G_{F}$ is the
Fermi constant.

The secondary neutrino fluxes produced by light particle decays on Earth are
given by

$\displaystyle
F_{\nu_{\beta}}=\sum_{\alpha=e,\mu,\tau}P_{\nu_{\alpha}\rightarrow\nu_{\beta}}F_{\nu_{\alpha}}^{0},$
(11)

where $F_{\nu_{\alpha}}^{0}$ is the initial neutrino flux produced by their
decays. To reduce computational complexity, we approximate that light particle
decays occur instantaneously when we estimate this effect,

$\displaystyle
P_{\nu_{\alpha}\rightarrow\nu_{\beta}}=P_{\nu_{\alpha}\rightarrow\nu_{\beta}}|_{r=r_{\rm
D}},\ \ \ \ r_{\rm D}=\tau\gamma\beta,$ (12)

where $\tau$, $\gamma=E/m$, and $\beta=p/E$ are the lifetime, Lorentz boost
factor, and velocity for the light particle, respectively. For
$\rho\rightarrow 0$, $P_{\nu_{\alpha}\rightarrow\nu_{\beta}}$ in eq. (4)
reduces the well-known time-averaged oscillation probability in vacuum.

In figure 1, we show the oscillation probabilities, eq. (12), as a function of
the matter density, $\rho$, with neutrino energy of $E_{\nu}=10\ {\rm MeV}$.
We only consider $P_{\nu_{\alpha}\rightarrow\nu_{e}}$ and
$P_{\bar{\nu}_{\alpha}\rightarrow\bar{\nu}_{e}}$ $(\alpha=e,\mu,\tau)$ because
water Cherenkov detectors mainly detect $\nu_{e}$ and $\bar{\nu}_{e}$ in the
10-100 MeV-energy region. For neutrino mass-squared differences and mixing
parameters in vacuum, we use their best-fit values [108] 111NuFIT 5.2 (2022),
www.nu-fit.org (see also ref. [109]). The left figure denotes the case of the
normal ordering (NO) of neutrino masses, $\Delta m^{2}_{31}>0$, while the
right figure denotes the inverted ordering (IO) of neutrino masses, $\Delta
m^{2}_{31}<0$. From this figure, we can expect the following general results
of neutrino oscillations on the secondary flux:

  * •

$P_{\nu_{\mu}\rightarrow\nu_{e}}$ and $P_{\nu_{\tau}\rightarrow\nu_{e}}$ are
very similar as well as $P_{\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{e}}$ and
$P_{\bar{\nu}_{\tau}\rightarrow\bar{\nu}_{e}}$. This indicates that the ratio
between $\nu_{\mu}$ and $\nu_{\tau}$, and their anti particles, produced by
light particle decays in the SN envelope does not affect the signal at Earth
approximately. Collectively, we would still denote $\nu_{\mu}$ and
$\nu_{\tau}$ as $\nu_{x}$ ($x=\mu,\tau$) as in the standard SN neutrinos.

  * •

The difference between the total signal at Earth in the NO and IO case would
not very significant because water Cherenkov detectors detect both $\nu_{e}$
and $\bar{\nu}_{e}$ (see e.g., FIG. S1 in ref. [40] for the detection cross
sections in a water Cherenkov detector). One exception would be the case light
particles decay into only $\nu_{e}$ or $\bar{\nu}_{e}$ in the high density
region.

Figure 1: The oscillation probabilities for $\nu_{\alpha}\rightarrow\nu_{e}$
and $\bar{\nu}_{\alpha}\rightarrow\bar{\nu}_{e}$,
$P_{\nu_{\alpha}\rightarrow\nu_{e}}$, ($\alpha=e,\mu,\tau$) as a function of
the matter density $\rho$ with $E_{\nu}$ in unit of $10\ {\rm MeV}$. The
vertical solid and dashed lines are the two resonance values,
$2E_{\nu}V_{e}=\Delta_{31}^{2}\cos 2\theta_{13}$ and $2E_{\nu}V_{e}=\Delta
m_{21}^{2}\cos 2\theta_{12}$, with $V_{e}=2\sqrt{2}G_{F}n_{e}$ ($n_{e}$:
electron density). Left: Normal Ordering of neutrino masses, $\Delta
m^{2}_{31}>0$. Right: Inverted Ordering of neutrino masses, $\Delta
m^{2}_{31}<0$.

##  4 Secondary neutrino fluxes by decays of light particles

###  4.1 General description

First, we introduce the secondary neutrino flux produced by the decays of
hypothetical light particles generally and schematically. We denote light
particles as $\chi$. To calculate the light particle spectrum produced in the
core, we solve the Boltzmann equation governing their production in region
small enough to be constant density and large enough for particle interactions
[110, 111],

$\displaystyle\frac{\partial f_{\chi}}{\partial t}=\mathcal{C}[f],$ (13)

where $f_{\chi}$ is the distribution of light particles and $\mathcal{C}$ is
the collision term summed over all possible collisional interactions. The
emitted spectrum of light particles is given by,

$\displaystyle\frac{dN_{\chi}}{dE_{\chi}}=\int dt\int
dVe^{-\Gamma_{\chi}r_{\nu}/(\gamma\beta)}\int
d\Omega\frac{E_{\chi}p_{\chi}}{(2\pi)^{3}}f_{\chi}=\int dt\int
dVe^{-\Gamma_{\chi}r_{\nu}/(\gamma\beta)}\frac{E_{\chi}p_{\chi}}{2\pi^{2}}f_{\chi},$
(14)

where $t$ and $V$ denote the emitting time and the core volume, respectively,
and $\Gamma_{\chi}$ is the total decay rate of light particles in the rest
frame. When light particles decay inside the neutrino sphere of $r_{\nu}\sim
20\ {\rm km}$ [112] (outside of which neutrinos decouple with the medium), the
secondary neutrino flux would be thermalized, being the same with the standard
flux 222In general, neutrinos with higher energy decouple from the medium at
lower medium density. Around the standard neutrino sphere, $r_{\nu}\sim 20\
{\rm km}$, the matter density rapidly decreases [112]. The neutrino sphere for
high energy neutrinos, which is produced by the light particle decays, would
be the same with the standard one.. We account for the secondary flux produced
only outside the neutrino sphere, plugging the exponential factor
$e^{-\Gamma_{\chi}r_{\nu}/(\gamma\beta)}$.

If light particles interact with the SN medium feebly, they freely escape from
the SN (called the free-streaming case). On the other hand, if they interact
with the SN medium strongly, they are absorbed and/or rescattered (called the
trapping case). In the trapping case, light particles would be less likely to
both carry energy out of the SN and produce the secondary neutrino flux. This
effect for their decays is already included as
$e^{-\Gamma_{\chi}r_{\nu}/\gamma\beta}$ in eq. (14). The effects for
absorption and rescattering are characterized by the optical depth, which is
given by [9, 10],

$\displaystyle\tau_{\rm
abs}(E,r)=\int_{r}^{\infty}dr^{\prime}\Gamma(E,r^{\prime})/\beta,$ (15)

where $\Gamma$ is the effective absorption (or rescattering) rate and $r$ is
the production point for light particles. The effective absorption rate is
defined by $\Gamma\equiv\Gamma_{\rm abs}\mp\Gamma_{\rm prod}$ for bosonic
(fermionic) $\chi$ with $-\ (+)$ sign. $\Gamma_{\rm abs}$ and $\Gamma_{\rm
prod}$ are the absorption and production rates, respectively, which are
defined as $\mathcal{C}[f]=\Gamma_{\rm prod}(1\pm f_{\chi})-\Gamma_{\rm
abs}f_{\chi}-\Gamma_{\chi}f_{\chi}$ in eq. (13). If the medium is in thermal
equilibrium, the absorption and production rates are related as $\Gamma_{\rm
prod}=\Gamma_{\rm abs}e^{-E_{\chi}/T}$ by the principle of the detailed
balance. We use this relation when evaluating eq. (15). The emitted spectrum
from SNe is given by

$\displaystyle\frac{dN_{\chi}}{dE_{\chi}}=\int dt\int dV\langle e^{-\tau_{\rm
abs}(E_{\chi},r)}\rangle
e^{-\Gamma_{\chi}r_{\nu}/(\gamma\beta)}\frac{E_{\chi}p_{\chi}}{2\pi^{2}}f_{\chi},$
(16)

where $\langle e^{-\tau_{\rm abs}(E_{\chi},r)}\rangle$ is a directional
average of the absorption factor [10],

$\displaystyle\langle e^{-\tau_{\rm
abs}(E_{\chi},r)}\rangle=\frac{1}{2}\int^{1}_{-1}d\mu\
e^{-\int^{\infty}_{0}ds\ \tau_{\rm abs}(E,\sqrt{r^{2}+s^{2}-2rs\mu})}.$ (17)

Typically, the trapping case is already excluded by other cosmological
observations and experiments. In addition, due to the exponential suppression,
we would not improve significantly the constraints on the trapping regime from
the SN energy loss argument.

Light particles will also be trapped by the gravitational potential of the
core when the kinetic energy of $\chi$ satisfies [113, 12]

$\displaystyle E_{\rm kin}\leq K_{\rm tr}\equiv\frac{G_{N}M_{r}m_{\chi}}{r},$
(18)

where $G_{N}$ is the Newton constant, $M_{r}$ is the enclosed mass of the core
inside a radius $r$. One can take into account this effect by modifying the
energy spectrum for light particles,

$\displaystyle\frac{d^{2}N_{\chi}}{dE_{\chi}dr}\rightarrow\frac{d^{2}N_{\chi}}{dE_{\chi}dr}\Theta(E_{\chi}-K_{\rm
tr}-m_{\chi}),$ (19)

where $\Theta(x)$ is the Heaviside step function.

All emitted $\chi$ with masses and couplings of our interest decay en route to
Earth. The emitted number of $\nu_{\alpha}$ (per volume in the SN core) per
$t_{D}$ and $\omega_{\nu}$, where $t_{D}$ and $\omega_{\nu}$ are the time at
the $\chi$ decay in the laboratory frame and the neutrino energy in the rest
frame of $\chi$, is given by [114, 12, 39]

$\displaystyle\frac{d^{2}N_{\nu_{\alpha}}}{dt_{D}d\omega_{\nu}}={\rm
Br}_{\nu_{\alpha}}\bar{N}_{\nu_{\alpha}}\int\cos\theta\int
dE_{\chi}\frac{1}{\tau\gamma}\exp\left(-\frac{t_{D}}{\tau\gamma}\right)\frac{dN_{\chi}}{dE_{\chi}}(t_{D}=0)f(\omega_{\nu},\cos\theta),$
(20)

where $\theta$ is the emitted angle of $\nu_{\alpha}$ in the rest frame of
$\chi$, $\tau=\Gamma_{\chi}^{-1}$ the lifetime of $\chi$ in the rest frame,
${\rm Br}_{\nu_{\alpha}}$ denotes the branching ratio for a decay process, and
$\bar{N}_{\nu_{\alpha}}$ is the total number of the emitted $\nu_{\alpha}$ by
a $\chi$ decay. ${\rm
Br}_{\nu_{\alpha}}\bar{N}_{\nu}/(\tau\gamma)\exp(-t_{D}/(\tau\gamma))$ is the
production rate for $\nu_{\alpha}$ per time, i.e.,
$(dN_{\chi}/dt_{D})/N_{\chi}(t_{D}=0)$. $dN_{\chi}/dE_{\chi}(t_{D}=0)$ is
given by eq. (14). $f(\omega_{\nu},\cos\theta)$ is a distribution function for
$\omega_{\nu}$ and $\theta$ by a $\chi$ decay, normalized to be 1.

Neutrinos by the $\chi$ decays are emitted at an angle relative to the $\chi$
momentum. Then they travel a distance until they arrive at Earth, which is
longer than the travel distance of the standard SN neutrinos. The time delay
between neutrinos by the $\chi$ decays and the standard ones arriving at Earth
is [12]

$\displaystyle t_{\rm delay}=\frac{t_{D}}{\gamma^{2}(1+\beta\cos\theta)},$
(21)

where the time delay is suppressed by the Lorentz boost factor, $\gamma$. If
$t_{\rm delay}$ is longer than the data-taking time in neutrino telescopes,
the neutrinos by the light particle decays cannot be observed. We take into
account this by modifying

$\displaystyle\frac{d^{3}N_{\nu_{\alpha}}}{dt_{D}d\omega_{\nu}d\cos\theta}\rightarrow\frac{d^{3}N_{\nu_{\alpha}}}{dt_{D}d\omega_{\nu}d\cos\theta}\Theta(t_{\rm
data}-t_{\rm delay}),$ (22)

where $t_{\rm data}$ is the data-taking time after the first SN neutrino
arrives at Earth. We find the secondary neutrino fluence (time-integrated
flux) with energy in the laboratory frame,
$E_{\nu}=\gamma(1+\beta\cos\theta)\omega_{\nu}$,

$\displaystyle\frac{dN_{\nu_{\alpha}}}{dE_{\nu}}$ $\displaystyle={\rm
Br}_{\nu_{\alpha}}\bar{N}_{\nu_{\alpha}}\int^{1}_{\cos\theta_{\rm
min}}\cos\theta\int^{\infty}_{E_{\rm
min}}dE_{\chi}\frac{1}{\gamma(1+\beta\cos\theta)}\left[1-\exp\left(-\frac{\gamma(1+\beta\cos\theta)t_{\rm
data}}{\tau}\right)\right]\frac{dN_{\chi}}{dE_{\chi}}$ $\displaystyle\ \ \ \ \
\ \ \ \times
f\left(\frac{E_{\nu}}{\gamma(1+\beta\cos\theta)},\cos\theta\right),$ (23)

where $E_{\rm min}$ is the minimal energy of $\chi$ to produce neutrinos with
$E_{\nu}$ and $\cos\theta_{\rm min}$ is determined by the kinematically
allowed region for $\nu$ in the rest frame of $\chi$.

###  4.2 Heavy neutral lepton

For heavy neutral leptons mixing with $\nu_{\mu,\tau}$, the range of mixing
angle, $10^{-9}\lesssim U_{\mu N,\tau N}\lesssim 10^{-4}$, for masses of
$\mathcal{O}(100)\ {\rm MeV}$ is currently unconstrained. On the other hand,
this parameter region is severely constrained for HNLs mixing with $\nu_{e}$
by the terrestrial experiments [115, 116, 117, 118, 119] and Big Bang
Nucleosynthesis (BBN) [120, 121, 122]. We consider a single mass eigenstate of
heavy neutral leptons, $N$, mixing with $\nu_{\mu,\tau}$, with mass of $m_{N}$
and mixing angle $U_{\mu N,\tau N}$.

The Boltzmann equation for the production of HNLs is given by [123]

$\displaystyle\frac{\partial f_{N}(p_{N})}{\partial t}=\mathcal{C}[f]$ (24)

with

$\displaystyle\mathcal{C}[f]=\frac{1}{2E_{N}}\int\prod_{i=1,2,3}\frac{d^{3}p_{i}}{(2\pi)^{3}2E_{i}}S|\mathcal{M}|^{2}(2\pi)^{4}\delta^{(4)}(p_{3}+p_{N}-p_{1}-p_{2})f_{1}f_{2}(1-f_{3})(1-f_{N}),$
(25)

where $f_{i}$ is the thermalized distribution function, which is the Fermi-
Dirac distribution $f_{i}=[e^{(E_{i}-\mu_{i})/T}+1]^{-1}$ in the core, $f_{N}$
is the distribution function of HNLs, $S$ is the symmetry factor,
$|\mathcal{M}|^{2}$ is the squared matrix element summed over all spins except
for HNLs. In the core, $e^{-}$ are highly degenerate while
$n,p,\mu^{-},\nu_{e},\bar{\nu}_{e}$ are partially degenerate [8, 10, 40, 13].
Pauli-blocking suppresses the production of HNLs via the weak interactions
with $e^{-}$ while the Pauli-blocking effects for
$n,p,\mu^{-},\nu_{e},\bar{\nu}_{e}$ are sub-dominant. The dominant production
processes of HNLs in the core are neutrino scattering with nucleons via
neutral current interactions, $\overset{(-)}{\nu}n\rightarrow nN$ and
$\overset{(-)}{\nu}p\rightarrow pN$ [13]. For HNLs mixing with $\nu_{\mu}$,
the other dominant production processes are muon scattering with nucleons via
charged current interactions, $\mu^{-}p\rightarrow Nn$. Their squared matrix
elements are listed in table 1. Since the produced HNLs are not trapped in the
core for small mixing of our interest, we assume $(1-f_{N})\simeq 1$ and
neglect the inverse processes in the Boltzmann equation. The spectrum of HNLs
per volume and time are given by eq. (16) after numerically solving the
Boltzmann equation (24). The 9-dimensional integrals in the collision term of
eq. (25) can be analytically reduced to 2-dimensional integrals, following
ref. [124].

The decay rates of HNLs to SM particles are given by refs. [125, 126, 127,
128, 129]. For HNLs with masses of $\mathcal{O}(100)\ {\rm MeV}$ and
$U_{\tau}\neq 0$, the dominant decay processes are
$N\rightarrow\overset{(-)}{\nu}_{\tau}\pi_{0}$,
$N\rightarrow\overset{(-)}{\nu}_{\tau}\nu_{\alpha}\bar{\nu}_{\alpha}\
(\alpha=e,\mu,\tau)$ and
$N\rightarrow\overset{(-)}{\nu}_{\tau}\ell^{+}\ell^{-}\ (\ell=e,\mu)$. For
HNLs with masses of $\mathcal{O}(100)\ {\rm MeV}$ and $U_{\mu}\neq 0$, the
dominant decay processes are $N\rightarrow\overset{(-)}{\nu}_{\mu}\pi_{0}$,
$N\rightarrow\overset{(-)}{\nu}_{\mu}\nu_{\alpha}\bar{\nu}_{\alpha}\
(\alpha=e,\mu,\tau)$,
$N\rightarrow\overset{(-)}{\nu}_{\mu}\ell^{+}\ell^{\prime-}\
(\ell,\ell^{\prime}=e,\mu)$ and $N\rightarrow\mu^{\pm}\pi^{\mp}$. The 2-body
decay is the most dominant but the 3-body decay channels can produce electron
neutrinos and antineutrinos, which can be detected most efficiently (and
affect the branching ratio).

For the 2-body decays of $N\rightarrow\nu_{\tau}\pi_{0}$, the distribution
function in eq. (23) is given by $f_{\nu}(E)=(1/2)\delta(E-\bar{E})$ with
$\bar{E}=(m_{N}^{2}-m_{\pi}^{2})/2m_{N}$. We obtain the secondary neutrino
fluence for $N\rightarrow\nu_{\tau}\pi_{0}$,

$\displaystyle\frac{dN_{\nu_{\tau}}}{dE_{\nu}}\biggl{|}_{\rm
2-body}=\frac{m_{N}^{2}}{m_{N}^{2}-m_{\pi}^{2}}{\rm
Br}_{\nu_{\tau}}\left[1-\exp\left(-\frac{E_{\nu}t_{\rm
data}}{\bar{E}\tau}\right)\right]\int^{\infty}_{E_{\rm
min}}dE_{N}\frac{1}{p_{N}}\frac{dN_{N}}{dE_{N}},$ (26)

where

$\displaystyle\bar{E}=\frac{m_{N}^{2}-m_{\pi}^{2}}{2m_{N}},\ \ \ \ E_{\rm
min}=m_{N}\frac{E_{\nu}^{2}+\bar{E}^{2}}{2E_{\nu}\bar{E}}.$ (27)

For 3-body decays, we have 2 types of the distribution function
$f(E,\cos\theta)$ in eq. (23) in the rest frame of HNLs 333These formulae are
slightly different from those presented in refs. [12, 41].,

$\displaystyle f_{I}(E)$
$\displaystyle=\frac{1}{2}16\frac{E^{2}}{m_{N}^{3}}\left(3-4\frac{E}{m_{N}}\right),$
(28) $\displaystyle f_{II}(E)$
$\displaystyle=\frac{1}{2}96\frac{E^{2}}{m_{N}^{3}}\left(1-2\frac{E}{m_{N}}\right).$
(29)

Eqs. (28) and (29) are normalized to be 1/2 in the integral from 0 to
$m_{N}/2$, which is the kinematically allowed range for $\nu$. $f_{I}$
corresponds to the spectrum for both neutrinos in the decays
$N\rightarrow\nu_{\beta}\nu_{\alpha}\bar{\nu}_{\alpha}$ while $f_{II}$ is the
one for antineutrinos in this process. For the charge-conjugated process
$N\rightarrow\bar{\nu}_{\beta}\bar{\nu}_{\alpha}\nu_{\alpha}$, $f_{I}$ is the
spectrum for antineutrinos while $f_{II}$ is the one for neutrinos. We also
use eqs. (28) and (29) for $N\rightarrow\nu_{\beta}\ell^{-}\ell^{\prime+}$,
neglecting the charged lepton mass for simplicity.

Process | $S|\mathcal{M}|^{2}/(16G_{F}^{2}V_{ud}^{2}U_{\mu N,\tau N}^{2})$
---|---
$\overset{(-)}{\nu}n\rightarrow nN$ | $(G_{V}^{n}+G_{A}^{n})^{2}(p_{1}\cdot p_{2})(p_{3}\cdot p_{4})+(G_{V}^{n}-G_{A}^{n})^{2}(p_{1}\cdot p_{3})(p_{2}\cdot p_{4})-(G_{V}^{n2}+G_{A}^{n2})m_{n}^{2}(p_{1}\cdot p_{4})$
$\overset{(-)}{\nu}p\rightarrow pN$ | $(G_{V}^{p}+G_{A}^{p})^{2}(p_{1}\cdot p_{2})(p_{3}\cdot p_{4})+(G_{V}^{p}-G_{A}^{p})^{2}(p_{1}\cdot p_{3})(p_{2}\cdot p_{4})-(G_{V}^{p2}+G_{A}^{p2})m_{p}^{2}(p_{1}\cdot p_{4})$
$\mu^{-}p\rightarrow nN$ | $(g_{V}+g_{A})^{2}(p_{1}\cdot p_{2})(p_{3}\cdot p_{4})+(g_{V}-g_{A})^{2}(p_{1}\cdot p_{3})(p_{2}\cdot p_{4})-(g_{V}^{2}+g_{A}^{2})m_{n}m_{p}(p_{1}\cdot p_{4})$
Table 1: The squared matrix elements $S|\mathcal{M}|^{2}$ for the dominant
production processes of HNLs mixing with $\nu_{\mu,\tau}$ in the SN core.
Subscript momentum number correspond to reactions labeled as $1+2\rightarrow
3+4$. The last process is valid only for HNLs mixing with $\nu_{\mu}$.
$g_{V}=1$, $g_{A}=1.27$, $G_{V}^{n}=1/2$, $G_{V}^{p}=1/2-2\sin^{2}\theta_{W}$
and $G_{A}^{n}=G_{A}^{p}=g_{A}/2$ [130]. $V_{ud}$ and $\sin\theta_{W}$ are the
up-down element of the Cabibbo-Kobayashi-Maskawa matrix and the Weinberg
angle, respectively.

###  4.3 ${\rm U(1)}_{L_{\mu}-L_{\tau}}$ gauge boson

We consider the extension of the SM to include a spontaneously broken ${\rm
U(1)}_{L_{\mu}-L_{\tau}}$ gauge symmetry. Then a massive ${\rm
U(1)}_{L_{\mu}-L_{\tau}}$ gauge boson, denoted by $Z^{\prime}$, is introduced
and the effective Lagrangian for ${\rm U(1)}_{L_{\mu}-L_{\tau}}$ gauge bosons
is described by

$\displaystyle\mathcal{L}_{L_{\mu}-L_{\tau}}$
$\displaystyle=\mathcal{L}_{SM}-\frac{1}{4}Z^{\prime}_{\mu\nu}Z^{\prime\mu\nu}-\frac{\varepsilon}{2}Z^{\prime}_{\mu\nu}F^{\mu\nu}+\frac{m_{Z^{\prime}}^{2}}{2}Z^{\prime}_{\mu}Z^{\prime\mu}$
$\displaystyle\ \ \ \
+g_{\mu-\tau}Z^{\prime}_{\alpha}(\bar{\mu}\gamma^{\alpha}\mu+\bar{\nu}_{\mu}\gamma^{\alpha}P_{L}\nu_{\mu}-\bar{\tau}\gamma^{\alpha}\tau-\bar{\nu}_{\tau}\gamma^{\alpha}P_{L}\nu_{\tau}),$
(30)

where $g_{\mu-\tau}$ is the gauge coupling, $P_{L}=\frac{1}{2}(1-\gamma_{5})$
is the left chirality projector, $\varepsilon$ is the kinetic mixing
parameter, and $F_{\alpha\beta}$ and $Z^{\prime}_{\alpha\beta}$ are the field
strength for the photon and $Z^{\prime}$, respectively. The decay rates for
$Z^{\prime}\rightarrow\bar{\nu}_{\alpha}\nu_{\alpha}\ (\alpha=\mu,\tau)$ and
$Z^{\prime}\rightarrow\ell^{+}\ell^{-}\ (\ell=\mu,\tau)$ in the rest frame are

$\displaystyle\Gamma_{Z^{\prime}\rightarrow\bar{\nu}_{\alpha}\nu_{\alpha}}=\frac{g_{\mu-\tau}^{2}}{24\pi}m_{Z^{\prime}},\
\ \ \
\Gamma_{Z^{\prime}\rightarrow\ell^{+}\ell^{-}}=\frac{g_{\mu-\tau}^{2}}{12\pi}m_{Z^{\prime}}\left(1+\frac{2m_{l}^{2}}{m_{Z^{\prime}}^{2}}\right)\sqrt{1-\frac{4m_{\ell}^{2}}{m_{Z^{\prime}}^{2}}}\,.$
(31)

In the canonically diagonalized basis, we obtain the coupling of $Z^{\prime}$
to the electromagnetic current

$\displaystyle\mathcal{L}\supset\varepsilon Z^{\prime}_{\alpha}J_{\rm
EM}^{\alpha},\ \ \ \ J_{\rm
EM}^{\alpha}=e\sum_{f}Q_{f}\bar{f}\gamma^{\alpha}f,$ (32)

where $e$ is the electric charge and $f$ is a SM fermion with charge $Q_{f}$.
We assume at the tree level this kinetic mixing vanishes. Nevertheless, the
muon and tau loops lead to this kinetic mixing. Such irreducible contributions
to $\varepsilon$ below energy scales of the muon mass are given by [131]

$\displaystyle\varepsilon=-\frac{eg_{\mu-\tau}}{2\pi^{2}}\int^{1}_{0}dx\
x(1-x)\log\left[\frac{m_{\tau}^{2}-x(1-x)q^{2}}{m_{\mu}^{2}-x(1-x)q^{2}}\right]\xrightarrow[m_{\mu}\gg
q]{}-\frac{eg_{\mu-\tau}}{12\pi^{2}}\log\frac{m_{\tau}^{2}}{m_{\mu}^{2}}\simeq-\frac{g_{\mu-\tau}}{70}.$
(33)

Note that model-dependent contributions could potentially arise if there exist
exotic particles charged under both ${\rm U}(1)_{\rm EM}$ and ${\rm
U(1)}_{L_{\mu}-L_{\tau}}$ [132].

The kinetic mixing can induce a decay channel of $Z^{\prime}\rightarrow
e^{+}e^{-}$ whose decay rate is given by

$\displaystyle\Gamma_{Z^{\prime}\rightarrow e^{+}e^{-}}=\frac{(\varepsilon
e)^{2}}{12\pi^{2}}m_{Z^{\prime}}\left(1+\frac{2m_{e}^{2}}{m_{Z^{\prime}}^{2}}\sqrt{1-\frac{4m_{e}^{2}}{m_{Z^{\prime}}^{2}}}\right)\,,$
(34)

and the branching ratio of $Z^{\prime}\rightarrow e^{+}e^{-}$ for $m_{e}\ll
m_{Z^{\prime}}\ll m_{\mu}$ is

$\displaystyle{\rm Br}_{Z^{\prime}\rightarrow
e^{+}e^{-}}=\frac{\Gamma_{Z^{\prime}\rightarrow
e^{+}e^{-}}}{\Gamma_{Z^{\prime}\rightarrow\bar{\nu}_{\mu}\nu_{\mu}}+\Gamma_{Z^{\prime}\rightarrow\bar{\nu}_{\tau}\nu_{\tau}}}\simeq\left(\frac{\varepsilon
e}{g_{\mu-\tau}}\right)^{2}\simeq 2\times 10^{-5},$ (35)

where we use the value of $\varepsilon$ in Eq. (33). Effects of the relatively
small kinetic mixing of $|\varepsilon|\lesssim g_{\mu-\tau}/70$ on supernovae
would be sub-dominant [9]. In the following we consider the kinetic mixing of
$|\varepsilon|\lesssim g_{\mu-\tau}/70$.

The main production processes for $Z^{\prime}$ in the SN core are the
neutrino-pair coalescence ($\bar{\nu}\nu\rightarrow Z^{\prime}$) and the semi-
Compton scattering ($\mu\gamma\rightarrow\mu Z^{\prime}$) [9]. Subsequently,
the decay processes of $Z^{\prime}\rightarrow\bar{\nu}\nu$ (and
$Z^{\prime}\rightarrow\mu^{+}\mu^{-}$ if $m_{Z}^{\prime}>2m_{\mu}$) occur. The
emitted $Z^{\prime}$ spectrum per volume and time from SNe is given by

$\displaystyle\frac{dn_{Z^{\prime}}}{dE_{Z^{\prime}}dt}$
$\displaystyle=\sum_{\nu=\nu_{\mu},\nu_{\tau}}\frac{g_{\mu-\tau}^{2}m_{Z^{\prime}}^{2}}{48\pi^{3}}\int_{E_{\nu}^{\rm
min}}^{E_{\nu}^{\rm max}}dE_{\nu}\
f_{\nu}(E_{\nu})f_{\bar{\nu}}(E_{Z^{\prime}}-E_{\nu})+\frac{E_{Z^{\prime}}p_{Z^{\prime}}}{2\pi^{2}}\mathcal{C}_{\mu\gamma\rightarrow\mu
Z^{\prime}}[f],$ (36)

where $f_{\nu}(E)$ is the Fermi-Dirac distribution for $\nu$, $E_{\nu}^{\rm
max}=\frac{1}{2}(E_{Z^{\prime}}+p_{Z^{\prime}})$ and $E_{\nu}^{\rm
min}=\frac{1}{2}(E_{Z^{\prime}}-p_{Z^{\prime}})$.
$\mathcal{C}_{\mu\gamma\rightarrow\mu Z^{\prime}}[f]\propto g_{\mu-\tau}^{2}$
is the collision term for $\mu\gamma\rightarrow\mu Z^{\prime}$, which is given
in the same manner of Eq. (25) except for the spin of particles. The first
term in Eq. (36) stems from neutrino-pair coalescence, while the second term
stems from the semi-Compton scattering. To calculate the second term, we use
an analogy of the Compton scattering ($\gamma\mu\rightarrow\gamma\mu$),
assuming $m_{Z^{\prime}}\simeq 0$. This approximation is valid since as
production of $Z^{\prime}$, the semi-Compton scattering only dominates at
small $m_{Z^{\prime}}$. The squared matrix element for
$\mu\gamma\rightarrow\mu Z^{\prime}$ with $m_{Z^{\prime}}=0$ is listed in
table 2. The 9-dimensional integrals in the collision term can be analytically
reduced to 3-dimensional integrals, following ref. [123].

The secondary neutrino fluence for $\nu_{\alpha}$ ($\alpha=\mu,\tau$) produced
by their decays, $Z^{\prime}\rightarrow\nu\bar{\nu}$, is given by, as in eq.
(26),

$\displaystyle\frac{dN_{\nu_{\alpha}}}{dE_{\nu}}={\rm
Br}_{\nu_{\alpha}}\left[1-\exp\left(-\frac{2E_{\nu}t_{\rm
data}}{m_{Z^{\prime}}\tau}\right)\right]\int^{\infty}_{E_{\rm
min}}dE_{Z^{\prime}}\frac{1}{p_{Z^{\prime}}}\frac{dN_{Z^{\prime}}}{dE_{Z^{\prime}}}\simeq{\rm
Br}_{\nu_{\alpha}}\int^{\infty}_{E_{\rm
min}}dE_{Z^{\prime}}\frac{1}{p_{Z^{\prime}}}\frac{dN_{Z^{\prime}}}{dE_{Z^{\prime}}},$
(37)

where $E_{\rm min}=E_{\nu}+m_{Z^{\prime}}^{2}/(4E_{\nu})$. The time delay in
eq. (21) is well suppressed by the Lorentz boost factor for
$Z^{\prime}\rightarrow\nu\bar{\nu}$ and the suppression factor by the time
delay in eq. (37) is negligible in the parameter region of our interest with
$t_{\rm data}\gtrsim 10\ {\rm s}$ for SN 1987A neutrino observations and
$t_{\rm data}\gtrsim 10^{3}\ {\rm s}$ for HK.

Process | $S|\mathcal{M}|^{2}/(8e^{2}g_{\mu-\tau}^{2})$
---|---
$\mu\gamma\rightarrow\mu Z^{\prime}$ | $\frac{p_{1}\cdot p_{4}}{p_{3}\cdot p_{4}}+\frac{p_{3}\cdot p_{4}}{p_{1}\cdot p_{4}}+2m_{\mu}^{2}\left(\frac{1}{p_{3}\cdot p_{4}}-\frac{1}{p_{1}\cdot p_{4}}\right)+m_{\mu}^{4}\left(\frac{1}{p_{3}\cdot p_{4}}-\frac{1}{p_{1}\cdot p_{4}}\right)^{2}$
Table 2: The squared matrix elements $S|\mathcal{M}|^{2}$ for
$\mu\gamma\rightarrow\mu Z^{\prime}$ assuming $m_{Z^{\prime}}\simeq 0$.
Notation as in Table 1.

###  4.4 ${\rm U(1)}_{B-L}$ gauge boson

The other $Z^{\prime}$ model of our interest is the ${\rm U(1)}_{B-L}$ gauge
extension of the SM. The Lagrangian for the ${\rm U(1)}_{B-L}$ model is

$\displaystyle\mathcal{L}_{B-L}$
$\displaystyle=\mathcal{L}_{SM}-\frac{1}{4}Z^{\prime}_{\mu\nu}Z^{\prime\mu\nu}-\frac{\varepsilon}{2}Z^{\prime}_{\mu\nu}F^{\mu\nu}+\frac{m_{Z^{\prime}}^{2}}{2}Z_{\mu}^{\prime}Z^{\prime\mu}+g_{B-L}Z^{\prime}_{\mu}\left(\frac{1}{3}\bar{q}\gamma^{\mu}q-\bar{\ell}\gamma^{\mu}\ell\right),$
(38)

where $g_{B-L}$ is the ${\rm U(1)}_{B-L}$ gauge coupling, and $q$ and $\ell$
denote all the quark and lepton fields, respectively; see refs. [16, 17] for
the details of the low-energy effective interactions in a medium. Right-handed
neutrinos have to be introduced to satisfy the anomaly cancellation condition.
In this analysis, we assume that they are heavy enough to be integrated out,
thus it gives no phenomenological consequence in the low-energy regime. The
decay rates for $Z^{\prime}\rightarrow\bar{\nu}_{\alpha}\nu_{\alpha}$
($\alpha=e,\mu,\tau$) and $Z^{\prime}\rightarrow\ell^{+}\ell^{-}$
($\ell=e,\mu,\tau$) are the same with eq. (31). Their decay to neutrinos
occurs in a flavor-universal way, so that the effects of neutrino oscillations
are averaged out.

The dominant production processes for $Z^{\prime}$ in the core are nucleon-
nucleon bremsstrahlung mediated by strong interactions ($pp\rightarrow
ppZ^{\prime}$, $nn\rightarrow nnZ^{\prime}$ and $np\rightarrow npZ^{\prime}$)
[16] and pion-induced reactions ($\pi^{-}p\rightarrow Z^{\prime}n$) [17].
However, the production of $Z^{\prime}$ via nucleon bremsstrahlung is only
formulated under the one-pion exchange approximation, in which nucleons
interact themselves via a one-pion mediator, leading to theoretical
uncertainties of an order of magnitude (see ref. [6] for the detailed
discussions in the axion case). In the following, we only consider pion-
induced reactions as the production of $Z^{\prime}$, following ref. [17]. The
secondary neutrino fluence for $\nu_{\alpha}$ ($\alpha=e,\mu,\tau$) by their
decays is described in the same way as in Eq. (37).

###  4.5 Majoron

We consider a (pseudo-)scalar boson $\phi$ coupled to neutrinos. Its
interactions can be described by

$\displaystyle\mathcal{L}_{\rm
int}=\frac{1}{2}g_{\alpha\beta}\bar{\nu}_{\beta}\nu_{\alpha}\phi,$ (39)

where $\nu$ is a 4-component field, and $\alpha,\beta=e,\mu,\tau$. In the
following, we assume Majorana neutrinos. For Majorana neutrinos, $\nu$
consists of left-handed active neutrinos and right-handed active neutrinos.

The dominant production processes in the SN core are
$\nu_{\alpha}\nu_{\beta}\rightarrow\phi$ and
$\bar{\nu}_{\alpha}\bar{\nu}_{\beta}\rightarrow\phi$ [26] 444As mentioned in
ref. [26], $\nu\nu\rightarrow\phi$ might stall because neutrinos might be
depleted in the core due to this process. However, the SM weak interaction
processes would supply neutrinos. In our analysis, we neglect any feedback to
the core by hypothesis particle productions as the typical analysis of SN
energy loss argument.. Subsequently, the decay processes
$\phi\rightarrow\nu_{\alpha}\nu_{\beta}$ and
$\phi\rightarrow\bar{\nu}_{\alpha}\bar{\nu}_{\beta}$ occur for massive $\phi$.
The decay rate for $\phi\rightarrow\nu_{\alpha}\nu_{\beta}$ in the rest frame
of $\phi$ in vacuum are

$\displaystyle\Gamma_{\phi\rightarrow\nu_{\alpha}\nu_{\beta}}$
$\displaystyle=\frac{g_{\alpha\beta}^{2}}{32\pi}m_{\phi}$ (40)

with $m_{\phi}$ the mass of $\phi$. The same expressions hold for the
processes $\phi\rightarrow\bar{\nu}_{\alpha}\bar{\nu}_{\beta}$.

The emitted spectrum per volume and time for $\phi$ from SNe is, following
refs. [26, 39],

$\displaystyle\frac{dn_{\phi}}{dE_{\phi}dt}=\sum\frac{g_{\alpha\beta}^{2}}{64\pi^{3}}\int
dE_{\nu}\
[m_{\phi}^{2}-2E_{\nu}(V_{\nu_{\alpha}}+V_{\nu_{\beta}})]f_{\nu_{\alpha}}(E_{\nu})f_{\nu_{\beta}}(E_{\phi}-E_{\nu}),$
(41)

where the summation is over all flavor neutrinos and antineutrinos,
$V_{\nu_{\alpha}}=-V_{\bar{\nu}_{\alpha}}$ is the effective potential for
neutrinos, which is $V_{\nu_{\alpha}}=V_{\rm CC}\delta_{\alpha e}+V_{\rm NC}$
with $V_{\rm CC}=\sqrt{2}G_{F}n_{e}$ and $V_{\rm
NC}=-\frac{1}{2}\sqrt{2}G_{F}n_{n}$. The contribution from protons is canceled
out by the contribution from electrons in an electrically neutral medium.
Typical values of $V_{\nu_{\alpha}}$ in the SN core are $V_{\nu_{e}}\simeq-10\
{\rm eV}$ and $V_{\nu_{\mu,\tau}}\simeq-20\ {\rm eV}$. The integration range
must satisfy the kinematically allowed relation with $V_{\nu_{\alpha}}\ll
E_{\nu}$ as follows

$\displaystyle 0\leq
1-\cos\theta=\frac{m_{\phi}^{2}-2E_{\phi}(V_{\nu_{\alpha}}+V_{\nu_{\beta}})}{2E_{\nu}(E_{\phi}-E_{\nu})}\leq
2,$ (42)

where $\theta$ is the relative angle between the momentums of a majoron and a
neutrino.

The secondary neutrino fluence for $\nu_{\alpha}$ produced by their decays is
given by, as in eq. (26),

$\displaystyle\frac{dN_{\nu_{\alpha}}}{dE_{\nu}}={\rm
Br}_{\nu_{\alpha}}\bar{N}_{\nu_{\alpha}}\int^{\infty}_{E_{\rm
min}}dE_{\phi}\frac{1}{p_{\phi}}\frac{dN_{\phi}}{dE_{\phi}},$ (43)

where $E_{\rm min}=E_{\nu}+m_{\phi}^{2}/4E_{\nu}$. The suppression factor by
the time delay as in eq. (37) is also negligible for $\phi\rightarrow\nu\nu$
in the parameter region of our interest.

##  5 Limits from SN 1987A and a future galactic supernova

###  5.1 Event rates and statistical analysis

As observations of SN 1987A neutrinos, we use the data from Kamiokande-II
(with a fiducial mass of 0.78 kton) [40] and IMB (6.8 kton) [133]. For future
galactic SN neutrino observations, we consider the detectors of Hyper-
Kamiokande (HK) (220 kton) [89]. SN neutrinos are dominantly captured by the
inverse beta decay (IBD), $\bar{\nu}_{e}+p\rightarrow e^{+}+n$, up to
$E_{\nu}\lesssim 70\ {\rm MeV}$ while for $E_{\nu}\gtrsim 70\ {\rm MeV}$, they
are captured by the charged current (CC) processes on oxygen,
$\bar{\nu}_{e}+\mathrm{O}\rightarrow e^{+}+\mathrm{X}$ and
$\nu_{e}+\mathrm{O}\rightarrow e^{-}+\mathrm{Y}$, where X and Y are nuclei at
the final state. We also account for the CC processes of $\nu_{\mu}$ and
$\bar{\nu}_{\mu}$ above the muon production threshold, $E_{\nu}\gtrsim
m_{\mu}=106\ {\rm MeV}$. The produced muons rapidly go to rest by ionization
and decay into $e^{\pm}$, whose spectrum is the known Michel $e^{\pm}$
spectrum ending at $m_{\mu}/2=53\ {\rm MeV}$. Below the muon Cherenkov
threshold of $E_{\mu^{\pm}}=160\ {\rm MeV}$, these events are called invisible
muons. This signal might be sub-dominant because its energy range overlaps
with that for the standard SN neutrinos. Above this threshold, visible
$\mu^{\pm}$ can contribute to the SN events. However, the Cherenkov threshold
behavior and the detection efficiency for visible $\mu^{\pm}$ and Michel
$e^{\pm}$ above $E_{\mu}=160\ {\rm MeV}$ in the IMB and Kamiokande-II are not
well known. We neglect the visible $\mu^{\pm}$ signal as in ref. [40], but the
result would not change significantly because the secondary neutrino fluence
is suppressed in the related energy region. We conservatively assume the data-
taking time for high-energy SN neutrino events in Kamiokande-II and IMB is
$t_{\rm data}=10\ {\rm s}$. For HK, we assume $t_{\rm data}=10^{3}\ {\rm s}$
and neglect background except for the standard SN neutrino events for
simplicity. We expect the other backgrounds such as atmospheric neutrino
events in the energy region of $1\textbf{--}1000\ {\rm MeV}$ are negligible
for $t_{\rm data}=10^{3}\ {\rm s}$ by rough rescaling a rate of 2/day in the
region of $20\textbf{--}2000\ {\rm MeV}$ for atmospheric neutrino events [86].
As a result, the value of $t_{\rm data}$ affects only the limit on HNLs in
this work.

The number of the events in an earth-based detector is approximately given by

$\displaystyle\frac{dN_{e}}{dE_{e}}(E_{e})$ $\displaystyle=\frac{1}{4\pi
d_{SN}^{2}}\sum_{x=p,\mathrm{O}}N_{x}$ $\displaystyle\ \ \ \
\times\left[\sum_{\nu=\nu_{e},\bar{\nu}_{e}}\frac{dN_{\nu}}{dE_{\nu}}(E_{\nu})\sigma_{\nu-x}(E_{e},E_{\nu})+\frac{dn_{e}}{dE_{e}}\sum_{\nu^{\prime}=\nu_{\mu},\bar{\nu}_{\mu}}\int
dE_{\nu}^{\prime}\frac{dN_{\nu^{\prime}}}{dE_{\nu^{\prime}}}(E_{\nu^{\prime}})\sigma_{\nu^{\prime}-x}(E_{e},E_{\nu^{\prime}})\right],$
(44)

where $N_{p}=6.69\times 10^{31}$ and $N_{\mathrm{O}}=3.35\times 10^{31}$ are
respectively the number of free proton and oxygen targets in 1 kton of the
water Cerenkov detector. $1/(4\pi d_{SN}^{2})$ is the geometrical dilution
factor with $d_{SN}=50\ {\rm kpc}$ for SN 1987A. We assume $d_{SN}=10\ {\rm
kpc}$ for a future galactic SN, which is the average distance [134].
$dN_{\nu}/dE_{\nu}$ is the neutrino spectrum from a supernova. $dn_{e}/dE_{e}$
is the Michel $e^{\pm}$ spectrum given by

$\displaystyle\frac{dn_{e}}{dE_{e}}=\frac{4}{m_{\mu}}\left(\frac{2E_{e}}{m_{\mu}}\right)^{2}\left(3-\frac{4E_{e}}{m_{\mu}}\right)$
(45)

with a cutoff at $E_{e}=m_{\mu}/2=53\ {\rm MeV}$. $\sigma_{\nu-x}$ is the
cross section between a neutrino and a nucleus. We take the IBD cross sections
for $\bar{\nu}_{e}$ from ref. [81] and for $\bar{\nu}_{\mu}$ from ref. [82],
and the CC cross sections for $\overset{(-)}{\nu}_{e}$ with
$\mathrm{{}^{16}O}$ from ref. [83] and for $\overset{(-)}{\nu}_{\mu}$ from
ref. [84] (see also FIG. S1 in ref. [40]). The number of events with detected
(reconstructed) energy with an energy resolution $\delta(E)$ and a detection
efficiency $\epsilon(E)$ is given by

$\displaystyle\frac{dN_{e}}{dE_{e}^{\text{det}}}(E_{e}^{\text{det}})=\int
dE_{e}\
\epsilon(E_{e})\frac{1}{\sqrt{2\pi}\delta(E_{e})}\text{exp}\bigg{[}-\frac{\big{(}E_{e}-E_{e}^{\text{det}}\big{)}^{2}}{2\delta^{2}(E_{e})}\bigg{]}\frac{dN_{e}}{dE_{e}}(E_{e}).$
(46)

The energy resolution $\delta(E)$ is taken as $\sqrt{(1.35\,{\rm MeV})E}$ for
IMB and $\sqrt{(0.75\,{\rm MeV})E}$ for Kamiokande-II [135]. For HK, we assume
that $\delta(E)$ is the same with the SK-IV observation,
$\delta(E)/E=0.0397(E/{\rm MeV})+0.349\sqrt{(E/{\rm MeV})}-0.0839$ [136]. The
detection efficiencies $\epsilon(E)$ are taken from ref. [135, 40] for IMB and
Kamiokande-II. For HK, we take $\epsilon(E)=1$.

To impose limits on light particle decays from SN 1987A and future galactic
SNe, we perform a maximum likelihood analysis, following refs. [135, 137, 40].
For water Cherenkov detectors, we define a likelihood,

$\displaystyle\mathcal{L}$ $\displaystyle\propto\exp\left[-\int^{E_{\rm
high}}_{E_{\rm low}}\frac{dN_{e}}{dE_{e}^{\rm det}}dE_{e}^{\rm
det}\right]\prod_{i}^{N_{\rm bin}}N_{e}^{N_{i}}$
$\displaystyle\propto\exp\left[-\int^{E_{\rm high}}_{E_{\rm
low}}\frac{dN_{e}}{dE_{e}^{\rm det}}dE_{e}^{\rm det}\right]\prod_{i}^{N_{\rm
obs}}\frac{dN_{e}}{dE_{e}^{\rm det}}(E_{i})$ (47)

where $N_{e}$ is the number of events with an energy bin, $E_{i}$ and $N_{i}$
are the actually observed energies and number of events with an energy bin.
$N_{\rm bin}$ and $N_{\rm obs}$ are the number of energy bins and the total
number of events. The normalization constant is not important because we will
consider only the likelihood ratio. For the SN 1987A, we consider the product
of the individual likelihoods of IMB and Kamiokande-II. According to the SN
1987A data, the relevant regions are $(E_{\rm low},E_{\rm high})=(7.5,50)\
{\rm MeV}$ for Kamiokande-II and $(19,75)\ {\rm MeV}$ for IMB [135, 40].
However, no other triggers except for muon-trigger events were observed in the
IMB [40]. We do not impose a upper cutoff of energy for the IMB. For HK, we
consider $(E_{\rm low},E_{\rm high})=(10,1000)\ {\rm MeV}$ and an energy bin
of 10 MeV. As the actual observational data for SN 1987A, we use FIG. S2 in
ref. [40]. As the expected data for a galactic SN, we use eq. (1) with
$(E_{\rm tot},E_{0},\alpha)=(1.98\times 10^{53}\ {\rm erg},\ 12.7\ {\rm MeV},\
2.39)$ for all $\nu$ and $\bar{\nu}$ for simplicity.

SN simulations with spherical symmetry have a tension with the SN 1987A data
[138]. In addition, the data is not informative about the pinching parameter
$\alpha$ [137]. So we fix $\alpha=2.39$ based on the SN model SFHo-18.8 and
conservatively marginalize the likelihood over the other SN model parameters,
$E_{0}$ and $E_{\rm tot}$,

$\displaystyle\tilde{\mathcal{L}}(g,m)=\max_{E_{0},E_{\rm
tot}}\mathcal{L}(g,m,E_{0},E_{\rm tot}),$ (48)

where $g$ and $m$ are the coupling and mass for hypothesis particles. We
confirm the marginalized likelihood $\tilde{L}$ is maximum at $g=0$, i.e., the
new signal is not preferred as in ref. [40]. Then we define a test statistics,

$\displaystyle\chi^{2}=2\left[\log\tilde{L}(0,m)-\log\tilde{L}(g,m_{\phi})\right].$
(49)

We compute threshold values of $(g,m)$ for 95% C.L. exclusion limits at
$\chi^{2}=2.7$.

###  5.2 Heavy neutral lepton

Figure 2 presents our limits (blue solid line) on HNLs mixing with $\nu_{\mu}$
(upper panel) and $\nu_{\tau}$ (bottom panel) from no observations for the
high-energy events of SN 1987A neutrinos in Kamiokande-II and IMB. We also
show the future sensitivity (blue dashed line) on HNLs from the observations
of galactic supernova neutrinos in HK with $d_{\rm SN}=10\ {\rm kpc}$. We
confirm the results in the NO and IO cases are the same within 20 $\%$ level
and show the conservative NO case in figure 2. In this figure we include other
relevant constraints from the other supernova arguments (green shaded region
enclosed by dotted line) [13], BBN & CMB assuming the standard cosmology
(light-gray shaded region) [120, 121, 122, 139, 13] and experiments (dark-gray
shaded region) [140]. We should note that the BBN constraints can be
significantly weakened in presence of a large neutrino asymmetry [141] or with
low reheating temperatures [142, 143].

Our constraints are one or two orders of magnitude weaker than the other SN
limits (in particular, the observations of $\gamma$-ray by HNL decays and the
excessive energy deposition by HNL decays in the SN envelope) [13]. We should
note that the SN constraints on HNLs highly depend on the SN core model
(temperature, chemical potentials, etc.) because the production rate of HNLs
in the core includes 9 phase-space integrals as in eq. (24). The difference
between the SN models would induce a large uncertainty and our reference SN
core model is different from ref. [13] 555We find that the constraints on HNLs
mixing with $\nu_{\mu}$ or $\nu_{\tau}$ from the SN energy loss with our
reference core model are a several factor weaker than those with the SN model
in ref. [13].. Compared with the constraint from SN 1987A energy loss with our
reference SN core model (not the result in ref. [13]), our new SN limit on the
mixing squared, $|U_{\mu N,\tau N}|^{2}$, is improved by a factor of 10–100,
depending on masses of HNLs. Observations of future galactic SN neutrinos in
HK with $d_{\rm SN}=10\ {\rm kpc}$ would have an order stronger sensitivities
on HNLs mixing with $U_{\tau N,\tau N}$ than the current SN limits in the mass
region of $300\ {\rm MeV}\lesssim m_{N}\lesssim 800\ {\rm MeV}$. These
sensitivities would improve our limits on $|U_{\mu N,\tau N}|^{2}$ by about 3
orders of magnitude in this mass region, roughly corresponding to the ratios
of the dilution factor due to the SN distance and the detector mass, $(10\
{\rm kpc}/50\ {\rm kpc})^{-2}\times(220\ {\rm kton}/6.8\ {\rm kton})\simeq
8\times 10^{2}$. For $m_{N}\lesssim 100\ {\rm MeV}$, the limit is suppressed
because the time delay of the secondary flux events from the standard ones
begins to exceed the data-taking time. For $m_{N}\gtrsim 100\ {\rm MeV}$, the
limit is also suppressed because the production of heavy $N$ is suppressed by
the Boltzmann factor.

Figure 2: SN1987A constraint on heavy neutral leptons mixing with $\nu_{\mu}$
(upper panel) and $\nu_{\tau}$ (bottom panel) from no observations of high
energy neutrino (high $\nu$) events (blue solid line) and future sensitivity
from observations of galactic high energy SN neutrinos in HK with $d_{\rm
SN}=10\ {\rm kpc}$ (blue dashed line). The results in the NO and IO cases are
the same within a few $10\%$ level. We also show the other SN limits (green
shaded region enclosed by dotted line) [13] and other relevant constraints
(see text for details).

###  5.3 ${\rm U(1)}_{L_{\mu}-L_{\tau}}$ gauge boson

Figure 3 presents our limit (blue solid lines) and future sensitivity in HK
with $d_{\rm SN}=10\ {\rm kpc}$ (blue dashed lines) on ${\rm
U(1)}_{L_{\mu}-L_{\tau}}$ gauge bosons, $Z^{\prime}$. Our constraint and
future sensitivity are applicable for the small kinetic mixing up to a natural
one, $|\varepsilon|\lesssim g_{\mu-\tau}/70$ in eq. (33). We confirm the
results in the NO and IO cases are the same within 20% level. In this figure
we include other relevant constraints from the other SN 1987A arguments (green
shaded region enclosed by dotted line) [9, 18], CMB [144], contribution of
$Z^{\prime}$ bosons to the effective number of neutrino species $\Delta N_{\rm
eff}$ [145] and stellar cooling, which is rescaled from [146, 147]. We also
show the preferred regions to explain the Hubble constant tension (yellow
band) [145]. Note that constraints by the stellar cooling stem from kinetic
mixing of $Z^{\prime}$ with photons and electrons and depend on the value of
the kinetic mixing, $\varepsilon$. In figure 3, the constraints by the
Borexino and stellar cooling with a natural value of
$\varepsilon=-g_{\mu-\tau}/70$ are shown. For $\varepsilon=0$, these
constraints are invalid.

Our lower constraint is the strongest in the mass region of $0.1\ {\rm
MeV}\lesssim m_{Z^{\prime}}\lesssim 400\ {\rm MeV}$. Compared with the
constraint from SN 1987A energy loss [9], our constraint on the coupling
$g_{\mu-\tau}$ is improved by a factor of 3–10. Observations of future
galactic SN neutrinos in HK with $d_{SN}=10\ {\rm kpc}$ would improve our
limit on the coupling by a factor of 26, corresponding to the improvements of
the SN distance squared times the detector mass. The results do not depend on
the data-taking time $t_{\rm data}$ because of the short time delay of the
secondary fluxes. For $m_{Z^{\prime}}\gtrsim 10\ {\rm MeV}$, the production
processes of neutrino-pair coalescence ($\nu\bar{\nu}\rightarrow Z^{\prime}$)
are dominant while the production process of semi-Compton scattering
($\mu\gamma\rightarrow\mu Z^{\prime}$) is dominant below $Z^{\prime}$ mass of
10 MeV. This is because the production rate for neutrino-pair coalescence in
eq. (36) is proportional to $m_{Z^{\prime}}^{2}$ while that for constant semi-
Compton scattering with light $Z^{\prime}$ does not depend on
$m_{Z^{\prime}}$.

Our upper limit is relaxed by a factor of a few compared with the SN limit in
ref. [9]. This is because our treatment of $Z^{\prime}\rightarrow\nu\bar{\nu}$
in the trapping case is different from that in ref. [9]. Note that in the
estimation of the production rate of $\mu\gamma\rightarrow\mu Z^{\prime}$ in
ref. [9], there is a missing factor of 2 and the authors assume the rest frame
of muons, leading to an enhancement of the production rate [10]. These
differences from our treatment also affect our improvement on ${\rm
U(1)}_{L_{\mu}-L_{\tau}}$ gauge bosons.

Figure 3: SN1987A constraint on ${\rm U(1)}_{L_{\mu}-L_{\tau}}$ gauge bosons
from no observations of high energy neutrino (high $\nu$) events (blue solid
line) and future sensitivity from observations of galactic high energy SN
neutrinos in HK with $d_{\rm SN}=10\ {\rm kpc}$ (blue dashed line). The
results are the same in the NO and IO cases within a few 10% level. Our
results are applicable for the kinetic mixing up to a natural one
$|\varepsilon|\lesssim g_{\mu-\tau}/70$. We show the constraint from SN 1987A
energy loss argument (green shaded region enclosed by dotted line) [9, 18] and
other constraints in the case of a natural kinetic mixing,
$\varepsilon=-g_{\mu-\tau}/70$ (see text for details). The preferred regions
to explain the $H_{0}$ tension (yellow band) [145] are also shown. Figure 4:
SN1987A constraints on ${\rm U(1)}_{B-L}$ gauge bosons from no observations of
high energy neutrino (high $\nu$) events (solid lines) and future sensitivity
from observations of galactic high energy SN neutrinos in HK with $d_{\rm
SN}=10\ {\rm kpc}$ (dashed lines). We also show the other SN limits (green
shaded region enclosed by dotted lines) [17], BBN [148, 15] and stellar
cooling [149, 147, 150].

###  5.4 ${\rm U(1)}_{B-L}$ gauge boson

Figure 4 presents our limits (blue solid lines) and future sensitivity in HK
with $d_{\rm SN}=10\ {\rm kpc}$ (blue dashed lines) on ${\rm U(1)}_{B-L}$
gauge bosons. In this case, $Z^{\prime}$ bosons decay to neutrinos flavor-
universally so that there is no effect of neutrino oscillations and our
results are independent of the neutrino mass ordering. In this figure we
include other relevant constraints from the other SN arguments (green shaded
region enclosed by dotted line) [17], BBN [148, 15] and stellar cooling [149,
147, 150]. We do not show the SN limit by the $Z^{\prime}$ production from
nucleon-nucleon bremsstrahlung because of the large uncertainty of the
production rate. See ref. [16] for interested readers. For $m_{Z^{\prime}}\geq
2m_{e}$, the observations of $\gamma$-ray by $Z^{\prime}$ decays produced in
the SN core and the excessive energy deposition by $Z^{\prime}$ decays in the
SN envelope have already imposed stringent constraints on $Z^{\prime}$ [17].
For $m_{Z^{\prime}}\geq T$, our constraint becomes weaker by the Boltzmann
suppression $e^{-m_{Z^{\prime}}/T}$. We only compute and show our constraints
and future sensitivity for $m_{Z^{\prime}}\lesssim 2m_{e}$ 666For
$m_{Z^{\prime}}\lesssim 2m_{e}$, the kinetic mixing coupling $\varepsilon$ is
constrained much more stringently by cosmology than SN 1987A observations
[151, 17]. We assume $\varepsilon=0$ in the ${\rm U(1)}_{B-L}$ model and
impose a supernova limit on $g_{B-L}$ in this work.. However, the future
sensitivity for $m_{Z^{\prime}}\geq T$ in HK may still reach the unconstrained
region. For the strong coupling in a trapping regime, our constraint is only
connected to the upper bound for the trapping regime in ref. [17].

Our lower limits are the most stringent in the region of $0.1\ {\rm
MeV}\lesssim m_{Z^{\prime}}\lesssim 1\ {\rm MeV}$. Compared with the SN energy
loss argument [17], our constraint on the coupling $g_{B-L}$ is improved by a
factor of 14. Observations of future galactic SN neutrinos in HK with
$d_{SN}=10\ {\rm kpc}$ would improve our limit on the coupling by a factor of
30, corresponding to the improvements of the SN distance squared and the
detector mass. The results are approximately independent of $t_{\rm data}$.
The previous SN limits are weaker for $m_{Z^{\prime}}\lesssim 2m_{e}\sim 1\
{\rm MeV}$. This is because $Z^{\prime}$ bosons with $m_{Z^{\prime}}\lesssim
2m_{e}$ cannot decay to charged particles and $\gamma$-ray telescopes cannot
observe a secondary flux from the produced charged particles.

Our improvement on ${\rm U(1)}_{B-L}$ gauge bosons is stronger than that on
${\rm U(1)}_{L_{\mu}-L_{\tau}}$ gauge bosons. This might be because the
produced ${\rm U(1)}_{B-L}$ gauge bosons would have a typical energy of $\sim
m_{\pi^{-}}+3T+\alpha_{n,p}\mu_{n,p}$ due to the production process of
$\pi^{-}p\rightarrow Z^{\prime}n$ while ${\rm U(1)}_{L_{\mu}-L_{\tau}}$ gauge
bosons would have a typical energy of $\sim 3T+\alpha_{\mu}\mu_{\mu}$.
$\alpha_{n,p,\mu}$ is a suppression factor of the chemical potential due to
the Pauli-blocking effects, but we could not estimate $\alpha_{n,p,\mu}$
quantitatively. Since the detection rate in the Water Cherenkov detectors is
roughly $\sigma_{\nu}\propto E_{\nu}^{2}$, neutrinos produced by the decays of
${\rm U(1)}_{B-L}$ gauge bosons might be detected more efficiently.

###  5.5 Majoron

Figure 5 presents our limits (solid lines) and future sensitivities in HK with
$d_{\rm SN}=10\ {\rm kpc}$ (dashed lines) on flavored majorons coupled to
neutrinos only with $g_{ee}$ (blue), $g_{ex}$ (yellow) or $g_{xx}$
($x=\mu,\tau$) (red). We confirm the constraints in the NO and IO cases are
the same within 20% level and the constraints for $x=\mu$ and $\tau$ are also
the same within 20% level. We include other relevant constraints from the SN
1987A energy loss argument (dotted lines) [27, 26] (we recalculate this
argument based on the SN model SFHo 18.8.), CMB [152, 144] and BBN [145, 153].
At $g\sim 10^{-4}$, majorons would enter a trapping regime by rescattering
$\nu\phi\rightarrow\nu\phi$ [27]. Our constraints are connected to the upper
limits for the trapping case in ref. [27]. For $m_{\phi}\lesssim 10\ {\rm
keV}$, the effective potential for neutrinos induced by the electron and
nucleus background in the core dominates the production rates of majorons,
changing the dispersion relation for neutrinos. For $m_{\phi}\lesssim 10\ {\rm
keV}$, the production rates and lower bounds do not depend on the majoron mass
$m_{\phi}$.

Our limits are the strongest in the region of $1\ {\rm eV}\lesssim
m_{\phi}\lesssim 500\ {\rm MeV}$. Compared with the constraint from SN 1987A
energy loss, our constraints on the couplings are improved by a factor of 13
for $g_{ee}$, 10 for $g_{ex}$, and 7 for $g_{xx}$. Observations of future
galactic SN neutrinos in HK with $d_{SN}=10\ {\rm kpc}$ would improve our
limits on the couplings by a factor of 30, corresponding to the improvements
of the SN distance squared and the detector mass. The results do not depend on
the data-taking time $t_{\rm data}$ because of the short time delay of the
secondary fluxes. In each flavor, the constraint on $g_{ee}$ is stronger than
those on $g_{ex}$ and $g_{xx}$ because a huge number of $\nu_{e}$ is emitted
in the SN core due to the large chemical potential for $\mu_{e}$ in the
duration of 2 s. The improvement for $g_{ee}$ is also stronger than those for
$g_{ex}$ and $g_{xx}$. This would be because the detection cross sections are
roughly proportional to $\sigma_{\nu}\propto E_{\nu}^{2}$ and the process
$\nu_{e}\nu_{e}\rightarrow\phi\rightarrow\nu_{e}\nu_{e}$ produces the most
energetic secondary neutrinos due to the large chemical potential of
$\nu_{e}$. Our constraint on $g_{ee}$ is very similar to the constraint on
majoron coupled to neutrinos flavor-universally [40] because in both cases,
the dominant production process of majoron in the core is
$\nu_{e}\nu_{e}\rightarrow\phi$.

Figure 5: SN1987A constraints on flavored majorons with couplings of $g_{ee}$
(blue), $g_{ex}$ (yellow) or $g_{xx}$ $(x=\mu,\tau)$ (red) from no
observations of high energy neutrino (high $\nu$) events (solid lines) and
future sensitivity from observations of galactic high energy SN neutrinos in
HK with $d_{\rm SN}=10\ {\rm kpc}$ (dashed lines). The results are the same in
the NO and IO cases and the cases of $x=\mu$ and $x=\tau$ within a few 10%
level. We also show constraints from the SN 1987A energy loss argument (dotted
lines) [27, 26], BBN [145, 153] and CMB [152, 144].

##  6 Summary

The secondary neutrino fluxes by decays of light hypothetical particles
produced in the SN core can modify the high energy tail of the standard SN
neutrino flux. The lack of high energy SN 1987A neutrino events imposes the
strong limit on flavor-universal neutrino non-standard interactions with light
bosons from SN 1987A neutrino observations [40] and can significantly improve
the sensitivities on the interactions from future galactic SN neutrino
observations [39].

In this work, we extend refs. [39, 40] to more diverse particle physics
models, studying neutrino oscillation effects in the SN envelope on the
secondary neutrino fluxes for the first time. In particular, we obtain the
strong constraints and future sensitivity for HNLs, ${\rm
U(1)}_{L_{\mu}-L_{\tau}}$ and ${\rm U(1)}_{B-L}$ gauge bosons, and flavored
majorons. In refs. [12, 41], the authors considered the secondary fluxes by
decays of HNLs and the future sensitivity from observations of galactic
supernovae. We improve many calculations such as the production rate of HNLs
in the SN core, neutrino oscillations and the event rates at the detector. We
also obtain the current limit on HNLs from SN 1987A observations using the
above argument for the first time.

For HNLs, we find limits as in figure 2 from the absence of high energy SN
1987A neutrino events but these limits are one or two orders of magnitude
weaker than the current other SN limits [13]. Future observations of galactic
SN neutrino events in HK with $d_{\rm SN}=10\ {\rm kpc}$ will have an order
stronger sensitivities on mixing with $\nu_{\tau}$ than the current SN
constraints in the mass region of $300\ {\rm MeV}\lesssim m_{N}\lesssim 800\
{\rm MeV}$. Note that the SN constraints on HNLs highly depend on the SN model
because of the high dependence of the production of HNLs on the energy of
incoming and outgoing particles.

For ${\rm U(1)}_{L_{\mu}-L_{\tau}}$ gauge bosons, the absence of high energy
SN 1987A neutrinos imposes the strongest limit as in figure 3. We exclude the
parameter region with masses of $\lesssim 400\ {\rm MeV}$ and coupling down to
$g_{\mu-\tau}\sim 10^{-10}$. For $m_{Z^{\prime}}\lesssim 1\ {\rm MeV}$, our
limit is an order of magnitude stronger than the limit from observations of
$\Delta N_{\rm eff}$ [80]. Future observations of galactic SN neutrinos in HK
will improve the limit on coupling down to $g_{\mu-\tau}\sim 10^{-11}$ in this
mass range.

For ${\rm U(1)}_{B-L}$ gauge bosons, the absence of high energy SN 1987A
neutrinos imposes the strongest limit as in figure 4. We newly exclude the
parameter region with masses of $\lesssim 1\ {\rm MeV}$ and coupling down to
$g_{B-L}=4\times 10^{-12}$. Future observations of galactic SN neutrinos in HK
will improve the limit on coupling down to $g_{B-L}=10^{-13}$ in this mass
range.

For flavored majorons, the absence of high energy SN 1987A neutrinos imposes
the strongest limit as in figure 5. We exclude the parameter region with
$g\sim 10^{-7}$ for $m_{\phi}\lesssim 10\ {\rm keV}$ and $m_{\phi}g\sim
10^{-9}\ {\rm MeV}$ for $10\ {\rm keV}\lesssim m_{\phi}\lesssim 500\ {\rm
MeV}$. Future observations of galactic SN neutrinos in HK will improve the
limit down to $g\sim 10^{-8}$ for $m_{\phi}\lesssim 10\ {\rm keV}$ and
$m_{\phi}g\sim 10^{-11}\ {\rm MeV}$ for $10\ {\rm keV}\lesssim
m_{\phi}\lesssim 800\ {\rm MeV}$.

## Acknowledgments

This work was supported by IBS under the project code, IBS-R018-D1. MM also
acknowledges support from the grant NRF-2022R1A2C1009686.

## References

  * [1] G.G. Raffelt, _Stars as laboratories for fundamental physics: The astrophysics of neutrinos, axions, and other weakly interacting particles_ , University of Chicago Press (1996).
  * [2] G.G. Raffelt, _Neutrinos and the stars_ , _Proc. Int. Sch. Phys. Fermi_ 182 (2012) 61 [1201.1637].
  * [3] G. Raffelt and D. Seckel, _Bounds on exotic-particle interactions from sn1987a_ , _Phys. Rev. Lett._ 60 (1988) 1793.
  * [4] W. Keil, H.-T. Janka, D.N. Schramm, G. Sigl, M.S. Turner and J.R. Ellis, _A Fresh look at axions and SN-1987A_ , _Phys. Rev. D_ 56 (1997) 2419 [astro-ph/9612222].
  * [5] J.H. Chang, R. Essig and S.D. McDermott, _Supernova 1987A Constraints on Sub-GeV Dark Sectors, Millicharged Particles, the QCD Axion, and an Axion-like Particle_ , _JHEP_ 09 (2018) 051 [1803.00993].
  * [6] P. Carenza, T. Fischer, M. Giannotti, G. Guo, G. Martínez-Pinedo and A. Mirizzi, _Improved axion emissivity from a supernova via nucleon-nucleon bremsstrahlung_ , _JCAP_ 10 (2019) 016 [1906.11844].
  * [7] P. Carenza, B. Fore, M. Giannotti, A. Mirizzi and S. Reddy, _Enhanced Supernova Axion Emission and its Implications_ , _Phys. Rev. Lett._ 126 (2021) 071102 [2010.02943].
  * [8] R. Bollig, W. DeRocco, P.W. Graham and H.-T. Janka, _Muons in Supernovae: Implications for the Axion-Muon Coupling_ , _Phys. Rev. Lett._ 125 (2020) 051104 [2005.07141].
  * [9] D. Croon, G. Elor, R.K. Leane and S.D. McDermott, _Supernova Muons: New Constraints on $Z$’ Bosons, Axions and ALPs_, _JHEP_ 01 (2021) 107 [2006.13942].
  * [10] A. Caputo, G. Raffelt and E. Vitagliano, _Muonic boson limits: Supernova redux_ , _Phys. Rev. D_ 105 (2022) 035022 [2109.03244].
  * [11] K. Choi, H.J. Kim, H. Seong and C.S. Shin, _Axion emission from supernova with axion-pion-nucleon contact interaction_ , _JHEP_ 02 (2022) 143 [2110.01972].
  * [12] L. Mastrototaro, A. Mirizzi, P.D. Serpico and A. Esmaili, _Heavy sterile neutrino emission in core-collapse supernovae: Constraints and signatures_ , _JCAP_ 01 (2020) 010 [1910.10249].
  * [13] P. Carenza, G. Lucente, L. Mastrototaro, A. Mirizzi and P.D. Serpico, _Comprehensive constraints on heavy sterile neutrinos from core-collapse supernovae_ , 2311.00033.
  * [14] J.H. Chang, R. Essig and S.D. McDermott, _Revisiting Supernova 1987A Constraints on Dark Photons_ , _JHEP_ 01 (2017) 107 [1611.03864].
  * [15] S. Knapen, T. Lin and K.M. Zurek, _Light Dark Matter: Models and Constraints_ , _Phys. Rev. D_ 96 (2017) 115021 [1709.07882].
  * [16] C.S. Shin and S. Yun, _Dark gauge boson production from neutron stars via nucleon-nucleon bremsstrahlung_ , _JHEP_ 02 (2022) 133 [2110.03362].
  * [17] C.S. Shin and S. Yun, _Dark gauge boson emission from supernova pions_ , _Phys. Rev. D_ 108 (2023) 055014 [2211.15677].
  * [18] D.G. Cerdeño, M. Cermeño and Y. Farzan, _Constraints from the duration of supernova neutrino burst on on-shell light gauge boson production by neutrinos_ , _Phys. Rev. D_ 107 (2023) 123012 [2301.00661].
  * [19] K. Choi, C.W. Kim, J. Kim and W.P. Lam, _Constraints on the Majoron Interactions From the Supernova SN1987A_ , _Phys. Rev. D_ 37 (1988) 3225.
  * [20] Z. Berezhiani and A. Smirnov, _Matter-induced neutrino decay and supernova 1987a_ , _Physics Letters B_ 220 (1989) 279.
  * [21] K. Choi and A. Santamaria, _Majorons and Supernova Cooling_ , _Phys. Rev. D_ 42 (1990) 293.
  * [22] S. Chang and K. Choi, _Constraints from nucleosynthesis and SN1987A on majoron emitting double beta decay_ , _Phys. Rev. D_ 49 (1994) 12 [hep-ph/9303243].
  * [23] M. Kachelriess, R. Tomas and J.W.F. Valle, _Supernova bounds on Majoron emitting decays of light neutrinos_ , _Phys. Rev. D_ 62 (2000) 023004 [hep-ph/0001039].
  * [24] R. Tomas, H. Pas and J.W.F. Valle, _Generalized bounds on Majoron - neutrino couplings_ , _Phys. Rev. D_ 64 (2001) 095005 [hep-ph/0103017].
  * [25] S. Hannestad, P. Keranen and F. Sannino, _A Supernova constraint on bulk Majorons_ , _Phys. Rev. D_ 66 (2002) 045002 [hep-ph/0204231].
  * [26] Y. Farzan, _Bounds on the coupling of the Majoron to light neutrinos from supernova cooling_ , _Phys. Rev. D_ 67 (2003) 073015 [hep-ph/0211375].
  * [27] L. Heurtier and Y. Zhang, _Supernova Constraints on Massive (Pseudo)Scalar Coupling to Neutrinos_ , _JCAP_ 02 (2017) 042 [1609.05882].
  * [28] T. Brune and H. Päs, _Massive Majorons and constraints on the Majoron-neutrino coupling_ , _Phys. Rev. D_ 99 (2019) 096005 [1808.08158].
  * [29] J.A. Grifols, E. Masso and R. Toldra, _Gamma-rays from SN1987A due to pseudoscalar conversion_ , _Phys. Rev. Lett._ 77 (1996) 2372 [astro-ph/9606028].
  * [30] J.W. Brockway, E.D. Carlson and G.G. Raffelt, _SN1987A gamma-ray limits on the conversion of pseudoscalars_ , _Phys. Lett. B_ 383 (1996) 439 [astro-ph/9605197].
  * [31] A. Payez, C. Evoli, T. Fischer, M. Giannotti, A. Mirizzi and A. Ringwald, _Revisiting the SN1987A gamma-ray limit on ultralight axion-like particles_ , _JCAP_ 02 (2015) 006 [1410.3747].
  * [32] M. Meyer, M. Giannotti, A. Mirizzi, J. Conrad and M.A. Sánchez-Conde, _Fermi Large Area Telescope as a Galactic Supernovae Axionscope_ , _Phys. Rev. Lett._ 118 (2017) 011103 [1609.02350].
  * [33] J. Jaeckel, P.C. Malta and J. Redondo, _Decay photons from the axionlike particles burst of type II supernovae_ , _Phys. Rev. D_ 98 (2018) 055032 [1702.02964].
  * [34] F. Calore, P. Carenza, M. Giannotti, J. Jaeckel and A. Mirizzi, _Bounds on axionlike particles from the diffuse supernova flux_ , _Phys. Rev. D_ 102 (2020) 123005 [2008.11741].
  * [35] R.Z. Ferreira, M.C.D. Marsh and E. Müller, _Strong supernovae bounds on ALPs from quantum loops_ , _JCAP_ 11 (2022) 057 [2205.07896].
  * [36] M. Diamond, D.F.G. Fiorillo, G. Marques-Tavares and E. Vitagliano, _Axion-sourced fireballs from supernovae_ , _Phys. Rev. D_ 107 (2023) 103029 [2303.11395].
  * [37] W. DeRocco, P.W. Graham, D. Kasen, G. Marques-Tavares and S. Rajendran, _Observable signatures of dark photons from supernovae_ , _JHEP_ 02 (2019) 171 [1901.08596].
  * [38] F. Calore, P. Carenza, M. Giannotti, J. Jaeckel, G. Lucente, L. Mastrototaro et al., _511 keV line constraints on feebly interacting particles from supernovae_ , _Phys. Rev. D_ 105 (2022) 063026 [2112.08382].
  * [39] K. Akita, S.H. Im and M. Masud, _Probing non-standard neutrino interactions with a light boson from next galactic and diffuse supernova neutrinos_ , _JHEP_ 12 (2022) 050 [2206.06852].
  * [40] D.F.G. Fiorillo, G.G. Raffelt and E. Vitagliano, _Strong Supernova 1987A Constraints on Bosons Decaying to Neutrinos_ , _Phys. Rev. Lett._ 131 (2023) 021001 [2209.11773].
  * [41] V. Syvolap, _Testing heavy neutral leptons produced in the supernovae explosions with future neutrino detectors_ , 2301.07052.
  * [42] S.W. Falk and D.N. Schramm, _Limits From Supernovae on Neutrino Radiative Lifetimes_ , _Phys. Lett. B_ 79 (1978) 511.
  * [43] A. Sung, H. Tu and M.-R. Wu, _New constraint from supernova explosions on light particles beyond the Standard Model_ , _Phys. Rev. D_ 99 (2019) 121305 [1903.07923].
  * [44] A. Caputo, H.-T. Janka, G. Raffelt and E. Vitagliano, _Low-Energy Supernovae Severely Constrain Radiative Particle Decays_ , _Phys. Rev. Lett._ 128 (2022) 221103 [2201.09890].
  * [45] G. Chauhan, S. Horiuchi, P. Huber and I.M. Shoemaker, _Low-Energy Supernovae Bounds on Sterile Neutrinos_ , 2309.05860.
  * [46] P.-W. Chang, I. Esteban, J.F. Beacom, T.A. Thompson and C.M. Hirata, _Toward Powerful Probes of Neutrino Self-Interactions in Supernovae_ , _Phys. Rev. Lett._ 131 (2023) 071002 [2206.12426].
  * [47] D.F.G. Fiorillo, G. Raffelt and E. Vitagliano, _Supernova Emission of Secretly Interacting Neutrino Fluid: Theoretical Foundations_ , 2307.15122.
  * [48] D.F.G. Fiorillo, G. Raffelt and E. Vitagliano, _Large Neutrino Secret Interactions, Small Impact on Supernovae_ , 2307.15115.
  * [49] A.M. Abdullahi et al., _The present and future status of heavy neutral leptons_ , _J. Phys. G_ 50 (2023) 020501 [2203.08039].
  * [50] R. Foot, _New Physics From Electric Charge Quantization?_ , _Mod. Phys. Lett. A_ 6 (1991) 527.
  * [51] X.G. He, G.C. Joshi, H. Lew and R.R. Volkas, _NEW Z-prime PHENOMENOLOGY_ , _Phys. Rev. D_ 43 (1991) 22.
  * [52] X.-G. He, G.C. Joshi, H. Lew and R.R. Volkas, _Simplest Z-prime model_ , _Phys. Rev. D_ 44 (1991) 2118.
  * [53] A. Davidson, _$B-L$ as the fourth color within an $\mathrm{SU}(2)_{L}\times\mathrm{U}(1)_{R}\times\mathrm{U}(1)$ model_, _Phys. Rev. D_ 20 (1979) 776.
  * [54] R.N. Mohapatra and R.E. Marshak, _Local B-L Symmetry of Electroweak Interactions, Majorana Neutrinos and Neutron Oscillations_ , _Phys. Rev. Lett._ 44 (1980) 1316.
  * [55] C. Wetterich, _Neutrino Masses and the Scale of B-L Violation_ , _Nucl. Phys. B_ 187 (1981) 343.
  * [56] W. Buchmuller, C. Greub and P. Minkowski, _Neutrino masses, neutral vector bosons and the scale of B-L breaking_ , _Phys. Lett. B_ 267 (1991) 395.
  * [57] Y. Chikashige, R.N. Mohapatra and R.D. Peccei, _Are There Real Goldstone Bosons Associated with Broken Lepton Number?_ , _Phys. Lett. B_ 98 (1981) 265.
  * [58] G.B. Gelmini and M. Roncadelli, _Left-Handed Neutrino Mass Scale and Spontaneously Broken Lepton Number_ , _Phys. Lett. B_ 99 (1981) 411.
  * [59] J. Schechter and J.W.F. Valle, _Neutrino Decay and Spontaneous Violation of Lepton Number_ , _Phys. Rev. D_ 25 (1982) 774.
  * [60] J. Asaadi, E. Church, R. Guenette, B.J.P. Jones and A.M. Szelc, _New light Higgs boson and short-baseline neutrino anomalies_ , _Phys. Rev. D_ 97 (2018) 075021 [1712.08019].
  * [61] B. Chauhan and S. Mohanty, _Signature of light sterile neutrinos at IceCube_ , _Phys. Rev. D_ 98 (2018) 083021 [1808.04774].
  * [62] A.Y. Smirnov and V.B. Valera, _Resonance refraction and neutrino oscillations_ , _JHEP_ 09 (2021) 177 [2106.13829].
  * [63] M. Dentler, I. Esteban, J. Kopp and P. Machado, _Decaying Sterile Neutrinos and the Short Baseline Oscillation Anomalies_ , _Phys. Rev. D_ 101 (2020) 115013 [1911.01427].
  * [64] A. de Gouvêa, O.L.G. Peres, S. Prakash and G.V. Stenico, _On The Decaying-Sterile Neutrino Solution to the Electron (Anti)Neutrino Appearance Anomalies_ , _JHEP_ 07 (2020) 141 [1911.01447].
  * [65] Y.S. Jeong, S. Palomares-Ruiz, M.H. Reno and I. Sarcevic, _Probing secret interactions of eV-scale sterile neutrinos with the diffuse supernova neutrino background_ , _JCAP_ 06 (2018) 019 [1803.04541].
  * [66] W. Abdallah, R. Gandhi and S. Roy, _Requirements on common solutions to the LSND and MiniBooNE excesses: a post-MicroBooNE study_ , _JHEP_ 06 (2022) 160 [2202.09373].
  * [67] Muon g-2 collaboration, _Final Report of the Muon E821 Anomalous Magnetic Moment Measurement at BNL_ , _Phys. Rev. D_ 73 (2006) 072003 [hep-ex/0602035].
  * [68] T. Araki, F. Kaneko, T. Ota, J. Sato and T. Shimomura, _MeV scale leptonic force for cosmic neutrino spectrum and muon anomalous magnetic moment_ , _Phys. Rev. D_ 93 (2016) 013014 [1508.07471].
  * [69] S. Borsanyi et al., _Leading hadronic contribution to the muon magnetic moment from lattice QCD_ , _Nature_ 593 (2021) 51 [2002.12347].
  * [70] Muon g-2 collaboration, _Measurement of the Positive Muon Anomalous Magnetic Moment to 0.46 ppm_ , _Phys. Rev. Lett._ 126 (2021) 141801 [2104.03281].
  * [71] J.A. Carpio, K. Murase, I.M. Shoemaker and Z. Tabrizi, _High-energy cosmic neutrinos as a probe of the vector mediator scenario in light of the muon g-2 anomaly and Hubble tension_ , _Phys. Rev. D_ 107 (2023) 103057 [2104.15136].
  * [72] L.G. van den Aarssen, T. Bringmann and C. Pfrommer, _Is dark matter with long-range interactions a solution to all small-scale problems of \Lambda CDM cosmology?_ , _Phys. Rev. Lett._ 109 (2012) 231301 [1205.5809].
  * [73] F.-Y. Cyr-Racine and K. Sigurdson, _Limits on Neutrino-Neutrino Scattering in the Early Universe_ , _Phys. Rev. D_ 90 (2014) 123533 [1306.1536].
  * [74] J.F. Cherry, A. Friedland and I.M. Shoemaker, _Neutrino Portal Dark Matter: From Dwarf Galaxies to IceCube_ , 1411.1071.
  * [75] X. Chu, B. Dasgupta and J. Kopp, _Sterile neutrinos with secret interactions—lasting friendship with cosmology_ , _JCAP_ 10 (2015) 011 [1505.02795].
  * [76] L. Lancaster, F.-Y. Cyr-Racine, L. Knox and Z. Pan, _A tale of two modes: Neutrino free-streaming in the early universe_ , _JCAP_ 07 (2017) 033 [1704.06657].
  * [77] X. Chu, B. Dasgupta, M. Dentler, J. Kopp and N. Saviano, _Sterile neutrinos with secret interactions—cosmological discord?_ , _JCAP_ 11 (2018) 049 [1806.10629].
  * [78] C.D. Kreisch, F.-Y. Cyr-Racine and O. Doré, _Neutrino puzzle: Anomalies, interactions, and cosmological tensions_ , _Phys. Rev. D_ 101 (2020) 123505 [1902.00534].
  * [79] E. Grohs, G.M. Fuller and M. Sen, _Consequences of neutrino self interactions for weak decoupling and big bang nucleosynthesis_ , _JCAP_ 07 (2020) 001 [2002.08557].
  * [80] M. Escudero, D. Hooper, G. Krnjaic and M. Pierre, _Cosmology with A Very Light L μ $-$ Lτ Gauge Boson_, _JHEP_ 03 (2019) 071 [1901.02010].
  * [81] A. Strumia and F. Vissani, _Precise quasielastic neutrino/nucleon cross-section_ , _Phys. Lett. B_ 564 (2003) 42 [astro-ph/0302055].
  * [82] J.A. Formaggio and G.P. Zeller, _From eV to EeV: Neutrino Cross Sections Across Energy Scales_ , _Rev. Mod. Phys._ 84 (2012) 1307 [1305.7513].
  * [83] E. Kolbe, K. Langanke and P. Vogel, _Estimates of weak and electromagnetic nuclear decay signatures for neutrino reactions in Super-Kamiokande_ , _Phys. Rev. D_ 66 (2002) 013007.
  * [84] J. Marteau, J. Delorme and M. Ericson, _Neutrino oxygen interactions: Role of nuclear physics in the atmospheric neutrino anomaly_ , in _34th Rencontres de Moriond: Electroweak Interactions and Unified Theories_ , pp. 121–126, 1999 [hep-ph/9906449].
  * [85] Kamiokande-II collaboration, _Observation of a Neutrino Burst from the Supernova SN 1987a_ , _Phys. Rev. Lett._ 58 (1987) 1490.
  * [86] R.M. Bionta et al., _Observation of a Neutrino Burst in Coincidence with Supernova SN 1987a in the Large Magellanic Cloud_ , _Phys. Rev. Lett._ 58 (1987) 1494.
  * [87] E.N. Alekseev, L.N. Alekseeva, V.I. Volchenko and I.V. Krivosheina, _Possible Detection of a Neutrino Signal on 23 February 1987 at the Baksan Underground Scintillation Telescope of the Institute of Nuclear Research_ , _JETP Lett._ 45 (1987) 589.
  * [88] K. Rozwadowska, F. Vissani and E. Cappellaro, _On the rate of core collapse supernovae in the milky way_ , _New Astron._ 83 (2021) 101498 [2009.03438].
  * [89] Hyper-Kamiokande collaboration, _Hyper-Kamiokande Design Report_ , 1805.04163.
  * [90] _Garching core-collapse supernova research archive_ , _https://wwwmpa.mpa-garching.mpg.de/ccsnarchive/_ .
  * [91] T. Fischer, P. Carenza, B. Fore, M. Giannotti, A. Mirizzi and S. Reddy, _Observable signatures of enhanced axion emission from protoneutron stars_ , _Phys. Rev. D_ 104 (2021) 103012 [2108.13726].
  * [92] B. Fore and S. Reddy, _Pions in hot dense matter and their astrophysical implications_ , _Phys. Rev. C_ 101 (2020) 035809 [1911.02632].
  * [93] D. Page, M.V. Beznogov, I. Garibay, J.M. Lattimer, M. Prakash and H.-T. Janka, _NS 1987A in SN 1987A_ , _Astrophys. J._ 898 (2020) 125 [2004.06078].
  * [94] M. Hempel, _Nucleon self-energies for supernova equations of state_ , _Phys. Rev. C_ 91 (2015) 055807 [1410.6337].
  * [95] R. Bollig, H.T. Janka, A. Lohs, G. Martinez-Pinedo, C.J. Horowitz and T. Melson, _Muon Creation in Supernova Matter Facilitates Neutrino-driven Explosions_ , _Phys. Rev. Lett._ 119 (2017) 242702 [1706.04630].
  * [96] M.T. Keil, G.G. Raffelt and H.-T. Janka, _Monte Carlo study of supernova neutrino spectra formation_ , _Astrophys. J._ 590 (2003) 971 [astro-ph/0208035].
  * [97] I. Tamborra, B. Muller, L. Hudepohl, H.-T. Janka and G. Raffelt, _High-resolution supernova neutrino spectra represented by a simple fit_ , _Phys. Rev. D_ 86 (2012) 125031 [1211.3920].
  * [98] A.S. Dighe and A.Y. Smirnov, _Identifying the neutrino mass spectrum from the neutrino burst from a supernova_ , _Phys. Rev. D_ 62 (2000) 033007 [hep-ph/9907423].
  * [99] H. Duan, G.M. Fuller and Y.-Z. Qian, _Collective Neutrino Oscillations_ , _Ann. Rev. Nucl. Part. Sci._ 60 (2010) 569 [1001.2799].
  * [100] A. Mirizzi, I. Tamborra, H.-T. Janka, N. Saviano, K. Scholberg, R. Bollig et al., _Supernova Neutrinos: Production, Oscillations and Detection_ , _Riv. Nuovo Cim._ 39 (2016) 1 [1508.00785].
  * [101] S. Chakraborty, R. Hansen, I. Izaguirre and G. Raffelt, _Collective neutrino flavor conversion: Recent developments_ , _Nucl. Phys. B_ 908 (2016) 366 [1602.02766].
  * [102] I. Tamborra and S. Shalgar, _New Developments in Flavor Evolution of a Dense Neutrino Gas_ , _Ann. Rev. Nucl. Part. Sci._ 71 (2021) 165 [2011.01948].
  * [103] T.-K. Kuo and J.T. Pantaleone, _Neutrino Oscillations in Matter_ , _Rev. Mod. Phys._ 61 (1989) 937.
  * [104] M. Blennow and A.Y. Smirnov, _Neutrino propagation in matter_ , _Adv. High Energy Phys._ 2013 (2013) 972485 [1306.2903].
  * [105] S.P. Mikheyev and A.Y. Smirnov, _Resonance Amplification of Oscillations in Matter and Spectroscopy of Solar Neutrinos_ , _Sov. J. Nucl. Phys._ 42 (1985) 913.
  * [106] S.E. Woosley and T.A. Weaver, _The Evolution and explosion of massive stars. 2. Explosive hydrodynamics and nucleosynthesis_ , _Astrophys. J. Suppl._ 101 (1995) 181.
  * [107] K. Nakazato, K. Sumiyoshi, H. Suzuki, T. Totani, H. Umeda and S. Yamada, _Supernova Neutrino Light Curves and Spectra for Various Progenitor Stars: From Core Collapse to Proto-neutron Star Cooling_ , _Astrophys. J. Suppl._ 205 (2013) 2 [1210.6841].
  * [108] I. Esteban, M.C. Gonzalez-Garcia, M. Maltoni, T. Schwetz and A. Zhou, _The fate of hints: updated global analysis of three-flavor neutrino oscillations_ , _JHEP_ 09 (2020) 178 [2007.14792].
  * [109] P.F. de Salas, D.V. Forero, S. Gariazzo, P. Martínez-Miravé, O. Mena, C.A. Ternes et al., _2020 global reassessment of the neutrino oscillation picture_ , _JHEP_ 02 (2021) 071 [2006.11237].
  * [110] G. Sigl and G. Raffelt, _General kinetic description of relativistic mixed neutrinos_ , _Nucl. Phys. B_ 406 (1993) 423.
  * [111] K. Akita and M. Yamaguchi, _A Review of Neutrino Decoupling from the Early Universe to the Current Universe_ , _Universe_ 8 (2022) 552 [2210.10307].
  * [112] T. Fischer, S. Chakraborty, M. Giannotti, A. Mirizzi, A. Payez and A. Ringwald, _Probing axions with the neutrino signal from the next galactic supernova_ , _Phys. Rev. D_ 94 (2016) 085012 [1605.08780].
  * [113] H.K. Dreiner, C. Hanhart, U. Langenfeld and D.R. Phillips, _Supernovae and light neutralinos: SN1987A bounds on supersymmetry revisited_ , _Phys. Rev. D_ 68 (2003) 055004 [hep-ph/0304289].
  * [114] L. Oberauer, C. Hagner, G. Raffelt and E. Rieger, _Supernova bounds on neutrino radiative decays_ , _Astropart. Phys._ 1 (1993) 377.
  * [115] A. de Gouvêa and A. Kobach, _Global Constraints on a Heavy Neutrino_ , _Phys. Rev. D_ 93 (2016) 033005 [1511.00683].
  * [116] J. Beacham et al., _Physics Beyond Colliders at CERN: Beyond the Standard Model Working Group Report_ , _J. Phys. G_ 47 (2020) 010501 [1901.09966].
  * [117] T2K collaboration, _Search for heavy neutrinos with the T2K near detector ND280_ , _Phys. Rev. D_ 100 (2019) 052006 [1902.07598].
  * [118] NA62 collaboration, _Search for heavy neutral lepton production in $K^{+}$ decays to positrons_, _Phys. Lett. B_ 807 (2020) 135599 [2005.09575].
  * [119] S. Alekhin et al., _A facility to Search for Hidden Particles at the CERN SPS: the SHiP physics case_ , _Rept. Prog. Phys._ 79 (2016) 124201 [1504.04855].
  * [120] A. Boyarsky, O. Ruchayskiy and M. Shaposhnikov, _The Role of sterile neutrinos in cosmology and astrophysics_ , _Ann. Rev. Nucl. Part. Sci._ 59 (2009) 191 [0901.0011].
  * [121] O. Ruchayskiy and A. Ivashko, _Restrictions on the lifetime of sterile neutrinos from primordial nucleosynthesis_ , _JCAP_ 10 (2012) 014 [1202.2841].
  * [122] N. Sabti, A. Magalich and A. Filimonova, _An Extended Analysis of Heavy Neutral Leptons during Big Bang Nucleosynthesis_ , _JCAP_ 11 (2020) 056 [2006.07387].
  * [123] S. Hannestad and J. Madsen, _Neutrino decoupling in the early universe_ , _Phys. Rev. D_ 52 (1995) 1764 [astro-ph/9506015].
  * [124] A.D. Dolgov, S.H. Hansen and D.V. Semikoz, _Nonequilibrium corrections to the spectra of massless neutrinos in the early universe_ , _Nucl. Phys. B_ 503 (1997) 426 [hep-ph/9703315].
  * [125] D. Gorbunov and M. Shaposhnikov, _How to find neutral leptons of the $\nu$MSM?_, _JHEP_ 10 (2007) 015 [0705.1729].
  * [126] A. Atre, T. Han, S. Pascoli and B. Zhang, _The Search for Heavy Majorana Neutrinos_ , _JHEP_ 05 (2009) 030 [0901.3589].
  * [127] J.C. Helo, S. Kovalenko and I. Schmidt, _Sterile neutrinos in lepton number and lepton flavor violating decays_ , _Nucl. Phys. B_ 853 (2011) 80 [1005.1607].
  * [128] K. Bondarenko, A. Boyarsky, D. Gorbunov and O. Ruchayskiy, _Phenomenology of GeV-scale Heavy Neutral Leptons_ , _JHEP_ 11 (2018) 032 [1805.08567].
  * [129] P. Coloma, E. Fernández-Martínez, M. González-López, J. Hernández-García and Z. Pavlovic, _GeV-scale neutrinos: interactions with mesons and DUNE sensitivity_ , _Eur. Phys. J. C_ 81 (2021) 78 [2007.03701].
  * [130] S.W. Bruenn, _Stellar core collapse: Numerical model and infall epoch_ , _Astrophys. J. Suppl._ 58 (1985) 771.
  * [131] A. Kamada and H.-B. Yu, _Coherent Propagation of PeV Neutrinos and the Dip in the Neutrino Spectrum at IceCube_ , _Phys. Rev. D_ 92 (2015) 113004 [1504.00711].
  * [132] B. Holdom, _Two U(1)’s and Epsilon Charge Shifts_ , _Phys. Lett. B_ 166 (1986) 196.
  * [133] IMB collaboration, _Angular Distribution of Events From Sn1987a_ , _Phys. Rev. D_ 37 (1988) 3361.
  * [134] A. Mirizzi, G.G. Raffelt and P.D. Serpico, _Earth matter effects in supernova neutrinos: Optimal detector locations_ , _JCAP_ 05 (2006) 012 [astro-ph/0604300].
  * [135] B. Jegerlehner, F. Neubig and G. Raffelt, _Neutrino oscillations and the supernova SN1987A signal_ , _Phys. Rev. D_ 54 (1996) 1194 [astro-ph/9601111].
  * [136] Super-Kamiokande collaboration, _Solar Neutrino Measurements in Super-Kamiokande-IV_ , _Phys. Rev. D_ 94 (2016) 052010 [1606.07538].
  * [137] A. Mirizzi and G.G. Raffelt, _New analysis of the sn 1987a neutrinos with a flexible spectral shape_ , _Phys. Rev. D_ 72 (2005) 063001 [astro-ph/0508612].
  * [138] D.F.G. Fiorillo, M. Heinlein, H.-T. Janka, G. Raffelt and E. Vitagliano, _Supernova Simulations Confront SN 1987A Neutrinos_ , 2308.01403.
  * [139] L. Mastrototaro, P.D. Serpico, A. Mirizzi and N. Saviano, _Massive sterile neutrinos in the early Universe: From thermal decoupling to cosmological constraints_ , _Phys. Rev. D_ 104 (2021) 016026 [2104.11752].
  * [140] P.D. Bolton, F.F. Deppisch and P.S. Bhupal Dev, _Neutrinoless double beta decay versus other probes of heavy sterile neutrinos_ , _JHEP_ 03 (2020) 170 [1912.03058].
  * [141] G.B. Gelmini, M. Kawasaki, A. Kusenko, K. Murai and V. Takhistov, _Big Bang Nucleosynthesis constraints on sterile neutrino and lepton asymmetry of the Universe_ , _JCAP_ 09 (2020) 051 [2005.06721].
  * [142] G. Gelmini, S. Palomares-Ruiz and S. Pascoli, _Low reheating temperature and the visible sterile neutrino_ , _Phys. Rev. Lett._ 93 (2004) 081302 [astro-ph/0403323].
  * [143] G. Gelmini, E. Osoba, S. Palomares-Ruiz and S. Pascoli, _MeV sterile neutrinos in low reheating temperature cosmological scenarios_ , _JCAP_ 10 (2008) 029 [0803.2735].
  * [144] S. Sandner, M. Escudero and S.J. Witte, _Precision CMB constraints on eV-scale bosons coupled to neutrinos_ , _Eur. Phys. J. C_ 83 (2023) 709 [2305.01692].
  * [145] M. Escudero and S.J. Witte, _A CMB search for the neutrino mass mechanism and its relation to the Hubble tension_ , _Eur. Phys. J. C_ 80 (2020) 294 [1909.04044].
  * [146] H. An, M. Pospelov and J. Pradler, _New stellar constraints on dark photons_ , _Phys. Lett. B_ 725 (2013) 190 [1302.3884].
  * [147] E. Hardy and R. Lasenby, _Stellar cooling bounds on new light particles: plasma mixing effects_ , _JHEP_ 02 (2017) 033 [1611.05852].
  * [148] J. Heeck, _Unbroken B – L symmetry_ , _Phys. Lett. B_ 739 (2014) 256 [1408.6845].
  * [149] H. An, M. Pospelov, J. Pradler and A. Ritz, _Direct Detection Constraints on Dark Photon Dark Matter_ , _Phys. Lett. B_ 747 (2015) 331 [1412.8378].
  * [150] D.K. Hong, C.S. Shin and S. Yun, _Cooling of young neutron stars and dark gauge bosons_ , _Phys. Rev. D_ 103 (2021) 123031 [2012.05427].
  * [151] D. Kazanas, R.N. Mohapatra, S. Nussinov, V.L. Teplitz and Y. Zhang, _Supernova Bounds on the Dark Photon Using its Electromagnetic Decay_ , _Nucl. Phys. B_ 890 (2014) 17 [1410.0221].
  * [152] M. Archidiacono and S. Hannestad, _Updated constraints on non-standard neutrino interactions from Planck_ , _JCAP_ 07 (2014) 046 [1311.3873].
  * [153] G.-y. Huang, T. Ohlsson and S. Zhou, _Observational Constraints on Secret Neutrino Interactions from Big Bang Nucleosynthesis_ , _Phys. Rev. D_ 97 (2018) 075009 [1712.04792].