# Experimental Searches For Heavy Neutral Leptons

Sophie C. Middleton  California Institute of Technology, Pasadena, California
91125, USA

###### Abstract

on behalf of the BABAR Collaboration

The highly successful Standard Model is not complete. It does not explain the
baryonic asymmetry in the Universe, the existence of dark matter or the non-
zero masses of the neutrinos. Extensions of the Standard Model that propose
the existence of additional Heavy Neutral Leptons (HNLs) are well-motivated
and can explain several of these phenomena. In addition, light sterile
neutrinos of $\mathcal{O}(\text{eV}/c^{2})$ can explain experimentally
observed oscillation anomalies. The Neutrino Minimal Standard Model proposes
HNLs with masses $\mathcal{O}(\text{keV}/c^{2}-\text{GeV}/c^{2})$, while more
exotic models predict very large masses, up to the GUT scale. Due to the
multitude of models which hypothesize HNLs, the mass range to be explored by
experiments is large. Experimental searches for HNLs can be conducted at
existing neutrino, beam dump and collider-based experiments, and, depending on
the signature, can constrain mixing between additional neutrinos and any of
the three active neutrinos.

##  I Motivations

Heavy Neutral Leptons (HNLs) are predicted by many extensions to the Standard
Model (SM). They possess mass and therefore interact via gravity, but have no
electric charge, no weak hyper-charge, no weak isospin, and no color charge.
They are singlets under all gauge interactions and are subsequently referred
to as ‘‘sterile neutrinos." Sterile neutrinos have long been used to explain
the apparent smallness of the SM neutrino masses [1]. HNLs interact only with
the active neutrinos via mixing, and possibly the Higgs boson. Current
experimental limits are outlined in Refs. [2, 3]. The following sub-sections
document how HNLs can help resolve several known issues with the SM and how
they are ubiquitous in models that explain neutrino mass.

###  I.1 Inconsistencies between Standard Model and Observation

Standard Model predictions have been extensively tested and are found to be in
good agreement with experimental results; in some cases with extraordinary
precision of up to 1 part in a trillion [4]. Nevertheless, there remains a
need to extend the SM to explain several observational phenomena including:
the baryon asymmetry in the Universe (BAU), the existence of dark matter, and
the non-zero mass of the neutrinos.

The Neutrino Minimal Standard Model ($\nu$-MSM) [5] is one extension which
requires the existence of three HNLs. It is capable of explaining the origins
of neutrino masses, dark matter [6] and the BAU [7, 8]. The HNLs have Majorana
masses below the electroweak scale and realise the see-saw mechanism with SM
neutrinos. Two of the HNLs have masses in the $\mathcal{O}$(MeV/$c^{2}$ \-
GeV/$c^{2})$ range and a third, the dark matter candidate, has mass of
$\mathcal{O}$(keV/$c^{2}$). $\nu$-MSM is compatible with all current
measurements. Models with GeV/$c^{2}$ scale HNLs are the subject of intensive
theoretical study. Heavy Neutral Leptons with a masses from $\mathcal{O}(100$
MeV/$c^{2}$) up to a few $\mathcal{O}($GeV/$c^{2}$) can be produced in decays
of SM particles, while heavier states, of a few $\mathcal{O}($GeV/$c^{2}$),
can be directly produced at colliders.

###  I.2 Incorporating Neutrino Mass into the Standard Model

Observation of neutrino oscillations has established the non-zero mass of at
least two of the SM neutrinos. Since the discovery of neutrino oscillations a
global experimental program has measured, with precision, most of the
oscillation parameters using neutrinos from solar, accelerator, reactor and
atmospheric origins. Absolute values of these masses are yet to be determined,
but experiments have measured the mass squared differences, with current
bounds detailed in Ref.[9]. It appears that mixing in the lepton-sector has a
very different structure compared to that of the quarks. There is no current
explanation to why that might be. In addition, there are still several
remaining questions regarding neutrinos such as: establishing whether
neutrinos are Dirac or Majorana particles, their absolute masses and whether
they exhibit CP violation. Accurately answering these questions, and
accounting for neutrino masses in a consistent framework, is the focus of
global experimental and theoretical effort.

The most intuitive way to extend the SM to include the observed neutrino
masses would be to add a term coupling neutrinos to the Higgs field, analogous
to that for charged leptons:

$Y^{\nu}_{\alpha\beta}\bar{L}^{\alpha}\cdot\tilde{H}\nu^{\beta}_{R}+h.c.,$ (1)

where $\tilde{H}=i\sigma_{2}H^{*}$; $\sigma_{2}$ is the second Pauli matrix;
$H$ is the Higgs field; $Y^{\nu}_{\alpha\beta}$ describes the couplings of the
neutrino flavor states to the Higgs field; $L^{\alpha}$ is a $SU(2)_{L}$
doublet of left-handed leptons; and $\alpha,\beta$ are flavor indices.
Incorporating neutrino mass in this way requires the Yukawa coupling to be
orders of magnitude smaller than the that for charged leptons and that right-
handed neutrinos have no Majorana mass, despite there being no symmetry
preventing it. In addition, the masses and mixing angles would be expected to
have a similar hierarchy as for quarks - which is not the case.

The SM admits a dimension-five operator, the Weinberg operator, that is gauge
invariant:

$\mathcal{L}_{5}=\frac{c^{[5]}}{\Lambda}L^{T}\cdot\tilde{H}^{*}C^{\dagger}\tilde{H}^{\dagger}\cdot
L+h.c.$ (2)

where $\Lambda$ is the scale at which the particles responsible for lepton
number violation become relevant degrees of freedom; $c^{[5]}$ is a flavor-
dependent Wilson coefficient; and $C$ is the charge-conjugation matrix.

The Weinberg operator leads to Majorana masses of the neutrinos after
electroweak symmetry breaking. If the neutrinos are in fact Majorana
particles, there are a multitude of models which can generate Majorana mass
terms for left-handed fermions below the electroweak symmetry breaking scale.
These models are collectively referred to as ‘‘See-Saw Models" and can also
account for the smallness of the neutrino masses without introducing an
extremely small Yukawa coupling, they are categorized in the three ‘‘types":

  1. 1.

Type I Seesaw [10, 11, 12, 13] with a singlet fermion;

  2. 2.

Type II Seesaw [13, 14, 15, 16] with heavy triplet scalars;

  3. 3.

Type III Seesaw [17, 18] with triplet fermions.

One common feature of models that explain neutrino masses is the existence of
a new HNL state. Severe constraints exist for eV-scale Seesaw [19] from cosmic
surveys and Big-Bang Nucleo-synthesis (BBN) [20]. However, more natural
solutions remain at the GeV or TeV scale. For a search at the GeV-scale the
Yukawa coupling is of $\mathcal{O}(10^{-5})$ and can be searched for at
existing experiments. At the TeV-scale direct searches become less effective,
since the Yukawa coupling is smaller ($\mathcal{O}(10^{-6})$).

###  I.3 Suggestions of Additional Sterile Neutrino States

The existence of light sterile fermions with masses
$\mathcal{O}(\text{eV}/c^{2})$ can also provide an explanation for observed
anomalies in very short baseline oscillation measurements and cosmological
data analyses [21]. Re-analysis of data from the GALLEX [22] and SAGE [23]
solar neutrino experiments has exposed an unexplained 14$\pm$5$\%$ deficit in
the number of recorded $\nu_{e}$ \- referred to as the ‘‘Gallium anomaly" . In
addition, numerous analyses of the flux of $\bar{\nu}_{e}$ from reactors have
suggested a deficit of $\bar{\nu}_{e}$ at the 98.6$\%$ C.L. [24] \- denoted as
the ‘‘reactor anti-neutrino anomaly." When combined, both anomalies disfavour
the no-oscillation hypothesis at 99.97$\%$ ($3.6\sigma$). A third anomaly -
the ‘‘accelerator anomaly" - stems from measurements at the LSND [25]
experiment which evaluated the oscillation $\nu_{e}\rightarrow\nu_{\mu}$ at a
baseline of $L$ = 30m. LSND measured an excess of neutrinos at the level of
3.8$\sigma$ which could be explained by the existence of a sterile neutrino
with a mass $\mathcal{O}(\text{eV}/c^{2})$. Further support for this excess
was presented by the MiniBooNE experiment at the 2.8$\sigma$ level [26].

##  II Searching for Heavy Neutral Lepton States

Mixing between the beyond SM (BSM) heavy neutrino mass eigenstate and the
active neutrino states can be parameterized by the extended Pontecorvo Maki
Nakagawa Sakata (PMNS) matrix. The additional elements, $U_{l,n}$, represent
the mixing strength between the active neutrino flavor state, $l$, and the BSM
$n$-th neutrino mass state:

$\begin{pmatrix}\nu_{e}\\\ \nu_{\mu}\\\ \nu_{\tau}\\\ \nu_{\text{\tiny L}}\\\
\vdots\end{pmatrix}=\begin{pmatrix}U_{e1}&U_{e2}&U_{e3}&U_{e4}&\cdots\\\
U_{\mu 1}&U_{\mu 2}&U_{\mu 3}&U_{\mu 4}&\cdots\\\ U_{\tau 1}&U_{\tau
2}&U_{\tau 3}&U_{\tau 4}&\cdots\\\ U_{\text{\tiny L}1}&U_{\text{\tiny
L}2}&U_{\text{\tiny L}3}&U_{\text{\tiny L}4}&\cdots\\\
\vdots&\vdots&\vdots&\vdots&\ddots\\\ \end{pmatrix}\begin{pmatrix}\nu_{1}\\\
\nu_{2}\\\ \nu_{3}\\\ \nu_{4}\\\ \vdots\end{pmatrix},$ (3)

where $L$ represents some hypothetical additional lepton flavor. Given that
analyses of cosmological data and measurements of $Z$ boson decays, summarized
in Ref.[9], are consistent with there being only three charged lepton flavors
it is widely assumed that any HNL, with mass less than the $Z$ mass, must be
sterile, and would have no associated charged lepton. Here the PMNS matrix is
extended for just one HNL, but others can be added in the same way. The PMNS
matrix for the anti-neutrinos is identical to that for neutrinos under CPT
symmetry.

Experimental results tend to be presented as upper limits on
$|U_{(l=e,\mu,\tau),4}|^{2}$ (usually at the 95 $\%$ confidence level) as a
function of $m_{N}$, the mass of the possible HNL being sought. The
probability of a fourth neutrino state interacting with the electron
($|U_{e4}|^{2}$) or muon ($|U_{\mu 4}|^{2}$) has tight constraints [9], limits
on $|U_{\tau 4}|^{2}$ are weaker, motivating the possibility that $|U_{\tau
4}|>>|U_{e4}|$, $|U_{\mu 4}|$.

This article presents a summary of current bounds on the couplings to each of
the active neutrino states. In Sec.III long-established limits on all three
mixing strengths are discussed. In Sec. IV recent results, which have been
published in the last few years, are presented. In Sec. V the most recent
result from BABAR , which places new limits on $|U_{\tau 4}|^{2}$ in the 100
MeV/$c^{2}$ $<m_{4}<$ 1300 MeV/$c^{2}$ mass range. Existing bounds in the
range $\sim$ 300 MeV/$c^{2}$ to $\sim$ 1 GeV/$c^{2}$ range are particularly
weak. In Sec. VI future projections for near-term projects are discussed and
compared to current limits.

##  III  Experimental Searches

There are several means of searching for HNLs, depending on the mass regime
and proposed coupling to the active neutrinos. Additional neutral leptons can
be searched for in cosmology, in colliders, or in high-intensity experiments.
Sterile neutrinos can be responsible for contributions to electric and
magnetic leptonic dipole moments, and facilitate many rare transitions and
decays. Consequently, searches for additional neutrinos cross all frontiers of
particle physics.

###  III.1 Established Limits, prior to 2020

Figure 1: Sensitivity to Heavy Neutral Leptons with coupling to the (top)
electron, (middle) muon and (bottom) tau lepton. Current bounds (filled areas)
and near (next 5 years) future physics reach. All plots taken from Ref. [2].

####  III.1.1 Limits on Mixing with Electron and Muon Neutrinos

Mixing of heavy neutrinos with both $\nu_{e}$ and $\nu_{\mu}$ can be probed by
searching for bumps in the missing-mass distribution of pions and kaons
leptonic decays. Existing bounds come from a range of experiments:

  * •

PS 191 at CERN: [27] designed to search for neutrino decays in a low-energy
neutrino beam. No such events were found and limits were placed.

  * •

CHARM at CERN: [28] which conducted a search for HNLs in a prompt neutrino
beam produced by dumping 400 GeV protons in a Cu target. Visible decays with
electrons or muon in the final state:
$e^{+}e^{-}\nu_{e},\mu^{+}e^{-}\nu_{e},e^{+}\mu^{-}\nu_{\mu}$ and
$\mu^{+}\mu^{-}\nu_{\mu}$ were sought. This search provided limits on
$|U_{e4}|^{2},|U_{\mu 4}|^{2}<10^{-7}$ for masses up to 1.5 GeV/$c^{2}$. HNLs
were also sought through assuming a neutral-current neutrino interaction in
the CHARM calorimeter, and looking for neutrinos decaying to muons and
hadrons. This search was sensitive to masses of 0.5-2.8GeV/$c^{2}$ and
provided limits of $|U_{\mu 4}|^{2}<10^{-4}$.

  * •

Belle at KEK: [29] used 772 million $B\bar{B}$ pairs and their leptonic and
semi-leptonic decays mesons decays, $B\rightarrow(X)l\nu_{R}$, where $l=e,\mu$
and $X$ was the semi-leptonic case and either a charmed meson $D$($D^{*}$), a
light meson (e.g. $\pi.\rho,\eta$). This search had sensitivity to
$|U_{e4}|^{2}$ and $|U_{\mu 4}|^{2}$ in a range of masses between the kaon and
the $B$ meson mass.

  * •

DELPHI at LEP: [30] the best limits on all three mixing strengths in the mass
range above the $B$ meson mass come from the DELPHI detector at LEP. There
$3.3\times 10^{6}$ hadronic $Z^{0}$ decays were analyzed. Four separate
searches were performed for short-lived HNL production resulting in jets and
for longer living HNLs giving detectable secondary vertices or calorimeter
clusters. An upper limit for the branching ratio
$BR(Z^{0}\rightarrow\text{HNL}+\bar{\nu})$ of about $1.3\times 10^{-6}$ at
95$\%$ confdence level was found for masses between 3.5 and 50 GeV/$c^{2}$.

  * •

CMS at LHC: [31] used their three prompt charged lepton samples in any
combination of electrons and muons collected at a center-of-mass energy of 13
TeV. Their first results, corresponding to an integrated luminosity of
35.9$fb^{-1}$. The search is performed in the HNL mass range between 1 GeV and
1.2 TeV.

The allowed range of couplings is bounded from below by the Big Bang Nucleo-
synthesis (BBN) constraint [32], which ensures HNLs wouldn’t live long enough
in the early Universe to cause an an overabundance of Helium-4, and the See-
saw limit [33], below which the mixing to active neutrinos is too weak to
produce the observed oscillations.

####  III.1.2 Limits on Mixing with Muon Neutrinos

In addition to these searches, there are a few additional experiments which
have helped to constrain the ($|U_{\mu 4}|^{2}$, $m_{N}$) parameter space:

  1. 1.

NuTeV at Fermilab: [34] used $2\times 10^{18}$ 800 GeV protons interacting
with a beryllium oxide target and a proton dump, to search for HNLs in the
0.25–2.0 GeV/$c^{2}$ decaying to muonic final states ($\mu\mu\nu,\mu
e\nu,\mu\pi$ and $\mu\rho$). Upper limits were placed on $|U_{\mu 4}|^{2}$
down to $\mathcal{O}(10^{-7})$.

  2. 2.

E949 at BNL: [35] searched for HNLs from the process
$K^{+}\rightarrow\mu^{+}\mu_{R}$ using $1.7\times 10^{12}$ stopped kaons. They
set limits on $|U_{\mu 4}|^{2}$ of the level $\mathcal{O}(10^{-9}-10^{-7})$
for HNLs of masses 175-300 MeV/$c^{2}$.

####  III.1.3 Limits on Mixing with Tau Neutrinos

Constraints on the ($|U_{\tau 4}|^{2}$, $m_{N}$) parameter space are much
weaker; the existing constraints come from DELPHI and two other experiments in
the low mass region:

  1. 1.

CHARM at CERN: [36] used the neutrino flux from $2\times 10^{18}$ 400 GeV
protons on a solid copper target to place limits in the 10-290 MeV$/c^{2}$
mass range by re-interpreting the null result of a search for events produced
by the decay of neutral particles into two electrons.

  2. 2.

NOMAD at CERN: [37] used $4.1\times 10^{19}$ 450 GeV protons on target
collected at the WANF facility at CERN. The process
$D_{s}\rightarrow\tau\nu_{R}$ followed by the decay
$\nu_{R}\rightarrow\nu_{\tau}e^{+}e^{-}$ were analyzed and limits were placed
in the 10 to 190 MeV$/c^{2}$ mass range.

##  IV New Results, published 2020 - 2022

In the past few years several new experimental results for couplings with
muons and electrons have been published:

###  IV.1 MicroBooNe: Sterile neutrinos $\mathcal{O}(eV/c^{2})$

Recent results from MicroBooNe are detailed in Ref. [38]. MicroBooNE has
developed 3 distinct $\nu_{e}$ searches targeting the MiniBooNE excess:

  1. 1.

an exclusive search for two-body $\nu_{e}$ charged current quasi-elastic
(CCQE) scattering;

  2. 2.

a semi-inclusive search for pion-less $\nu_{e}$ events;

  3. 3.

an inclusive $\nu_{e}$ search containing any hadronic final state.

The combined results were consistent with nominal electron neutrino rate
expectations; no excess of electron neutrino events was observed. More results
follow.

###  IV.2 Kaon Searches: Couplings of HNLs $\mathcal{O}(MeV/c^{2})$ mixing
with muons and electrons

The NA62 experiment at CERN recently presented a search for for
$K^{+}\rightarrow\mu^{+}+\text{HNL}$, using the 2016-2018 data set in Ref.
[39]. The analysis found limits of $\mathcal{O}(10^{-8}$) of the neutrino
mixing parameter $|U_{\mu 4}|^{2}$ for HNL masses in the range 200-384
MeV/$c^{2}$, with lifetime exceeding 50 ns.

###  IV.3 LHC Searches: Couplings of HNLs $\mathcal{O}(GeV/c^{2})$ mixing with
muons and electrons

The ATLAS and CMS detectors at CERN have conducted searches for HNLs, detailed
in Refs. [40, 41].

CMS performed a search for HNLs produced with displaced vertices using final
states with three charged leptons (electrons or muons), the idea being that
HNLs could be produced through mixing with SM neutrinos. The decay length of
these particles can be large enough so that the secondary vertex of the HNL
decay can be resolved with the CMS silicon tracker. The selected final state
would consist of one lepton emerging from the primary proton-proton collision
vertex, along with two leptons forming a displaced, secondary vertex. In this
most recent analysis, data totalling 136 fb-1 were analyzed. Improved limits
on $|U_{e4}|^{2}$ and $|U_{\mu 4}|^{2}$ down to $\mathcal{O}(10^{-7})$ were
presented in the 1-20 GeV/$c^{2}$ mass range.

ATLAS analyzed leptonic decays of $W$ bosons extracted using 32.9 - 36.1 fb-1
of 13 TeV proton-proton collisions. HNLs could produced through mixing with
muon or electron neutrinos. They looked for both prompt and displaced leptonic
decay signatures, where the prompt signature requires three leptons produced
at the interaction point and the displaced signature comprises a prompt muon
from the $W$ boson decay and the requirement of a displaced di-lepton vertex.
This search placed limits on $|U_{e4}|^{2}$ and $|U_{\mu 4}|^{2}$ down to
$\mathcal{O}(10^{-5})$ for HNL masses in the range 4.5–50 GeV/$c^{2}$.

##  V  Couplings of HNLs $\mathcal{O}(MeV/c^{2}-GeV/c^{2})$ mixing with taus
at BABAR

The latest analysis from BABAR [42] presents new limits on $|U_{\tau 4}|^{2}$
in the 100 $<m_{4}<$ 1300 MeV/$c^{2}$ mass range. An overview of the BABAR
detector can be found in Ref. [43]. The data sample used in this analysis
corresponds to an integrated luminosity of 424 fb-1. The average cross-section
for $\tau$-pair production of electron-positron annihilation is
$\sigma(e^{+}e^{-}\rightarrow\tau^{+}\tau^{-})=(0.919\pm 0.003)$ nb [9]; so
the data sample corresponds to $\sim 4\times 10^{8}$ produced $\tau$-pairs,
before applying any reconstruction or selection criteria.

###  V.1 Experimental Strategy

This analysis used the approach proposed in Ref. [44], the key principle being
that if a HNL is produced in $\tau$ decay, the kinematics of the visible
particles would be modified with respect to SM $\tau$ decay with a massless
neutrino.

This search studies the 3-prong, pionic $\tau$ decay, which gives access to
the region 300$<m_{4}<$1360 MeV$/c^{2}$, where current limits are less
stringent. Denoting the three charged pions as a hadronic system $h^{-}$, the
decay can be considered a two-bodied:

$\tau^{-}\rightarrow h^{-}(E_{h},\vec{p}_{h})+\nu(E_{\nu},\vec{p}_{\nu}),$ (4)

where $\nu$ describes the outgoing neutrino state. An analogous equation could
be written for the $\tau^{+}$ channel. The allowed phase space of the
reconstructed energy, $E_{h}$, and invariant mass, $m_{h}$, of the hadronic
system would vary as functions of the mass of the neutrino. As the HNL gets
heavier the proportion of the original $\tau$-lepton’s energy going to the
visible pions diminishes. A visualization is presented in Fig. 2.

Figure 2: Energy and invariant mass of the hadronic system as fractions of
that of the incoming $\tau$ for cases: $m_{4}$ = 0, 500, 1000 MeV/$c^{2}$ are
shown.

In the center-of-mass frame the $\tau$-lepton energy is assumed to be
$\sqrt{s}/2$, when there is no initial state radiation. Since the direction of
the decaying $\tau$-lepton is not known, we cannot compute the neutrino mass
directly, but we know that $E_{h}$ must fall between two extremes that define
the kinematically allowed values:

$E_{\tau}-\sqrt{m^{2}_{4}+q_{+}^{2}}<E_{h}<E_{\tau}-\sqrt{m^{2}_{4}+q_{-}^{2}},$
(5)

where

$q_{\pm}=\frac{m_{\tau}}{2}\bigg{(}\frac{m_{h}^{2}-m_{\tau}^{2}-m_{4}^{2}}{m_{\tau}^{2}}\bigg{)}\sqrt{\frac{E_{\tau}^{2}}{m_{\tau}^{2}}-1}\pm\frac{E_{\tau}}{2}\sqrt{\big{(}1-\frac{(m_{h}+m_{4})^{2}}{m_{\tau}^{2}}\big{)}\big{(}1-\frac{(m_{h}-m_{4})^{2}}{m_{\tau}^{2}}\big{)}};$
(6)

and $3m_{\pi^{\pm}}<m_{h}<m_{\tau}-m_{4}$. As the HNL mass increases, the
allowed phase space of the visible system is reduced in the $E_{h}$,$m_{h}$
plane.

A HNL signal is sought by comparing the observed event yield density in the
($m_{h},E_{h}$) plane to a set of template 2D histogram distributions for the
background, obtained by simulating all $\tau$ known decays as well as
non-$\tau$ background events, and the potential HNL signal for different
$m_{4}$ mass values. Although the invariant mass and outgoing hadronic energy
($m_{h}$ and $E_{h}$) are correlated, more information can be extracted by
considering both variables simultaneously.

###  V.2 Signal and Background Simulations

Three potential sources background are considered:

  1. 1.

$\tau^{-}\rightarrow\pi^{-}\pi^{-}\pi^{+}\nu_{\tau}$, with an outgoing SM
neutrino;

  2. 2.

other SM $\tau$ decays that have been misidentified as the 3-prong (3 charged
pion) decay;

  3. 3.

non-$\tau$ backgrounds that have been misidentified as the 3-prong decay.

All SM background yields are estimated from Monte Carlo (MC) simulations which
are passed through the same reconstruction and digitization routines as the
data. All $\tau$-pair events are simulated using the KK2F [45] generator and
TAUOLA [46] which uses the averaged experimentally measured $\tau$ branching
rates as listed in Ref. [9].

Several non-$\tau$ backgrounds are also studied, including:

  * •

$e^{+}e^{-}\rightarrow\Upsilon(4S)\rightarrow B^{+}B^{-}$ and
$e^{+}e^{-}\rightarrow\Upsilon(4S)\rightarrow B^{0}$ $\bar{B}^{0}$ which are
simulated using EvtGen [47];

  * •

$e^{+}e^{-}\rightarrow u\bar{u},d\bar{d},s\bar{s}$ and $e^{+}e^{-}\rightarrow
c\bar{c}$ which are simulated using JETSET [48] [49];

  * •

$e^{+}e^{-}\rightarrow\mu^{+}\mu^{-}(\gamma)$ which are simulated using KK2F
[50].

A total of 26 signal samples were simulated, one for each of the HNL masses
across the range 100 MeV/$C^{2}$ $<m_{4}<$ 1300 MeV/$C^{2}$, at 100
MeV/$C^{2}$ increments. For each of these HNL masses, both a $\tau^{+}$ and
$\tau^{-}$ signal channel were simulated. Signal samples were produced within
the BABAR software environment using KK2F and TAUOLA by changing the value of
the outgoing neutrino mass in TAUOLA. The generated signal was passed through
the same digitization and reconstruction model as the SM background and data
samples.

###  V.3 Analysis Procedure

A binned likelihood approach is taken in which it is assumed that the contents
of a given bin, $i,j$, in the $(m_{h},E_{h})$ data histogram are distributed
as a Poisson distribution and may contain events emanating from any of the SM
background process, and potentially HNL signal events. The likelihood to
observe the selected candidates in all the $(m_{h},E_{h})$ bins is the product
of the Poisson probability to observe the selected events in each bin:

$\displaystyle\mathcal{L}=\prod_{\text{charge}}^{+-}\bigg{(}\prod_{\text{channel}}^{e\mu}\bigg{(}\prod^{ij}_{\text{bin}}\bigg{(}\frac{1}{n_{\text{obs},ij}!}\bigg{[}N_{\tau,\text{gen}}\cdot|U_{\tau
4}|^{2}\cdot p_{\text{\tiny HNL},ij}+N_{\tau,\text{gen}}\cdot(1-|U_{\tau
4}|^{2})\cdot
p_{\tau-\text{SM},ij}+n^{\text{reco}}_{BKG,ij}\bigg{]}^{(n_{\text{obs}})_{ij}}\times$
$exp\bigg{[}-(N_{\tau,\text{gen}}\cdot|U_{\tau 4}|^{2}\cdot
p_{HNL,ij}+N_{\tau,\text{gen}}\cdot(1-|U_{\tau 4}|^{2})\cdot p_{\tau-
SM,ij}+n^{\text{reco}}_{BKG,ij})\bigg{]}\bigg{)}_{\text{bin}}\times\prod_{k}f(\theta_{k},\tilde{\theta}_{k})\bigg{)}_{\text{channel}}\bigg{)}_{\text{charge}},$
(7)

where $n_{\text{obs}}$ is the number of observed events in the bin $ij$,
$N_{\tau,\text{gen}}$ is the number of generated $\tau$’s, $p_{HNL(\tau-
SM),ij}$ is the probability of a reconstructed event being in a given bin in
the HNL ($\tau-SM$) 2D template and $n^{\text{reco}}_{BKG,ij}$ is the expected
number of non-$\tau$ background events. The final product is a set of Gaussian
nuisance parameters. The expression involves a product over all bins, $ij$,
over the two 1-prong channels, and over both $\tau$-lepton charges ($\pm$).

A test statistic, $q$, can be defined as:

$q=-2\text{ln}\bigg{(}\frac{\mathcal{L}_{H_{0}}(|U_{\tau
4}|_{0}^{2};\hat{\hat{\theta}}_{0},\text{data})}{\mathcal{L}_{H_{1}}(|\hat{U}_{\tau
4}|^{2};\hat{\theta},\text{data})}\bigg{)}=-2\text{ln}(\Delta\mathcal{L}),$
(8)

where $\mathcal{L}$ in both the numerator and denominator describes the
maximized likelihood for two instances. The denominator is the maximized
(unconditional) likelihood giving the maximum likelihood estimator of
$|U_{\tau 4}|^{2}$ and the set of nuisance parameters ($\hat{\theta}$);
$\hat{\theta}$ is a vector of nuisance parameters that maximize the
likelihood. In the numerator the nuisance parameters are maximized for a given
value of $|U_{\tau 4}|^{2}$, i.e it is the conditional maximum-likelihood. The
ratio, $LR$, is consequently a function of $|U_{\tau 4}|^{2}$ through the
numerator. It must be noted that the numerator denotes the hypothesis for any
given value of $|U_{\tau 4}|^{2}$ (including the background only case i.e.
$|U_{\tau 4}|^{2}=0$). Reference [51] provides more details on likelihood-
based tests. The analysis aims to find the value of $|U_{\tau 4}|^{2}$ that
minimizes this quantity at the 95 $\%$ confidence level.

###  V.4 Uncertainties

Table 1 lists the relative contribution for each of the systematic
uncertainties on the normalization. These are parameterized as Gaussian
nuisance parameters.

Table 1: Systematic uncertainty contribution to the event yield (in $\%$) from each source, based on comparisons between MC simulations and data.  Uncertainty | Yield Change ($\pm$)
---|---
Luminosity | $0.44\%$
$\sigma(ee\rightarrow\tau\tau)$ | $0.31\%$
Branching Fractions (1 prong) | e: 0.22$\%$
|  $\mu$: 0.22$\%$
Branching Fractions (3 prong) | 3$\pi$: 0.57$\%$
PID Efficiency | $e:2\%$
|  $\mu$: 1$\%$
|  $\pi$: 3$\%$
Bhabha Contamination | 0.2$\%$
$q\bar{q}$ Contamination (data) | 0.1$\%$
Tracking Efficiency | negligible
Detector Modeling | negligible
Beam Energy | negligible
Tau Mass | negligible

In addition to these yield uncertainties, which effect all bins uniformly,
‘‘shape" uncertainties must also be accounted for. These come from
inefficiency in the MC modelling. For many hadronic $\tau$ decay channels the
relative uncertainties from experimental results are large.

A $\tau$-lepton decay to three charged pions is mediated by the $a_{1}(1260)$
resonance which decays through the intermediate $\rho\pi$ state. In the MC
samples used in this analysis the PDG [9] average of $m_{a_{1}}=1230\pm 40$
MeV/$c^{2}$ and a Breit-Wigner averaged width of $\Gamma_{a_{1}}=420\pm 35$
are used; Ref. [9] quotes the estimated width to be between 250 - 600. The
uncertainty associated with the $a_{1}$ resonance represents the dominant
contribution to the systematic error in our measurement. In order to
understand the effects of the uncertainty on the $a_{1}$ mass on the final
results in this analysis several additional MC simulations were built, in
which the $m_{a_{1}}$ was varied to $\pm 1\sigma$ of the experimental average
(where $\sigma=40$ MeV/$c^{2}$ ).

###  V.5 Results

Figure 3 shows the upper limit at the 95$\%$ confidence level provided by this
analysis using the described binned likelihood technique. The magenta line
represents the upper limit when all systematic uncertainties are considered.
To characterize deviations due to the uncertainty on $\Gamma_{a_{1}}$ the more
conservative PDG estimates are used. The dominant systematic uncertainty is,
by far, that due to the assumptions made within our simulation, the main
contribution being uncertainty in the intermediate resonances for the $\tau$
3-prong channel, and the dominant $\tau$ backgrounds. The relative systematic
uncertainty decreases as the mass of the hypothetical HNL increases, this is
expected since the effects of the modeling uncertainty become less apparent at
higher HNL masses.

Figure 3: Upper limits at 95$\%$ C.L. on $|U_{\tau 4}|^{2}$. The magenta line
represents the result when uncertainties are included. The magenta line is
expected to be a very conservative upper limit.

##  VI Future projections

References [2, 3] provide comprehensive reviews of projections for near and
far term projects. Figure-1 show near-term projections (unfilled lines) for
two experiments expected to publish data in the next few years. The new limit
from BABAR is competitive with both the projected limits from FASER (at 150
fb-1) and NA62 (at $10^{18}$) POT. The BABAR technique is also applicable to
data at Belle-II and elsewhere.

There are several other experiments planned for the next decade that will also
improve limits on mixing to all three active neutrinos. The SHiP experiment
[52] aims to make improvements in all three parameter spaces, with projected
limits down to $10^{-9}$ in the electron and muon sector and $10^{-7}$ in the
tau sector, in the few GeV/$c^{2}$ region. Other limits are expected from DUNE
[53], CODEX-b [54] and MATHUSLA [55].

Looking much further ahead, the FCC-ee could provide very powerful limits,
down to $\mathcal{O}(10^{-9}-10^{-12})$ [56, 57] for HNL masses of 5-80
GeV/$c^{2}$.

##  VII  Conclusions

To conclude, the existence of Heavy Neutral Leptons can provide solutions to
many issues which exist within the Standard Model. Depending on the mass and
coupling of the HNL to active neutrinos these new particles can be produced in
a range of channels. Consequently the experimental program searching for them
extends well beyond the neutrino sector. This article has documented recent
results from searches at neutrino experiments, beam-dumps and collider-based
experiments. This includes a new upper limit on $|U_{\tau 4}|^{2}$ set by
BABAR. The technique presented in Sec. V can also be applied future searches
at Belle-II. the results presented are competitive with projections for
experiments coming online in the next few years.

The next decade brings with it new high-intensity searches for HNLs coupling
to all three active neutrinos, many experiments are able to access mass-scales
which are predicted by well-motivated models. There is no denying that we sit
on the precipice of a very interesting time in the search for these elusive,
beyond SM particles.

###### Acknowledgements.

We are grateful for the extraordinary contributions of our PEP-II colleagues
in achieving the excellent luminosity and machine conditions that have made
this work possible. The success of this project also relies critically on the
expertise and dedication of the computing organizations that support BABAR.
The collaborating institutions wish to thank SLAC for its support and the kind
hospitality extended to them. This work is supported by the US Department of
Energy and National Science Foundation, the Natural Sciences and Engineer- ing
Research Council (Canada), Institute of High Energy Physics(China),the
Commissariat al’ Energie Atomique and Institut National de Physique Nucleaire
et de Physique des Particules (France), the Bundesministerium fur Bildung und
Forschung and Deutsche Forschungsge meinschaft (Germany), the Istituto
Nazionale di Fisica Nucleare (Italy), the Foundation for Fundamental Research
on Matter (The Netherlands), the Research Council of Norway, the Ministry of
Science and Technology of the Russian Federation, and the Particle Physics and
Astronomy Research Council (United Kingdom). Individuals have received support
from CONACYT (Mexico), the A. P. Sloan Foundation, the Research Corporation,
and the Alexander von Humboldt Foundation.

## References

  * Mohapatra and Senjanovic [1981a] R. N. Mohapatra and G. Senjanovic, Phys. Rev. D 23, 165 (1981a).
  * Beacham et al. [2019] J. Beacham et al., Journal of Physics G: Nuc. and Part. Phys. 47, 010501 (2019), ISSN 1361-6471, URL http://dx.doi.org/10.1088/1361-6471/ab4cd2.
  * Abdullahi et al. [2022] A. M. Abdullahi et al., in _2022 Snowmass Summer Study_ (2022), eprint 2203.08039.
  * Hanneke et al. [2008] D. Hanneke, S. Fogwell, and G. Gabrielse, Phys. Rev. Lett. 100, 120801 (2008), URL https://link.aps.org/doi/10.1103/PhysRevLett.100.120801.
  * Asaka and Shaposhnikov [2005a] T. Asaka and M. Shaposhnikov, Phys. Lett. B 620, 17–26 (2005a), ISSN 0370-2693, URL http://dx.doi.org/10.1016/j.physletb.2005.06.020.
  * Asaka et al. [2005] T. Asaka, S. Blanchet, and M. Shaposhnikov, Phys. Lett. B 631, 151–156 (2005), ISSN 0370-2693, URL http://dx.doi.org/10.1016/j.physletb.2005.09.070.
  * Asaka and Shaposhnikov [2005b] T. Asaka and M. Shaposhnikov, Phys. Lett. B 620, 17–26 (2005b), ISSN 0370-2693, URL http://dx.doi.org/10.1016/j.physletb.2005.06.020.
  * Boyarsky et al. [2009] A. Boyarsky, O. Ruchayskiy, and M. Shaposhnikov, Annual Review of Nuclear and Particle Science 59, 191–214 (2009), ISSN 1545-4134, URL http://dx.doi.org/10.1146/annurev.nucl.010909.083654.
  * Zyla et al. [2020] P. Zyla et al. (Particle Data Group), Prog. Theor. Exp. Phys. p. 083 C01 (2020).
  * Yanagida [1979] T. Yanagida, Conf. Proc. C 7902131, 95 (1979).
  * Mohapatra and Senjanovic [1980] R. N. Mohapatra and G. Senjanovic, Phys. Rev. Lett. 44, 912 (1980).
  * Gell-Mann et al. [1979] M. Gell-Mann, P. Ramond, and R. Slansky, Conf. Proc. C 790927, 315 (1979), eprint 1306.4669.
  * Schechter and Valle [1980] J. Schechter and J. W. F. Valle, Phys. Rev. D 22, 2227 (1980).
  * Abela et al. [1981] R. Abela, M. Daum, G. H. Eaton, R. Frosch, B. Jost, P. R. Kettle, and E. Steiner, Phys. Lett. B 105, 263 (1981), [Erratum: Phys.Lett.B 106, 513 (1981)].
  * Minehart et al. [1984] R. C. Minehart, K. O. H. Ziock, R. Marshall, W. A. Stephens, M. Daum, B. Jost, and P. R. Kettle, Phys. Rev. Lett. 52, 804 (1984).
  * Cheng and Li [1980] T. P. Cheng and L.-F. Li, Phys. Rev. D 22, 2860 (1980).
  * Lazarides et al. [1981] G. Lazarides, Q. Shafi, and C. Wetterich, Nucl. Phys. B 181, 287 (1981).
  * Mohapatra and Senjanovic [1981b] R. N. Mohapatra and G. Senjanovic, Phys. Rev. D 23, 165 (1981b).
  * Canetti and Shaposhnikov [2010a] L. Canetti and M. Shaposhnikov, Journal of Cosmology and Astroparticle Physics 2010, 001 (2010a), URL https://doi.org/10.1088%2F1475-7516%2F2010%2F09%2F001.
  * Ruchayskiy and Ivashko [2012a] O. Ruchayskiy and A. Ivashko, JCAP 10, 014 (2012a), eprint 1202.2841.
  * Palazzo [2013] A. Palazzo, Mod. Phys. Lett. A 28, 1330004 (2013), URL https://doi.org/10.1142%2Fs0217732313300048.
  * Kaether et al. [2010] F. Kaether, W. Hampel, G. Heusser, J. Kiko, and T. Kirsten, Phys. Lett. B 685, 47 (2010), eprint 1001.2731.
  * Abdurashitov and Gothers [2009] J. N. Abdurashitov and Gothers (SAGE Collaboration), Phys. Rev. C 80, 015807 (2009), URL https://link.aps.org/doi/10.1103/PhysRevC.80.015807.
  * Mention et al. [2011] G. Mention et al., Phys. Rev. D 83, 073006 (2011), URL https://link.aps.org/doi/10.1103/PhysRevD.83.073006.
  * Aguilar et al. [2001] A. Aguilar et al. (LSND Collaboration), Phys. Rev. D 64, 112007 (2001), URL https://link.aps.org/doi/10.1103/PhysRevD.64.112007.
  * Aguilar-Arevalo et al. [2013] A. A. Aguilar-Arevalo et al. (MiniBooNE Collaboration), Phys. Rev. Lett. 110, 161801 (2013), URL https://link.aps.org/doi/10.1103/PhysRevLett.110.161801.
  * Bernardi et al. [1988] G. Bernardi et al., Phys. Lett. B 203, 332 (1988).
  * Capone [1985] A. Capone (CHARM), in _20th Rencontres de Moriond: Electroweak Interactions_ (1985), pp. 295–302.
  * Liventsev et al. [2013] D. Liventsev et al. (Belle), Phys. Rev. D 87, 071102 (2013), [Erratum: Phys.Rev.D 95, 099903 (2017)], eprint 1301.1105.
  * Abreu et al. [1997] P. Abreu et al. (DELPHI), Z. Phys. C 74, 57 (1997), [Erratum: Z.Phys.C 75, 580 (1997)].
  * Sirunyan et al. [2018] A. M. Sirunyan et al. (CMS), Phys. Rev. Lett. 120, 221801 (2018), eprint 1802.02965.
  * Ruchayskiy and Ivashko [2012b] O. Ruchayskiy and A. Ivashko, JCAP 10, 014 (2012b), eprint 1202.2841.
  * Canetti and Shaposhnikov [2010b] L. Canetti and M. Shaposhnikov, JCAP 09, 001 (2010b), eprint 1006.0133.
  * Vaitaitis et al. [1999] A. Vaitaitis et al. (NuTeV Collaboration), Phys. Rev. Lett. 83, 4943 (1999).
  * Artamonov et al. [2015] A. V. Artamonov et al. (E949 Collaboration), Phys. Rev. D 91, 052001 (2015), URL https://link.aps.org/doi/10.1103/PhysRevD.91.052001.
  * Orloff et al. [2002] J. Orloff, A. Rozanov, and C. Santoni, Phys. Lett. B 550, 8–15 (2002), ISSN 0370-2693, URL http://dx.doi.org/10.1016/S0370-2693(02)02769-7.
  * Astier et al. [2001] P. Astier et al. (NOMAD Collaboration), Phys. Lett. B 506, 27–38 (2001), ISSN 0370-2693, URL http://dx.doi.org/10.1016/S0370-2693(01)00362-8.
  * Abratenko et al. [2021] P. Abratenko et al. (MicroBooNE) (2021), eprint 2110.14054.
  * Cortina Gil et al. [2021] E. Cortina Gil et al. (NA62), Phys. Lett. B 816, 136259 (2021), eprint 2101.12304.
  * Aad et al. [2019] G. Aad et al. (ATLAS), JHEP 10, 265 (2019), eprint 1905.09787.
  * Tumasyan et al. [2022] A. Tumasyan et al. (CMS) (2022), eprint 2201.05578.
  * [42] _Not yet published_.
  * Aubert et al. [2013] B. Aubert et al. (BaBar Collaboration), Nucl. Instrum. Meth. A 729, 615 (2013).
  * Kobach and Dobbs [2015] A. Kobach and S. Dobbs, Phys. Rev. D 91, 053006 (2015), URL https://link.aps.org/doi/10.1103/PhysRevD.91.053006.
  * Jadach et al. [2000] S. Jadach, B. F. L. Ward, and Z. Was, Comput. Phys. Commun. 130, 260 (2000), eprint hep-ph/9912214.
  * Was and Jadach [1992] Z. Was and S. Jadach, in _26th International Conference on High-energy Physics_ (1992), pp. 1777–1780.
  * Lange [2001] D. J. Lange, Nucl. Instrum. Meth. A 462, 152 (2001).
  * Sjöstrand [1986] T. Sjöstrand, Comp. Phys. Commu. 39, 347 (1986), ISSN 0010-4655, URL https://www.sciencedirect.com/science/article/pii/0010465586900962.
  * Sjöstrand and Bengtsson [1987] T. Sjöstrand and M. Bengtsson, Comp. Phys. Commu. 43, 367 (1987), ISSN 0010-4655, URL https://www.sciencedirect.com/science/article/pii/0010465587900543.
  * Ward et al. [2003] B. Ward, S. Jadach, and Z. Was, Nuc. Phys. B - Proceedings Supplements 116, 73–77 (2003), ISSN 0920-5632, URL http://dx.doi.org/10.1016/S0920-5632(03)80147-0.
  * Cowan et al. [2011] G. Cowan, K. Cranmer, E. Gross, and O. Vitells, Eur. Phys. J. C 71, 1554 (2011), [Erratum: Eur.Phys.J.C 73, 2501 (2013)], eprint 1007.1727.
  * Gorbunov et al. [2020] D. Gorbunov, I. Krasnov, Y. Kudenko, and S. Suvorov, Phys. Lett. B 810, 135817 (2020), eprint 2004.07974.
  * Carbajal and Gago [2022] S. Carbajal and A. M. Gago (2022), eprint 2202.09217.
  * Gligorov et al. [2018] V. V. Gligorov, S. Knapen, M. Papucci, and D. J. Robinson, Phys. Rev. D 97, 015023 (2018), URL https://link.aps.org/doi/10.1103/PhysRevD.97.015023.
  * Lubatti et al. [2020] H. Lubatti et al. (MATHUSLA), JINST 15, C06026 (2020), eprint 1901.04040.
  * Blondel and Janot [2022] A. Blondel and P. Janot, Eur. Phys. J. Plus 137 (2022).
  * Shen et al. [2022] Y.-F. Shen, J.-N. Ding, and Q. Qin, Eur. Phys. J. C 82, 398 (2022), eprint 2201.05831.