compat=1.0.0

UMN-TH-4217/23

FTPI-MINN-23-10

Heavy Neutral Leptons
from Stopped Muons and Pions

Yohei Ema,a,b Zhen Liu,b Kun-Feng Lyu,b Maxim Pospelova,b

a | _William I. Fine Theoretical Physics Institute, School of Physics and Astronomy,_
---|---
| _University of Minnesota, Minneapolis, MN 55455, USA_
b | _School of Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA_

Stopped muons, which are generic in pion-at-rest experiments, can shed light
on heavy neutral leptons (HNLs) in unexplored parameter spaces. If the HNL is
lighter than the muon, the HNL can be produced from decays of muons and pions.
The HNL will travel from the production location and decay into visible
Standard Model (SM) modes, leaving signals inside downstream detectors. We
find that in the case that the HNL dominantly mixes with muon neutrinos, the
LSND constraint on the mixing angle squared is stronger than all the previous
constraints by more than an order of magnitude. In this study, we recast the
LSND measurement of the $\nu-e$ scattering. Future experiments such as PIP2-BD
could further improve the sensitivity, provided they can distinguish the HNL
events from backgrounds induced by the SM neutrinos.

###### Contents

  1. 1 Introduction
  2. 2 Preliminary
    1. 2.1 HNL production rate
    2. 2.2 HNL decay rate
  3. 3 Neutrino experiments with stopped muon
    1. 3.1 New constraint: LSND
    2. 3.2 Future sensitivitity: PIP2-BD
  4. 4 Summary

##  1 Introduction

The discovery of neutrino oscillation stands out as clear evidence of physics
beyond the Standard Model (SM) [1, 2, 3, 4]. Neutrinos are exactly massless
within the SM, and hence non-zero neutrino masses require an extension of the
SM in the neutrino sector.

Among possible extensions, the introduction of heavy neutral leptons (HNLs) is
particularly interesting due to its minimality (see e.g., Refs [5, 6, 7, 8, 9]
for representative reviews of the HNL). The HNL can couple to the SM lepton
and Higgs doublets at the renormalizable level. After the electroweak symmetry
breaking, the HNL interacts with the SM particles through mixing with the SM
neutrinos. The phenomenology of the HNL can be diverse, depending on the mass
of the HNL. If the HNL is light, $m_{N}\lesssim 1\,\mathrm{eV}$ where $m_{N}$
is the HNL mass, the HNL can participate in the neutrino oscillation, which
may be relevant for, _e.g._ , the anomalies observed at LSND [10] and
MiniBooNE [11].

Even if the HNL is heavy enough to be irrelevant for the neutrino oscillation,
it can still leave signals in different types of experiments. In particular,
if the HNL is lighter than the muon, the HNL can be produced from the decay of
muons and pions. Once produced, the HNL can travel downwards and decay into an
$e^{+}e^{-}$ pair and a neutrino, leaving signals inside downstream detectors,
as long as $m_{N}>2m_{e}$ with $m_{e}$ being the electron mass.

Typically, one can search for HNLs in the particle regime in the high-energy
beamdump experiments and as well at colliders. In the regime of HNL mass
around 10-100 MeV, the searches are dominated by HNL from pion decays and kaon
decays. In this paper, we point out that stopped muons, which commonly exist
in the pion-at-rest experiments, can provide leading productions and enables
us to probe new parameter spaces in HNLs. In this paper, we focus on neutrino
experiments with stopped muons and stopped pions as the source of the HNL,
such as LSND and PIP2-BD [12, 13]. In the past, the LSND experiment has been
used to set constraints on dark photons, dark scalars, and light dark matter
(see , e.g., [14, 15, 16, 17]). In the parameter region of our interest, the
HNL decay length is far longer than the typical size of the laboratory. The
small velocity of the HNL from the muon and pion decay-at-rest, in comparison
to that decay-in-flight, therefore enhances the probability of the HNL
decaying inside the detectors, strengthening the sensitivity. Together with
the large accumulation of proton-on-target (POT), we find that the LSND
experiment provides an upper limit on the mixing angle between the muon
neutrino and the HNL which is stronger by more than an order of magnitude than
the previous constraints from T2K [18, 19] and $\mu$BooNE [20, 21].
Interestingly, our LSND constraint still does not utilize all the properties
of the HNL decay, as LSND did not perform a search of the $e^{+}e^{-}$ events
and thus the HNL events are compared with the single $e$ events originating
from the SM neutrinos. Future experiments such as PIP2-BD potentially improve
the sensitivity significantly, provided that they can distinguish the HNL-
induced $e^{+}e^{-}$ events from backgrounds induced by, _e.g._ , the SM
neutrinos. While this part of the parameter space is disfavored by big bang
nucleosynthesis (BBN) in the minimal model of HNLs, a trivial model
modification, such as dark decay channels can be introduced to render this
model safe from the cosmological constraints.

The rest of this paper is organized as follows. In Sec. 2, we summarize the
production and decay rates of the HNL. This section provides basic formulas
that we use to derive the sensitivity in the subsequent section. We consider
both the muon mixing and electron mixing cases, with an emphasis on the former
case. In Sec. 3, we derive the current constraint on the mixing angle from
LSND and the future sensitivity of PIP2-BD. We summarize our results in Sec.
4.

##  2 Preliminary

In this section, we summarize the basic ingredients necessary for our
analysis. For definiteness, we consider one (pseudo-)Dirac HNL which
dominantly mixes with either muon or electron neutrinos. The relevant part of
the Lagrangian is given by

$\displaystyle\mathcal{L}$
$\displaystyle=\bar{N}\left(i\not{\partial}-m_{N}\right)N-\frac{g}{\sqrt{2}}U_{lN}^{*}\bar{l}\not{W}^{-}N-\frac{g}{2\cos\theta_{W}}U_{lN}^{*}\bar{\nu}_{l}\not{Z}N+(\mathrm{h.c.}),$
(2.1)

where $N$ is the HNL field with its mass $m_{N}$, $g$ is the SU(2) gauge
coupling, $\theta_{W}$ is the weak mixing angle, and $U_{lN}$ is the mixing
angle between the HNL and the SM neutrino. We mainly focus on the case
$m_{N}<m_{\mu}$, where $m_{\mu}$ is the muon mass, so that the HNL can be
produced by the decay of (stopped) muons and pions. We consider both the muon
mixing case $l=\mu$ and the electron mixing case $l=e$, with a particular
interest in the former case. Throughout this paper, we ignore the electron and
neutrino masses as the energy scale of our interest is above $\sim
10\,\mathrm{MeV}$.

Throughout this paper, we assume (pseudo-)Dirac mass $m_{N}$. While in
principle, both Dirac and Majorana masses are possible, the (pseudo-)Dirac
mass is easier to reconcile with neutrino phenomenology. In the inverse seesaw
models, the HNLs are dominated by their Dirac mass and receive a small
correction from their Majorana mass terms \natural1\natural1$\natural$1hence
dubbed pseudo-Dirac.. A single generation with Majorana mass would imply
$|U_{lN}|^{2}m_{N}$ contributions to the active neutrino masses, which is
above the neutrino mass limits for the interesting values of these parameters,
implying that the Majorana mass $m_{N}$ would require an additional finely
tuned contribution to $m_{\nu}$. In general multi-generation HNLs with all
mixings with different neutrinos turned on, an inverse seesaw-like mass
spectrum is generically found to be consistent with the SM neutrino mass
considerations, see, e.g., discussions in review [8].

###  2.1 HNL production rate

We first discuss the HNL production rates from the decay of muons and pions.
The muon mixing case and electron mixing case are studied separately.

#### Muon mixing case

If $m_{N}<m_{\mu}$, the HNL is predominantly produced from the decay of muons
at neutrino experiments with stopped muons such as LSND. The amplitude is
diagrammatically given by

$\displaystyle i\mathcal{M}(\mu\to e\nu_{e}N)=\leavevmode\hbox to0pt{\vbox
to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }
\feynman[inline = (base.b), horizontal=b to d]
\vertex[label=180:${\scriptstyle\mu}$](a); \vertex[right = of a] (b);
\vertex[above right = of b, label=360:${\scriptstyle N}$] (c); \vertex[below
right = of b] (d); \vertex[above right = of d, label=360:${\scriptstyle e}$]
(e); \vertex[below right = of d, label=360:${\scriptstyle\nu_{e}}$] (f);
\diagram*{ (a) -- [fermion] (b) -- [fermion] (c), (b) -- [photon] (d), (f) --
[fermion] (d) -- [fermion] (e), }; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys<EMAIL_ADDRESS>(2.2)

The decay rate is expressed as

$\displaystyle\Gamma(\mu\to
e\nu_{e}N)=\int_{m_{N}}^{(m_{\mu}^{2}+m_{N}^{2})/2m_{\mu}}dE_{N}\frac{d\Gamma(\mu\to
Ne\nu_{e})}{dE_{N}},$ (2.3)

with the differential decay rate given by

$\displaystyle\frac{d\Gamma(\mu\to
e\nu_{e}N)}{dE_{N}}=\frac{G_{F}^{2}\left|{U_{\mu
N}}\right|^{2}}{12\pi^{3}}\,\left(3E_{N}(m_{\mu}^{2}+m_{N}^{2})-4m_{\mu}E_{N}^{2}-2m_{\mu}m_{N}^{2}\right)\sqrt{E_{N}^{2}-m_{N}^{2}},$
(2.4)

where $G_{F}$ is the Fermi constant. If we set $m_{N}=0$, this formula
correctly reproduces the well-known muon decay rate up to the factor $|U_{\mu
N}|^{2}$.

Although subdominant, if $m_{N}<m_{\pi}-m_{\mu}$ with $m_{\pi}$ the pion mass,
the HNL can also be produced from the decay of pions, whose diagram is given
by

$\displaystyle i\mathcal{M}(\pi\to\mu N)$ $\displaystyle=\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }
\feynman[inline = (base.b), horizontal=b to d]
\vertex[label=180:${\scriptstyle\pi}$](a); \vertex[right = of a] (b);
\vertex[above right = of b, label=360:${\scriptstyle N}$] (c); \vertex[below
right = of b, label=360:${\scriptstyle\mu}$] (d); \diagram*{ (a) -- [scalar]
(b), (c) -- [fermion] (b) -- [fermion] (d), };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys<EMAIL_ADDRESS>(2.5)

In this case, the final-state HNL is monochromatic, and the differential decay
rate is given by

$\displaystyle\frac{d\Gamma(\pi\to\mu N)}{dE_{N}}$
$\displaystyle=\frac{G_{F}^{2}f_{\pi}^{2}\left|{V_{ud}}\right|^{2}\left|{U_{\mu
N}}\right|^{2}}{8\pi
m_{\pi}^{3}}\left((m_{\mu}^{2}+m_{N}^{2})m_{\pi}^{2}-(m_{\mu}^{2}+m_{N}^{2})^{2}+4m_{\mu}^{2}m_{N}^{2}\right)$
$\displaystyle\times\sqrt{m_{\pi}^{4}-2(m_{\mu}^{2}+m_{N}^{2})m_{\pi}^{2}+(m_{\mu}^{2}-m_{N}^{2})^{2}}\times\delta\left(E_{N}-\frac{m_{\pi}^{2}+m_{N}^{2}-m_{\mu}^{2}}{2m_{\pi}}\right),$
(2.6)

where $V_{ud}$ is the $ud$-component of the CKM matrix and $f_{\pi}$ is the
pion decay constant.

#### Electron mixing case

As in the muon mixing case, the HNL can be produced from the muon decay if
$m_{N}<m_{\mu}$, and the relevant amplitude is given diagrammatically as

$\displaystyle i\mathcal{M}(\mu\to e\nu_{\mu}N)=\leavevmode\hbox to0pt{\vbox
to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }
\feynman[inline = (base.b), horizontal=b to d]
\vertex[label=180:${\scriptstyle\mu}$](a); \vertex[right = of a] (b);
\vertex[above right = of b, label=360:${\scriptstyle\nu_{\mu}}$] (c);
\vertex[below right = of b] (d); \vertex[above right = of d,
label=360:${\scriptstyle e}$] (e); \vertex[below right = of d,
label=360:${\scriptstyle N}$] (f); \diagram*{ (a) -- [fermion] (b) --
[fermion] (c), (b) -- [photon] (d), (f) -- [fermion] (d) -- [fermion] (e), };
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys<EMAIL_ADDRESS>(2.7)

The decay rate is given by

$\displaystyle\Gamma(\mu\to
e\nu_{\mu}N)=\int_{m_{N}}^{(m_{\mu}^{2}+m_{N}^{2})/2m_{\mu}}dE_{N}\frac{d\Gamma(\mu\to
e\nu_{\mu}N)}{dE_{N}},$ (2.8)

with the differential rate given by

$\displaystyle\frac{d\Gamma(\mu\to e\nu_{\mu}N)}{dE_{N}}$
$\displaystyle=\frac{G_{F}^{2}\left|{U_{eN}}\right|^{2}}{2\pi^{3}}E_{N}(m_{\mu}^{2}+m_{N}^{2}-2m_{\mu}E_{N})\sqrt{E_{N}^{2}-m_{N}^{2}}.$
(2.9)

This formula again correctly reproduces the well-known muon decay rate in the
limit $m_{N}=0$.

In the electron mixing case, the pion produces the HNL if $m_{N}<m_{\pi}$. The
contribution of the pion decay is important at the relatively high $m_{N}$
region. This is in contrast to the muon mixing case, where the decay is
kinematically allowed only for $m_{N}<m_{\pi}-m_{\mu}$. The decay rate of pion
is obtained simply by replacing $\mu\to e$ in the muon mixing case as

$\displaystyle\frac{d\Gamma(\pi\to eN)}{dE_{N}}$
$\displaystyle=\frac{G_{F}^{2}f_{\pi}^{2}\left|{V_{ud}}\right|^{2}\left|{U_{eN}}\right|^{2}m_{N}^{2}}{8\pi
m_{\pi}^{3}}\left(m_{\pi}^{2}-m_{N}^{2}\right)^{2}\times\delta\left(E_{N}-\frac{m_{\pi}^{2}+m_{N}^{2}}{2m_{\pi}}\right).$
(2.10)

Note that the chirality flip is supplied by the HNL and this rate is not
suppressed by the electron mass.

###  2.2 HNL decay rate

We next discuss the decay rate of the HNL. We first note that, in the case of
our interest, the total decay rate of the HNL $\Gamma_{N}$ into the SM
particles is estimated as

$\displaystyle\Gamma_{N}\sim\Gamma_{\mu}\left(\frac{m_{N}}{m_{\mu}}\right)^{5}\left|{U_{lN}}\right|^{2},$
(2.11)

and thus the lifetime is estimated as

$\displaystyle c\tau_{N}\sim
10^{8}\,\mathrm{m}\times\left(\frac{m_{\mu}}{m_{N}}\right)^{5}\left(\frac{10^{-6}}{|{U_{lN}}|^{2}}\right).$
(2.12)

Thus $N$ is sufficiently long-lived at the laboratory scale and only a small
fraction of the HNL decays within the laboratory, even after including the
velocity of the HNL (smaller than the speed of light $c$). Therefore we focus
on the partial decay rate $\Gamma(N\to e^{+}e^{-}\nu)$ as the other decay
modes with only neutrinos in the final state are simply invisible.

#### Muon mixing case

In the muon mixing case, the HNL decays via the neutral current. The relevant
diagram is given by

$\displaystyle i\mathcal{M}(N\to e^{+}e^{-}\nu_{\mu})$
$\displaystyle=\leavevmode\hbox to0pt{\vbox
to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }
\feynman[inline = (base.b), horizontal=b to d]
\vertex[label=180:${\scriptstyle N,\,p_{N}}$](a); \vertex[right = of a] (b);
\vertex[above right = of b, label=360:${\scriptstyle\nu_{\mu},\,p_{\nu}}$]
(c); \vertex[below right = of b] (d); \vertex[above right = of d,
label=360:${\scriptstyle e^{-},\,p_{2}}$] (e); \vertex[below right = of d,
label=360:${\scriptstyle e^{+},\,p_{1}}$] (f); \diagram*{ (a) -- [fermion] (b)
-- [fermion] (c), (b) -- [photon] (d), (f) -- [fermion] (d) -- [fermion] (e),
}; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys<EMAIL_ADDRESS>(2.13)

The spin-averaged matrix element square is computed as

$\displaystyle\left|{\mathcal{M}(N\to
e^{+}e^{-}\nu_{\mu})}\right|^{2}=16G_{F}^{2}\left|{U_{\mu
N}}\right|^{2}\left[4\sin^{4}\theta_{W}(p_{2}\cdot p_{N})(p_{1}\cdot
p_{\nu})+\left(1-2\sin^{2}\theta_{W}\right)^{2}(p_{2}\cdot p_{\nu})(p_{1}\cdot
p_{N})\right],$ (2.14)

and the decay rate is given by

$\displaystyle\Gamma(N\to e^{+}e^{-}\nu_{\mu})$
$\displaystyle=\frac{1}{(2\pi)^{5}}\frac{1}{64m_{N}^{3}}\int_{0}^{m_{N}^{2}}dm_{2\nu}^{2}\int_{0}^{m_{N}^{2}-m_{2\nu}^{2}}dm_{1\nu}^{2}\int_{0}^{2\pi}d\psi\int_{-1}^{1}d\cos\theta\int_{0}^{2\pi}d\phi\left|{\mathcal{M}}\right|^{2}.$
(2.15)

Here we have defined

$\displaystyle
m_{1\nu}^{2}=(p_{1}+p_{\nu})^{2}=m_{N}^{2}-2m_{N}E_{-}^{N},\quad
m_{2\nu}^{2}=(p_{2}+p_{\nu})^{2}=m_{N}^{2}-2m_{N}E_{+}^{N},$ (2.16)

where $E_{\pm}^{N}$ is the energy of $e^{\pm}$ in the $N$-rest frame, and
$(\psi,\theta,\phi)$ are the Euler angle of the orientation of the final
particles (they are on a single plane) relative to the initial $N$. Note that
the angle between $\vec{p}_{1}$ and $\vec{p}_{2}$ in the $N$-rest frame is
fixed once we fix $E_{\pm}^{N}$ due to the energy-momentum conservation as

$\displaystyle\frac{\vec{p}_{1}\cdot\vec{p}_{2}}{|\vec{p}_{1}||\vec{p}_{2}|}=1-\frac{m_{N}}{E_{+}^{N}E_{-}^{N}}\left(E_{+}^{N}+E_{-}^{N}-\frac{m_{N}}{2}\right).$
(2.17)

Later we use this expression for the Monte Carlo (MC) simulation of the cut
efficiency at the LSND experiment. The total decay rate is easily obtained
from this expression as

$\displaystyle\Gamma(N\to
e^{+}e^{-}\nu_{\mu})=\frac{G_{F}^{2}m_{N}^{5}}{768\pi^{3}}\left|{U_{\mu
N}}\right|^{2}\left(1-4\sin^{2}\theta_{W}+8\sin^{4}\theta_{W}\right).$ (2.18)

This correctly reproduces the parameter dependence of Eq. (2.11).

#### Electron mixing case

In the electron mixing case, the HNL can decay into $e^{+}e^{-}\nu_{e}$ both
through the neutral and charged currents. The amplitude is given by

$\displaystyle i\mathcal{M}(N\to e^{+}e^{-}\nu_{e})$
$\displaystyle=\leavevmode\hbox to0pt{\vbox
to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }
\feynman[inline = (base.b), horizontal=b to d]
\vertex[label=180:${\scriptstyle N,\,p_{N}}$](a); \vertex[right = of a] (b);
\vertex[above right = of b, label=360:${\scriptstyle\nu_{e},\,p_{\nu}}$] (c);
\vertex[below right = of b] (d); \vertex[above right = of d,
label=360:${\scriptstyle e^{-},\,p_{2}}$] (e); \vertex[below right = of d,
label=360:${\scriptstyle e^{+},\,p_{1}}$] (f); \diagram*{ (a) -- [fermion] (b)
-- [fermion] (c), (b) -- [photon] (d), (f) -- [fermion] (d) -- [fermion] (e),
}; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}+\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{
}\nullfont\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }
\feynman[inline = (base.b), horizontal=b to d]
\vertex[label=180:${\scriptstyle N,\,p_{N}}$](a); \vertex[right = of a] (b);
\vertex[above right = of b, label=360:${\scriptstyle e^{-},\,p_{2}}$] (c);
\vertex[below right = of b] (d); \vertex[above right = of d,
label=360:${\scriptstyle\nu_{e},\,p_{\nu}}$] (e); \vertex[below right = of d,
label=360:${\scriptstyle e^{+},\,p_{1}}$] (f); \diagram*{ (a) -- [fermion] (b)
-- [fermion] (c), (b) -- [photon] (d), (f) -- [fermion] (d) -- [fermion] (e),
}; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys<EMAIL_ADDRESS>(2.19)

We may note in passing the relative minus sign between the two diagrams from
the exchange of the positions of $\nu_{e}$ and $e^{-}$. The amplitude square
after the spin average is computed as

$\displaystyle\left|{\mathcal{M}(N\to
e^{+}e^{-}\nu_{e})}\right|^{2}=\frac{g^{4}}{2m_{W}^{4}}\left|{U_{eN}}\right|^{2}\left[4\sin^{4}\theta_{W}(p_{2}\cdot
p_{N})(p_{1}\cdot p_{\nu})+\left(1+2\sin^{2}\theta_{W}\right)^{2}(p_{2}\cdot
p_{\nu})(p_{1}\cdot p_{N})\right].$ (2.20)

The rest of the computation proceeds in the same way as the muon mixing case.
The decay rate is given by

$\displaystyle\Gamma(N\to e^{+}e^{-}\nu_{e})$
$\displaystyle=\frac{1}{(2\pi)^{5}}\frac{1}{64m_{N}^{3}}\int_{0}^{m_{N}^{2}}dm_{2\nu}^{2}\int_{0}^{m_{N}^{2}-m_{2\nu}^{2}}dm_{1\nu}^{2}\int_{0}^{2\pi}d\psi\int_{-1}^{1}d\cos\theta\int_{0}^{2\pi}d\phi\left|{\mathcal{M}}\right|^{2},$
(2.21)

and after performing all the integrations, it is given by

$\displaystyle\Gamma(N\to e^{+}e^{-}\nu_{e})$
$\displaystyle=\frac{G_{F}^{2}m_{N}^{5}}{768\pi^{3}}\left|{U_{eN}}\right|^{2}\left(1+4\sin^{2}\theta_{W}+8\sin^{4}\theta_{W}\right).$
(2.22)

In the electron mixing case, if $m_{\mu}<m_{N}<m_{\pi}$, the HNL produced from
the pion decay can also decay as $N\to\mu e\nu_{e}$. However, as we will see,
the LSND constraint for this mass range in the electron mixing case is anyway
weaker than the other constraints, and thus we may ignore this decay channel
in the rest of the analysis.

##  3 Neutrino experiments with stopped muon

Armed with the basic formulas, we now derive the current constraint and future
sensitivity of neutrino experiments with stopped muon on the HNL mixing
parameter $U_{lN}$. Our signal is $e^{+}e^{-}$ produced from the HNL decay.
The total event number of such a decay inside a detector is estimated as

$\displaystyle N^{(i)}_{ee}$
$\displaystyle=N_{i}\times\epsilon_{\mathrm{det}}\times\frac{1}{\Gamma_{i}}\int
dE_{N}\frac{d\Gamma(i\to fN)}{dE_{N}}\times\frac{L_{\mathrm{det}}}{\gamma\beta
c\tau_{N\to ee\nu}},$ (3.1)

where $i=\mu$ if the HNL is produced from the stopped muon and $i=\pi$ if
produced from the stopped pion, and $f$ represents SM particles depending on
the decay channel of $i$. In this expression, $N_{i}$ is the total number of
$i$ that decays at the HNL production point, $\epsilon_{\mathrm{det}}$ is the
probability that a single HNL produced from $i$ passes through the detector,
and $\Gamma_{i}$ is the total decay width of $i$. The relativistic factor
$\gamma\beta=(E_{N}^{2}-m_{N}^{2})^{1/2}/m_{N}$ with $E_{N}$ being the energy
of the HNL, $L_{\mathrm{det}}$ denotes the length of the detector, and
$\tau_{N\to ee\nu}$ is the decay lifetime of the HNL to a pair of $e^{+}e^{-}$
at rest. The signal rate has no $N$ total width $\Gamma_{N}$ dependence so
long as the long lifetime limit is valid, as $1/c\tau_{N}\times{\rm
Br}(N\rightarrow ee\nu)=1/c\tau_{N\rightarrow ee\nu}$. Note that
$\gamma\beta<m_{\mu}/m_{N}$ for the HNL produced from the stopped muon. This
enhances the sensitivity to the HNL compared to neutrino experiments with
pions and muons decay-in-flight, such as the MiniBooNE experiment [22, 23].

Of course, experiments impose cuts on events to eliminate backgrounds, and
hence we compute the efficiency $\epsilon^{(i)}$ that the $e^{+}e^{-}$ event
from the HNL decay passes through the cuts. For this purpose, we perform a
Monte Carlo simulation to generate the distribution of the energy $E_{\pm}$
and the angle $\cos\theta_{\pm}$ \natural2\natural2$\natural$2with respect to
the direction of the production point to the center of the detector of the
$e^{\pm}$ in the laboratory frame. More specifically, we generate the HNL
energy, $E_{N}$, from the differential decay rates of $\mu$ and $\pi$ computed
in Sec. 2.1 , and the energy of $e^{\pm}$ in the $N$-rest frame $E_{\pm}^{N}$,
as well as the Euler angle ($\psi,\theta,\phi$) from Eqs. (2.15) and (2.21).
We then perform the Lorentz transformation to obtain $E_{\pm}$ and
$\cos\theta_{\pm}$ in the laboratory frame. The HNL event number after
imposing the cut is given by

$\displaystyle
N_{ee}=N_{ee}^{(\mu)}\times\epsilon^{(\mu)}+N_{ee}^{(\pi)}\times\epsilon^{(\pi)}.$
(3.2)

In the following, we study the current constraint on $U_{lN}$ from the LSND
experiment [24], and the future sensitivity of the PIP2-BD experiment [12,
13], respectively.

###  3.1 New constraint: LSND

Here we derive the current constraint on the mixing parameter $U_{lN}$ from
the LSND experiment. The Liquid Scintillator Neutrino Detector (LSND)
experiment [24] was located at the Los Alamos Neutron Science Center (LANSCE).
An 800 MeV proton beam is delivered onto the water target (and other materials
such as tungsten), producing a large flux of pions dumped at the copper beam
stop. Most of $\pi^{-}$ are absorbed by the materials and hence the abundance
of $\mu^{-}$ is suppressed. The main source of the neutrino flux is the decay
of stopped $\pi^{+}$ and $\mu^{+}$. The detector is located $\sim
30\,\mathrm{m}$ downstream and is well-shielded. In the case of our interest,
the HNL can be produced from the decay of stopped $\mu^{+}$ and $\pi^{+}$.
Once produced, a small portion of the HNL travels downstream and decays inside
the detector, providing $e^{+}e^{-}$ signatures.

The LSND experiment did not perform any search on the $e^{+}e^{-}$ events.
Therefore, we may use the single $e$ event search in [25] to put a constraint
on the HNL. The number of $\nu_{e}$ that passes through the detector is
$\simeq 3\times 10^{19}$ in total [25]. Since $\nu_{e}$ is dominantly produced
from stopped $\mu^{+}$, which is dominantly produced from stopped $\pi^{+}$,
we can set
$N_{\mu}\times\epsilon_{\mathrm{det}}=N_{\pi}\times\epsilon_{\mathrm{det}}=3\times
10^{19}$. Therefore we obtain the total number of the HNL decay events inside
the detector originating from $\mu^{+}$ as

$\displaystyle N_{ee}^{(\mu)}\simeq\begin{cases}\displaystyle 1.5\times
10^{5}\left(\frac{\left|{U_{\mu
N}}\right|^{2}}{10^{-6}}\right)^{2}\left(\frac{m_{N}}{m_{\mu}}\right)^{6}\left(1-\frac{m_{N}}{m_{\mu}}\right)^{4}\left(1+\frac{4m_{N}}{m_{\mu}}+\frac{m_{N}^{2}}{m_{\mu}^{2}}\right)&:\mu\mathchar
45\relax\mathrm{mixing},\vspace{1mm}\\\ \displaystyle 8.1\times
10^{5}\left(\frac{\left|{U_{eN}}\right|^{2}}{10^{-6}}\right)^{2}\left(\frac{m_{N}}{m_{\mu}}\right)^{6}\left(1-\frac{m_{N}}{m_{\mu}}\right)^{4}\left(1+\frac{4m_{N}}{m_{\mu}}+\frac{m_{N}^{2}}{m_{\mu}^{2}}\right)&:e\mathchar
45\relax\mathrm{mixing},\end{cases}$ (3.3)

and that originating from $\pi^{+}$ as

$\displaystyle N_{ee}^{(\pi)}\simeq\begin{cases}\displaystyle 8.7\times
10^{4}\left(\frac{\left|{U_{\mu
N}}\right|^{2}}{10^{-6}}\right)^{2}\left(\frac{m_{N}}{m_{\mu}}\right)^{6}\frac{m_{\pi}\left((m_{\mu}^{2}+m_{N}^{2})m_{\pi}^{2}-(m_{\mu}^{2}+m_{N}^{2})^{2}+4m_{\mu}^{2}m_{N}^{2}\right)}{m_{\mu}(m_{\pi}^{2}-m_{\mu}^{2})^{2}}&:\mu\mathchar
45\relax\mathrm{mixing},\vspace{1mm}\\\ \displaystyle 4.1\times
10^{5}\left(\frac{\left|{U_{eN}}\right|^{2}}{10^{-6}}\right)^{2}\left(\frac{m_{N}}{m_{\mu}}\right)^{7}\frac{m_{N}m_{\pi}(m_{\pi}^{2}-m_{N}^{2})}{(m_{\pi}^{2}-m_{\mu}^{2})^{2}}&:e\mathchar
45\relax\mathrm{mixing},\end{cases}$ (3.4)

where we take $L_{\mathrm{det}}=7.6\,\mathrm{m}$ following the fiducial cut
specified in [25]. The total number of the HNL decay events is thus sizable,
and we can indeed put a stronger constraint on the mixing parameter than the
other experiments for $m_{N}<m_{\mu}$, especially in the muon mixing case, as
we will see.

Figure 1:  The efficiency of the $e^{+}e^{-}$ events passing the cuts 1 and 2
with the HNL originating from the muon decay. The total efficiency
$\epsilon^{(\mu)}$ is the sum of the cuts 1 and 2. _Left_ : the muon mixing
case. _Right_ : the electron mixing case.

Our signal mimics the single $e$ event in [25] if either (1) one $e^{\pm}$ has
small energy so that it is not detected, while the other $e^{\mp}$ has enough
energy for the detection, or (2) $e^{+}$ and $e^{-}$ are collinear enough so
that they cannot be separated within the detector resolution. Note that
$e^{+}$ and $e^{-}$ are produced at the same time in our case so that it is
excluded by neither past nor future activity cuts imposed in [25]. We thus
impose two alternative cuts on our events as

$\displaystyle\begin{cases}~{}\mathrm{cut}~{}1:~{}18\,\mathrm{MeV}<E_{\pm}<50\,\mathrm{MeV},~{}E_{\mp}<E_{\mathrm{th}},~{}\cos\theta_{\pm}>0.9,~{}\cos\theta_{\mp}<0.9,\vspace{2mm}\\\
~{}\mathrm{cut}~{}2:~{}18\,\mathrm{MeV}<E_{+}+E_{-}<50\,\mathrm{MeV},~{}\cos\theta_{\pm}>0.9.\end{cases}$
(3.5)

Cut 1 captures the $e^{+}e^{-}$ events mimicking single-$e$ events through one
of them being soft, and cut 2 captures when they are being colinear. The
events passing either cut 1 or cut 2 will show up in the $\nu-e$ scattering
data sample of the LSND. Here the cut on the angle originates from the
requirement that the observed event has the direction consistent with the
neutrino from the production point, and $E_{\mathrm{th}}$ is the threshold
energy of the $e^{\pm}$ detection. We take $E_{\mathrm{th}}=15\,\mathrm{MeV}$
as there are huge backgrounds of gamma rays in this energy range from, _e.g._
, Carbon excited states and thus the analysis is expected to be insensitive to
this region. We also show the result with $E_{\mathrm{th}}=10\,\mathrm{MeV}$
as an even more conservative choice of the threshold energy. In Fig. 1, we
show the efficiency of the muon decay channel $\epsilon^{(\mu)}$ computed by
our MC simulation. As one can see, the cut 1 dominates for $m_{N}\gtrsim
20\,\mathrm{MeV}$, while the cut 2 is more important at the lower mass range
as the lighter HNL is more boosted in the forward direction.

Figure 2:  The constraints on $|U_{lN}|^{2}$ from the LSND experiment together
with the previous constraints (all at 90% C.L.) from PSI [26], TRIUMF [27],
PIENU [28, 29, 30], T2K [18, 19], and $\mu$BooNE [20, 21], with the data
adapted from Heavy-Neutrino-Limits [31]. The solid blue line corresponds to
the summation of cut 1 with $E_{\mathrm{th}}=15\,\mathrm{MeV}$ and cut 2,
while the dashed blue line corresponds to that with
$E_{\mathrm{th}}=10\,\mathrm{MeV}$. The future sensitivity of PIP2-BD is also
shown, where the solid (dashed) orange line corresponds to 3 (50) events of
the HNL decay inside the detector, together with expected sensitivities of
DUNE [32] and PIONEER [33, 34]. _Left_ : the muon mixing case. _Right_ : the
electron mixing case. The kink around $m_{N}\simeq 34\,\mathrm{MeV}$ in the
muon mixing case comes from the pion contribution. In the electron mixing
case, the pion contribution dominates over the muon contribution for
$m_{N}\gtrsim 40\,\mathrm{MeV}$.

To derive the constraint, we may follow the procedure used in [25] to put an
upper bound on the neutrino magnetic dipole moments. The LSND experiment
observed in total 242 single $e$ events (including all the $\nu_{e}$,
$\nu_{\mu}$, and $\bar{\nu}_{\mu}$ initiated signals) while 229 events are
expected within the SM. With the systematic error included, we impose at 90%
C.L.

$\displaystyle
N_{ee}=N_{ee}^{(\mu)}\times\epsilon^{(\mu)}+N_{ee}^{(\pi)}\times\epsilon^{(\pi)}<55.$
(3.6)

In Fig. 2, we show the new constraint derived in this study from the LSND
(with the future sensitivity of PIP2-BD, which we describe in the next
section) and the previous constraints from PSI [26], TRIUMF [27], PIENU [28,
29, 30], T2K [18, 19], and $\mu$BooNE [20, 21]. The solid blue line
corresponds to $E_{\mathrm{th}}=15\,\mathrm{MeV}$ while the dashed blue line
does to $10\,\mathrm{MeV}$. As one can see, in the muon mixing case, our
constraint is better than the previous ones by more than one order of
magnitude in terms of $|U_{\mu N}|^{2}$ for $35\,\mathrm{MeV}\lesssim
m_{N}\lesssim 70\,\mathrm{MeV}$. In the electron mixing case, the current
constraint from LSND is weaker than but close to the constraint from PIENU.

Here we comment on the JSNS2 experiment at J-PARC [35, 36]. JSNS2 aims at
measuring the neutrino oscillation $\bar{\nu}_{\mu}\to\bar{\nu}_{e}$ with the
neutrinos produced from stopped pions and muons created by a 3 GeV proton
beam, with the lower duty factor than LSND and the Gd-loaded liquid
scintillator detector. Thus, JSNS2 provides a direct test of the LSND neutrino
oscillation results [10]. In our case, the HNL decay would be
indistinguishable from the single $e$ events (without the final state neutron)
as the JSNS2 detector is not sensitive to the angular separation of $e^{+}$
and $e^{-}$ [37]. The current data is $1.45\times 10^{22}$ POT, 13 % of the
designed total POT [38], compared to $1.8\times 10^{23}$ POT at LSND with a
larger volume of the detector. Still, once enough data is accumulated, it
would be interesting to study the sensitivity of JSNS2 on the HNL. In
particular, JSNS2 can be sensitive to neutrinos from Kaon decay-at-rest, due
to higher proton beam energy, which may allow us to go beyond the muon and
pion mass thresholds.

Before closing this subsection, we briefly comment on cosmological constraints
on the mixing parameter. Indeed, in the parameter region of our interest, the
HNL decay rate into the SM particles is small so that it may decay after the
BBN, disturbing the consistency between the theory and observation [39, 40,
41, 42, 43]. However, this can be avoided if, _e.g._ , the HNL decays into
invisible particles before the BBN. As long as the HNL decay length is longer
than $\sim\mathcal{O}(10\,\mathrm{m})$, this additional decay channel does not
affect our analysis. In this sense, our constraint is independent of the one
from BBN.

###  3.2 Future sensitivitity: PIP2-BD

In this subsection, we derive the future sensitivity of the PIP2-BD experiment
[12, 13]. The Proton Improvement Project II (PIP-II) is a major upgrade of the
accelerator complex at Fermilab to meet the requirement of hosting the Deep
Underground Neutrino Experiment (DUNE) [44, 45]. On top of its main purpose,
PIP-II has the flexibility of supporting multiple experiments, and a beam dump
facility, PIP2-BD, was proposed as a candidate experiment to study sub-GeV
dark sectors.

In the case of our interest, PIP2-BD can probe the HNL in the same way as
LSND. The HNL is produced from stopped pions and muons at the target. We may
assume five years of physics run with the baseline PAR option [12, 13], which
results in an 800 MeV proton beam with $1.2\times 10^{23}$ POT. We take the
formation rate of stopped $\pi^{+}$ (and hence $\mu^{+}$) per proton to be
0.1, as PIP2-BD is designed to use a lighter target that has a larger
formation rate than Hg. For the mercury target, for instance, the COHERENT
experiment reported a formation rate of $(9.0\pm 0.9)\times 10^{-2}$ with
slightly higher proton beam energy [46]. This fixes the parameters in Eq.
(3.1) as $N_{\pi}=N_{\mu}=1.2\times 10^{22}$. We assume that the active volume
of the detector is cylindrical in shape with 4.5 m in height and 4.5 m in
diameter, located 18 m away from the HNL production point. We then estimate
from the solid angle coverage that $\epsilon_{\mathrm{det}}\simeq 3.9\times
10^{-3}$. The total event number of the HNL decay originating from $\mu^{+}$
is thus given as

$\displaystyle N_{ee}^{(\mu)}\simeq\begin{cases}\displaystyle 1.4\times
10^{5}\left(\frac{\left|{U_{\mu
N}}\right|^{2}}{10^{-6}}\right)^{2}\left(\frac{m_{N}}{m_{\mu}}\right)^{6}\left(1-\frac{m_{N}}{m_{\mu}}\right)^{4}\left(1+\frac{4m_{N}}{m_{\mu}}+\frac{m_{N}^{2}}{m_{\mu}^{2}}\right)&:\mu\mathchar
45\relax\mathrm{mixing},\vspace{1mm}\\\ \displaystyle 7.5\times
10^{5}\left(\frac{\left|{U_{eN}}\right|^{2}}{10^{-6}}\right)^{2}\left(\frac{m_{N}}{m_{\mu}}\right)^{6}\left(1-\frac{m_{N}}{m_{\mu}}\right)^{4}\left(1+\frac{4m_{N}}{m_{\mu}}+\frac{m_{N}^{2}}{m_{\mu}^{2}}\right)&:e\mathchar
45\relax\mathrm{mixing},\end{cases}$ (3.7)

and that originating from $\pi^{+}$ as

$\displaystyle N_{ee}^{(\pi)}\simeq\begin{cases}\displaystyle 8.0\times
10^{4}\left(\frac{\left|{U_{\mu
N}}\right|^{2}}{10^{-6}}\right)^{2}\left(\frac{m_{N}}{m_{\mu}}\right)^{6}\frac{m_{\pi}\left((m_{\mu}^{2}+m_{N}^{2})m_{\pi}^{2}-(m_{\mu}^{2}+m_{N}^{2})^{2}+4m_{\mu}^{2}m_{N}^{2}\right)}{m_{\mu}(m_{\pi}^{2}-m_{\mu}^{2})^{2}}&:\mu\mathchar
45\relax\mathrm{mixing},\vspace{1mm}\\\ \displaystyle 3.8\times
10^{5}\left(\frac{\left|{U_{eN}}\right|^{2}}{10^{-6}}\right)^{2}\left(\frac{m_{N}}{m_{\mu}}\right)^{7}\frac{m_{N}m_{\pi}(m_{\pi}^{2}-m_{N}^{2})}{(m_{\pi}^{2}-m_{\mu}^{2})^{2}}&:e\mathchar
45\relax\mathrm{mixing},\end{cases}$ (3.8)

where we take $L_{\mathrm{det}}=4.5\,\mathrm{m}$. Without any dedicated study
on possible backgrounds at this moment, we may simply draw the lines that
correspond to 3 and 50 events of the HNL decay inside the detector, assuming
75 % of event acceptance independent of $m_{N}$ following [13]. In Fig. 2, we
show the future sensitivity of PIP2-BD, together with the expected
sensitivities of DUNE [32] and PIONEER [33, 34]. The solid orange line
corresponds to 3 events while the dashed orange line corresponds to 50 events.
It is clear from the figure that PIP2-BD has the potential of exploring the
parameter region that is not yet and will not be covered by the other
experiments, especially in the muon mixing case. Here the key difference from
LSND is the event acceptance, thus it would be essential to distinguish the
$e^{+}e^{-}$ events from backgrounds, in particular, the single $e$ events
from the ordinary neutrinos, to realize this level of sensitivity.

##  4 Summary

In this paper, we studied the sensitivity of the neutrino experiments with
stopped muons and pions as the source, such as LSND and PIP2-BD, on the HNL.
With the help of the small velocity of the produced HNL as well as a large
number of POT, we find that the LSND constraint on the mixing angle squared is
stronger than the previous constraints by more than one order of magnitude
below the muon mass, if the HNL dominantly mixes with muon neutrinos. Since
LSND does not perform a search for $e^{+}e^{-}$ events, the HNL decay events
are compared with the single $e$ events induced by the SM neutrinos in our
analysis. Therefore, if future experiments such as PIP2-BD can distinguish the
$e^{+}e^{-}$ events from backgrounds induced by the SM neutrinos, the
sensitivity will be significantly improved.

Although we have focused on PIP2-BD as a future experiment, there are other
near-future experiments that are expected to be sensitive to the HNL. For
instance, the Coherent CAPTAIN-Mills (CCM) experiment aims at probing sterile
neutrinos and sub-GeV dark matters produced at the proton beam dump at LANSCE
[47]. Even though the planned number of the total POT as well as the detector
volume is smaller than LSND, CCM may potentially explore a new parameter
region of the HNL, given that the $e^{\pm}$ are efficiently resolved by,
_e.g._ , the Cherenkov light with their liquid Ar detector. The COHERENT
experiment [48, 49, 46] is another possibility once the detector volume is
upgraded to be large enough.

The stopped muon and pion sources have been shown to be sensitive probes of
sub-100-MeV physics for dark photons, dark scalars, light dark matter, and now
for the HNLs. Their sensitivity to axion-like particles (ALPs) could also be
interesting. Our initial estimations show that the LSND detector is capable of
setting competitive bounds on ALP-photon couplings, which also means that the
sensitivity of PIP2-BD to ALPs must also be investigated.

##### Acknowledgements

We would like to thank Drs. M. Hostert, P. Huber, W. Louis, T. Maruyama, and
R. Tayloe for helpful discussions. This work is supported in part by the DOE
grant DE-SC0011842. Z.L. and K.L. were supported in part by DOE grant DE-
SC0022345. The Feynman diagrams in this paper are generated by TikZ-Feynman
[50]. The data associated with the figures in this paper can be accessed via
GitHub.

## References

  * [1] R. Davis, Jr., D. S. Harmer, and K. C. Hoffman, “Search for neutrinos from the sun,” Phys. Rev. Lett. 20 (1968) 1205–1209.
  * [2] Super-Kamiokande Collaboration, Y. Fukuda et al., “Evidence for oscillation of atmospheric neutrinos,” Phys. Rev. Lett. 81 (1998) 1562–1567, arXiv:hep-ex/9807003.
  * [3] SNO Collaboration, Q. R. Ahmad et al., “Direct evidence for neutrino flavor transformation from neutral current interactions in the Sudbury Neutrino Observatory,” Phys. Rev. Lett. 89 (2002) 011301, arXiv:nucl-ex/0204008.
  * [4] 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, arXiv:2007.14792 [hep-ph].
  * [5] A. Boyarsky, O. Ruchayskiy, and M. Shaposhnikov, “The Role of sterile neutrinos in cosmology and astrophysics,” Ann. Rev. Nucl. Part. Sci. 59 (2009) 191–214, arXiv:0901.0011 [hep-ph].
  * [6] M. Drewes, “The Phenomenology of Right Handed Neutrinos,” Int. J. Mod. Phys. E 22 (2013) 1330019, arXiv:1303.6912 [hep-ph].
  * [7] B. Dasgupta and J. Kopp, “Sterile Neutrinos,” Phys. Rept. 928 (2021) 1–63, arXiv:2106.05913 [hep-ph].
  * [8] A. M. Abdullahi et al., “The present and future status of heavy neutral leptons,” J. Phys. G 50 (2023) no. 2, 020501, arXiv:2203.08039 [hep-ph].
  * [9] B. Batell et al., “Dark Sector Studies with Neutrino Beams,” in Snowmass 2021.  7, 2022.  arXiv:2207.06898 [hep-ph].
  * [10] LSND Collaboration, A. Aguilar et al., “Evidence for neutrino oscillations from the observation of $\bar{\nu}_{e}$ appearance in a $\bar{\nu}_{\mu}$ beam,” Phys. Rev. D 64 (2001) 112007, arXiv:hep-ex/0104049.
  * [11] MiniBooNE Collaboration, A. A. Aguilar-Arevalo et al., “Significant Excess of ElectronLike Events in the MiniBooNE Short-Baseline Neutrino Experiment,” Phys. Rev. Lett. 121 (2018) no. 22, 221801, arXiv:1805.12028 [hep-ex].
  * [12] W. Pellico, C. Bhat, J. Eldred, C. Johnstone, J. Johnstone, K. Seiya, C.-Y. Tan, M. Toups, P. deNiverville, and R. Van De Water, “FNAL PIP-II Accumulator Ring,” arXiv:2203.07339 [physics.acc-ph].
  * [13] M. Toups et al., “PIP2-BD: GeV Proton Beam Dump at Fermilab’s PIP-II Linac,” in Snowmass 2021.  3, 2022.  arXiv:2203.08079 [hep-ex].
  * [14] B. Batell, M. Pospelov, and A. Ritz, “Exploring Portals to a Hidden Sector Through Fixed Targets,” Phys. Rev. D 80 (2009) 095024, arXiv:0906.5614 [hep-ph].
  * [15] P. deNiverville, M. Pospelov, and A. Ritz, “Observing a light dark matter beam with neutrino experiments,” Phys. Rev. D 84 (2011) 075020, arXiv:1107.4580 [hep-ph].
  * [16] S. Foroughi-Abari and A. Ritz, “LSND Constraints on the Higgs Portal,” Phys. Rev. D 102 (2020) no. 3, 035015, arXiv:2004.14515 [hep-ph].
  * [17] Y. Ema, Z. Liu, K.-F. Lyu, and M. Pospelov, “Flavor-changing light bosons with accidental longevity,” JHEP 02 (2023) 135, arXiv:2211.00664 [hep-ph].
  * [18] T2K Collaboration, K. Abe et al., “Search for heavy neutrinos with the T2K near detector ND280,” Phys. Rev. D 100 (2019) no. 5, 052006, arXiv:1902.07598 [hep-ex].
  * [19] C. A. Argüelles, N. Foppiani, and M. Hostert, “Heavy neutral leptons below the kaon mass at hodoscopic neutrino detectors,” Phys. Rev. D 105 (2022) no. 9, 095006, arXiv:2109.03831 [hep-ph].
  * [20] MicroBooNE Collaboration, P. Abratenko et al., “Search for a Higgs Portal Scalar Decaying to Electron-Positron Pairs in the MicroBooNE Detector,” Phys. Rev. Lett. 127 (2021) no. 15, 151803, arXiv:2106.00568 [hep-ex].
  * [21] K. J. Kelly and P. A. N. Machado, “MicroBooNE experiment, NuMI absorber, and heavy neutral leptons,” Phys. Rev. D 104 (2021) no. 5, 055015, arXiv:2106.06548 [hep-ph].
  * [22] MiniBooNE Collaboration, A. A. Aguilar-Arevalo et al., “The Neutrino Flux prediction at MiniBooNE,” Phys. Rev. D 79 (2009) 072002, arXiv:0806.1449 [hep-ex].
  * [23] MiniBooNE Collaboration, A. A. Aguilar-Arevalo et al., “The MiniBooNE Detector,” Nucl. Instrum. Meth. A 599 (2009) 28–46, arXiv:0806.4201 [hep-ex].
  * [24] LSND Collaboration, C. Athanassopoulos et al., “The Liquid scintillator neutrino detector and LAMPF neutrino source,” Nucl. Instrum. Meth. A 388 (1997) 149–172, arXiv:nucl-ex/9605002.
  * [25] LSND Collaboration, L. B. Auerbach et al., “Measurement of electron - neutrino - electron elastic scattering,” Phys. Rev. D 63 (2001) 112001, arXiv:hep-ex/0101039.
  * [26] M. Daum, B. Jost, R. M. Marshall, R. C. Minehart, W. A. Stephens, and K. O. H. Ziock, “Search for Admixtures of Massive Neutrinos in the Decay $\pi^{+}\to\mu^{+}$ Neutrino,” Phys. Rev. D 36 (1987) 2624.
  * [27] D. I. Britton et al., “Improved search for massive neutrinos in pi+ —$>$ e+ neutrino decay,” Phys. Rev. D 46 (1992) R885–R887.
  * [28] PIENU Collaboration, A. Aguilar-Arevalo et al., “Improved search for heavy neutrinos in the decay $\pi\rightarrow e\nu$,” Phys. Rev. D 97 (2018) no. 7, 072012, arXiv:1712.03275 [hep-ex].
  * [29] PIENU Collaboration, A. Aguilar-Arevalo et al., “Search for heavy neutrinos in $\pi\to\mu\nu$ decay,” Phys. Lett. B 798 (2019) 134980, arXiv:1904.03269 [hep-ex].
  * [30] D. A. Bryman and R. Shrock, “Constraints on Sterile Neutrinos in the MeV to GeV Mass Range,” Phys. Rev. D 100 (2019) 073011, arXiv:1909.11198 [hep-ph].
  * [31] E. Fernández-Martínez, M. González-López, J. Hernández-García, M. Hostert, and J. López-Pavón, “Effective portals to heavy neutral leptons,” arXiv:2304.06772 [hep-ph].
  * [32] J. M. Berryman, A. de Gouvea, P. J. Fox, B. J. Kayser, K. J. Kelly, and J. L. Raaf, “Searches for Decays of New Particles in the DUNE Multi-Purpose Near Detector,” JHEP 02 (2020) 174, arXiv:1912.07622 [hep-ph].
  * [33] PIONEER Collaboration, W. Altmannshofer et al., “PIONEER: Studies of Rare Pion Decays,” arXiv:2203.01981 [hep-ex].
  * [34] PIONEER Collaboration, W. Altmannshofer et al., “Testing Lepton Flavor Universality and CKM Unitarity with Rare Pion Decays in the PIONEER experiment,” in Snowmass 2021.  3, 2022.  arXiv:2203.05505 [hep-ex].
  * [35] S. Ajimura et al., “Technical Design Report (TDR): Searching for a Sterile Neutrino at J-PARC MLF (E56, JSNS2),” arXiv:1705.08629 [physics.ins-det].
  * [36] JSNS2 Collaboration, S. Ajimura et al., “The JSNS2 detector,” Nucl. Instrum. Meth. A 1014 (2021) 165742, arXiv:2104.13169 [physics.ins-det].
  * [37] J. R. Jordan, Y. Kahn, G. Krnjaic, M. Moschella, and J. Spitz, “Signatures of Pseudo-Dirac Dark Matter at High-Intensity Neutrino Experiments,” Phys. Rev. D 98 (2018) no. 7, 075020, arXiv:1806.05185 [hep-ph].
  * [38] JSNS2, JSNS2 -II, Collaboration, T. Maruyama, “The status of JSNS2 and JSNS2-II,” PoS NuFact2021 (2022) 159.
  * [39] S. Sarkar, “Big bang nucleosynthesis and physics beyond the standard model,” Rept. Prog. Phys. 59 (1996) 1493–1610, arXiv:hep-ph/9602260.
  * [40] A. D. Dolgov, S. H. Hansen, G. Raffelt, and D. V. Semikoz, “Heavy sterile neutrinos: Bounds from big bang nucleosynthesis and SN1987A,” Nucl. Phys. B 590 (2000) 562–574, arXiv:hep-ph/0008138.
  * [41] O. Ruchayskiy and A. Ivashko, “Restrictions on the lifetime of sterile neutrinos from primordial nucleosynthesis,” JCAP 10 (2012) 014, arXiv:1202.2841 [hep-ph].
  * [42] M. Hufnagel, K. Schmidt-Hoberg, and S. Wild, “BBN constraints on MeV-scale dark sectors. Part I. Sterile decays,” JCAP 02 (2018) 044, arXiv:1712.03972 [hep-ph].
  * [43] A. Boyarsky, M. Ovchynnikov, O. Ruchayskiy, and V. Syvolap, “Improved big bang nucleosynthesis constraints on heavy neutral leptons,” Phys. Rev. D 104 (2021) no. 2, 023517, arXiv:2008.00749 [hep-ph].
  * [44] DUNE Collaboration, B. Abi et al., “Deep Underground Neutrino Experiment (DUNE), Far Detector Technical Design Report, Volume I Introduction to DUNE,” JINST 15 (2020) no. 08, T08008, arXiv:2002.02967 [physics.ins-det].
  * [45] DUNE Collaboration, B. Abi et al., “Deep Underground Neutrino Experiment (DUNE), Far Detector Technical Design Report, Volume II: DUNE Physics,” arXiv:2002.03005 [hep-ex].
  * [46] COHERENT Collaboration, D. Akimov et al., “First Measurement of Coherent Elastic Neutrino-Nucleus Scattering on Argon,” Phys. Rev. Lett. 126 (2021) no. 1, 012002, arXiv:2003.10630 [nucl-ex].
  * [47] CCM Collaboration, A. A. Aguilar-Arevalo et al., “First dark matter search results from Coherent CAPTAIN-Mills,” Phys. Rev. D 106 (2022) no. 1, 012001, arXiv:2105.14020 [hep-ex].
  * [48] COHERENT Collaboration, D. Akimov et al., “The COHERENT Experiment at the Spallation Neutron Source,” arXiv:1509.08702 [physics.ins-det].
  * [49] COHERENT Collaboration, D. Akimov et al., “Observation of Coherent Elastic Neutrino-Nucleus Scattering,” Science 357 (2017) no. 6356, 1123–1126, arXiv:1708.01294 [nucl-ex].
  * [50] J. Ellis, “TikZ-Feynman: Feynman diagrams with TikZ,” Comput. Phys. Commun. 210 (2017) 103–123, arXiv:1601.05437 [hep-ph].