# Constraints on Heavy Neutral Leptons interacting with a singlet scalar

James M. Cline  McGill University, Department of Physics, 3600 University St.,
Montréal, QC H3A2T8 Canada  Guillermo Gambini  McGill University, Department
of Physics, 3600 University St., Montréal, QC H3A2T8 Canada  Instituto de
Física Gleb Wataghin, UNICAMP, Rua Sérgio Buarque de Holanda 777, 13083-859,
Campinas-SP, Brasil

###### Abstract

Heavy neutral leptons (HNLs) are an attractive minimal extension of the
Standard Model, as is a singlet scalar $s$ mixing with the Higgs boson. If
both are present, it is natural for HNLs to interact with $s$. For a light
singlet, the decay $N\to s\nu$ can dominate over weak HNL decays. We
reinterpret existing constraints on HNL mixing from the DELPHI, CHARM and
Belle experiments for 0.5-100 GeV mass HNLs, taking into account the new decay
channel. Although the constraints are typically weakened, in some cases they
can become stronger, due to observable $s\to\ell^{+}\ell^{-}$ decays in the
detectors. The method presented here could be used to recast constraints from
other (older) experiments without resorting to computationally expensive Monte
Carlo simulations. In addition, we update and correct some errors in the
analysis of the original constraints, in the absence of the singlet.

##  I Introduction

Right-handed neutrinos $N_{i}$ are the minimal extension of the standard model
(SM) needed to explain neutrino masses. Although their masses might be
expected to far exceed the weak scale, in principle there is no restriction on
how light they could be, down to the scale $m_{\nu}$ of active neutrino
masses. For $N$ masses in the range $0.1-1000\,$GeV, it is common to refer to
them as heavy neutral leptons (HNLs). They mix with light neutrino flavors
$\nu_{i}$ with mixing angles of order $U\sim\sqrt{m_{\nu}/m_{N}}$. The
phenomenological implications of HNLs are therefore enhanced if they are
relatively light Cvetic:2013eza ; Cvetic:2014nla ; Cvetic:2020lyh . They can
be produced in $e^{+}e^{-}$ collisions Dittmar:1989yg ; L3:1992xaz ;
DELPHI:1996qcc ; L3:2001xsz ; L3:2001zfe , neutrino beams Pais:1975dg ;
NuTeV:1999kej ; Asaka:2012bb ; MicroBooNE:2019izn ; Breitbach:2021gvv , in
beam dump experiments (from the decays of mesons) Gronau:1984ct ;
Gninenko:2012anz ; Bonivento:2013jag ; Drewes:2018gkc , and in the early
Universe, leading to constraints from Big Bang Nucleosynthesis Dolgov:2000jw ;
Sabti:2020yrt ; Boyarsky:2020dzc ; Bondarenko:2021cpc . HNL oscillations could
be observed in $W$, $\tau$, and $B_{c}$ rare decays Cvetic:2018elt ;
Tapia:2019coy ; Tapia:2021gne  and its nature (Dirac or Majorana) could be
inferred from rare meson decays Cvetic:2010rw ; Cvetic:2015naa  or Z boson
decays Blondel:2021mss . A sufficiently weakly coupled HNL can be a viable
cold dark matter candidate Asaka:2005an ; Asaka:2005pn ; Boyarsky:2009ix ;
Cline:2020mdt . HNLs have been constrained by searches at the Large Hadron
Collider CMS:2018iaf ; ATLAS:2019kpx .

Another simple and highly motivated extension of the SM is a singlet scalar
field $s$ that mixes with the Higgs boson through the interaction
$|H|^{2}s^{2}$, if $s$ gets a vacuum expectation value. Such singlets are
constrained by collider searches acciarri1996search ; cms2012search ;
LHCb:2018cjc  and Higgs decays Falkowski:2015iwa ; Bojarski:2015kra , and
their possible enhancement of the electroweak phase transition in the early
universe could produce observable gravitational waves Leitao:2015fmj ;
Huang:2016cjm ; Hashino:2016xoj ; Vaskonen:2016yiu ; Cline:2021iff , and
facilitate electroweak baryogenesis Anderson:1991zb ; McDonald:1993ey ;
Choi:1993cv ; Espinosa:2011ax ; Cline:2012hg . In the present work we will be
interested in relatively light singlets. Recent constraints are summarized in
Ref. Winkler:2018qyg , for example.

If both HNLs and singlets exist in nature, they can interact with each other
via the Lagrangian term Hostert:2020xku ; deGouvea:2019qre ; Sanchez-
Vega:2014rka ; Alvarez-Salazar:2019cxw

${\cal L}\ni g_{s}s\bar{N}N\,,$ (1)

which is possible both for Dirac or Majorana HNLs. This scenario was proposed
in Ref. Cline:2020mdt  to enable a species of HNLs to be dark matter, with a
thermal relic density from annihilations $N\bar{N}\to ss$, or $N\bar{N}\to
s^{*}\to f\bar{f}$, where $f$ is any SM particle that couples to the Higgs
boson. In a generic theory of HNLs that mix with light neutrinos, the mixing
would generate the operator $s\bar{N}\nu$, opening the new decay channel
$N\to\nu s$ if $m_{s}<m_{N}$. In this case, the usual limits on the $N$-$\nu$
mixing angle $U$ will be modified, relative to the usual assumption that $N$
decays only through the weak interactions. In this way, regions of parameter
space in the $m_{N}$-$U$ plane, that are normally considered to be ruled out,
could be reopened—or in some cases constraints can become stronger, as we will
show. It is the purpose of this paper to estimate how the modified constraints
vary with the singlet mass $m_{s}$, its mixing with the Higgs $\theta_{s}$,
and the coupling $g_{s}$. The possible values of $g_{s}$ are constrained in
the special scenario of Ref. Cline:2020mdt  where they determine the dark
matter relic density. In this study we take a more generic approach and
consider $g_{s}$ to be a free parameter, since it is not otherwise
constrained.

This is a challenging task, since it requires the reinterpretation of
experimental limits that must be dealt with individually for each experiment.
We therefore choose to limit our investigation to the HNL mass range $0.3-100$
GeV, where current constraints are set by three experiments: DELPHI (from LEP)
DELPHI:1996qcc , CHARM charm , and Belle liventsev 111For lighter and heavier
HNLs see, for example, Arguelles:2021dqn  and Das:2015toa ; Das:2016hof ,
respectively. Limits from ongoing and future experiments can be found here
Chun:2019nwi .. Despite the fact that ATLAS ATLAS:2019kpx  and CMS CMS:2018iaf
set competitive limits for $|U_{e}|^{2}$ and $|U_{\mu}|^{2}$, these are not
much stronger than DELPHI on our region of interest (below 10 GeV); hence we
focus on the first three experiments.

Although we have tried to be as quantitative as possible, our results should
be considered as indicative of more definitive limits that would require a
dedicated reanalysis of data for each experiment (as opposed to recasting
published limits), which is beyond the scope of this study.

As a first step, we must be able to reproduce existing constraints in the
absence of the singlet coupling. This part of the exercise revealed that
published limits from CHARM and Belle change somewhat when updated branching
fractions for HNL production are employed, or other corrections that we
describe in the main body of the paper. Thus another result from the present
work is improved limits from these experiments, even in the absence of a
scalar singlet.

In section II we will present the main results for the three experiments,
describing the essential characteristics relevant to our study of each one.
Further details are relegated to the appendices. In section III we put our
results into the perspective of independent constraints on the singlet scalar,
to show the relation of those limits to the recasted ones derived in section
II. Conclusions are given in section IV.

Figure 1: DELPHI: Upper limits on $|U|^{2}$ for heavy neutral leptons that
couple to a light singlet scalar, which mixes with the Higgs boson, in the HNL
mass range from $2\,m_{\mu}$ up to the kinematical limit $m_{Z}$ at LEP. The
light blue contour is DELPHI’s limit on $|U|^{2}$ at the $95\%$ C.L. in the
case of no singlet, which is independent on the HNL flavor (see figure 2).

##  II Recasted constraints on HNL$-$active neutrino mixings

Our goal is to arrive at a reasonable approximation to how HNL mixing
constraints are changed by $N\to s\nu$ decays, without however doing a full
reanalysis of experimental data, which in any case would not be feasible given
our limited understanding of detector responses and systematic errors, and
lack of access to the data. Instead we will theoretically compute the number
of events expected to be produced in a given experiment, as a function of
$m_{N}$ and mixing $U$, and use the existing constraints to calibrate the
detection efficiency. As a preview, we present our new limits from the DELPHI
experiment in Fig. 1, discussed in more detail below, where the original
constraint is the light blue contour and modified ones depending upon the
singlet mass $m_{s}$, its mixing with the Higgs $\theta$ and its coupling
$g_{s}$ to the HNL are shown.

To obtain these results, we first compute the number $N_{N}$ of HNLs produced
in $Z\to N\nu$ decays at LEP, as a function of $m_{N}$ and $U$, relative to
the total number of $N_{Z}$ of $Z$ bosons. A branching fraction
$f=N_{N}/N_{Z}$ is excluded for $f$ greater than some value $f_{0}$ depending
upon the detection efficiency. By varying $f_{0}$, we can produce contours in
the $m_{N}$-$U$ plane. If one of those contours matches the existing DELPHI
limit, we can adopt the corresponding $f_{0}$ value and claim to sufficiently
understand how to reproduce the original constraint, and then investigate how
it changes in the presence of the new decay channel. The effectiveness of this
strategy is illustrated in Fig. 2 (left), which shows the good agreement
between the original LEP limit (dashed) and our reconstruction (solid). (For
the other two experiments, where the agreement seems less good, we will argue
below that our reconstructions more accurately reflect what the true limits
should be.)

Figure 2: (Left) Dashed curve: DELPHI limits at the 95% confidence level on
$|U|^{2}$ as a function of the HNL mass DELPHI:1996qcc . These limits apply
equally to all neutrino flavors. Light blue curve: our result using equation
4. (Center) Dashed curve: Belle limits at 90$\%$ C.L. on $|U_{e}|^{2}$
liventsev . Light blue curve: our result using the procedure described above.
The difference in the region $\gtrsim$ 2 GeV is due to the use of updated
formulas for HNL production and decay (Bondarenko:2018ptm, ), and the
treatment of HNL momenta (see appendix). (Right) Dashed curve: CHARM limits at
90$\%$ C.L. on $|U_{e}|^{2}$ charm . Light blue curve: our result for these
limits. A similar result is obtained in Boiarska:2021yho .

A possible effect of $N\to s\nu$ is that the HNL decays invisibly, since $s$
may be too long-lived to decay within the detector. The constraint on $U$ then
gets weakened according to the branching fraction for weak HNL decays versus
the $N\to s\nu$ channel. On the other hand, if $s$ is short-lived and decays
into $e^{+}e^{-}$ or $\mu^{+}\mu^{-}$, within the detector, those final states
might mimic weak decays of $N$, leading to new constrained regions of
parameter space, or they might be rejected by the search, depending on the
experiment. In the following, we will assume that electromagnetic decays of
$s$ could have mimicked weak decays of $N$ for DELPHI and CHARM, but not for
Belle, as will be explained. Therefore in some regions of parameter space,
$N\to s\nu$ can give rise to new excluded regions for DELPHI and CHARM, while
for Belle it can only relax the existing constraints.

Figure 3: Possible decays of the $Z$ boson including HNL intermediate states.

###  II.1 DELPHI

The DELPHI detector at LEP I collected $3.3\times 10^{6}$ hadronic $Z$ decays
from 1991 to 1994 DELPHI:1996qcc . In these decays, HNLs $N$ could be produced
via $Z\rightarrow N\bar{\nu}$ and $Z\rightarrow\bar{N}\nu$ through the mixing
with light neutrinos

$\displaystyle N^{\prime}$ $\displaystyle\cong$ $\displaystyle
N+U_{\alpha}\nu_{\alpha}$ $\displaystyle\nu^{\prime}_{\alpha}$
$\displaystyle\cong$ $\displaystyle\nu_{\alpha}-U_{\alpha}N$ (2)

in transforming between the Lagrangian and mass eigenstates. The branching
ratio is given by

$\displaystyle\text{BR}(Z\rightarrow N\bar{\nu}_{\alpha})$ $\displaystyle=$
$\displaystyle\text{BR}(Z\rightarrow\nu_{\alpha}\bar{\nu}_{\alpha})|U_{\alpha}|^{2}$
$\displaystyle\times$
$\displaystyle\left(1-\frac{m_{N}^{2}}{M_{Z}^{2}}\right)^{2}\left(1+\frac{1}{2}\frac{m_{N}^{2}}{M_{Z}^{2}}\right)\,,$

where $\text{BR}(Z\rightarrow\nu_{\alpha}\bar{\nu}_{\alpha})\approx 0.063$ for
any $\alpha=e,\mu,\tau$.222Unlike the other experiments, LEP limits apply
equally to all flavors of HNLs Dittmar:1989yg . The mean decay length of the
HNLs is $L\cong{3}|U|^{-2}({\text{GeV}}/{m_{N}})^{6}\ \text{cm}.$ The DELPHI
Collaboration studied three different decay topologies: $\nu\ell\bar{\ell}$,
$\nu q\bar{q}$, and $\ell q\bar{q}^{\prime}$ where $\ell=e,\mu,\tau$,
$q=u,d,s,c,b$, and $q\bar{q}^{\prime}=u\bar{d},c\bar{s}$ plus charge conjugate
states. These are illustrated in Fig. 3. The fraction of $Z$ bosons leading to
observed HNLs decaying inside the detector (via weak interactions) is

$f_{w}=\text{BR}(Z\rightarrow
N)\left(1-e^{-D_{L}/L}\right)\varepsilon(m_{N}),$ (4)

where $\text{BR}(Z\rightarrow N)=\text{BR}(Z\rightarrow
N\bar{\nu}_{\alpha})+\text{BR}(Z\rightarrow\bar{N}\nu_{\alpha})$, and the
reconstruction efficiency $\varepsilon(m_{N})$ is taken from Fig. 4 of Ref.
DELPHI:1996qcc . $D_{L}$ is the length of the region in which decays are
observed. As described in the appendix, we infer this parameter (obtaining
$D_{L}=200\,$cm) when reconstructing the published limit for weak HNL decays.

Four searches were performed covering HNL masses from $2\,m_{\mu}$ up to the
kinematic limit $m_{Z}$: (i) the decay products of HNLs with short lifetimes
and small masses (monojets) were searched for within 12 cm of the interaction
point; (ii) hadronic systems (acoplanar and acollinear jets) were produced by
HNLs with short lifetimes and large masses (40-80 GeV); (iii) a HNL with an
intermediate lifetime (decaying at radii from 12 to 110 cm) could have
produced an isolated set of charged particle tracks originating from the same
vertex; (iv) HNLs with long lifetimes would decay in the detection region
where charged particle tracks cannot be reconstructed (110 to 300 cm), so this
search had to rely on localized clusters of energy depositions and hits in the
outermost layers of the detector. DELPHI computed the detection efficiencies
from $10^{5}$ signal events for HNL masses from 1.5 to 85 GeV and mean decay
lengths from 0 to 2000 cm.

In our study, we calculated the fraction of HNLs decaying inside the DELPHI
detector considering detector length $D_{L}=$ 200 cm and using the global
efficiency shown in Fig. 4 of Ref. DELPHI:1996qcc . We scaled our result to
match the published DELPHI constraint; see figure 2, left.

The next step is to add the effects of the gauge singlet scalar $s$ that
couples to HNLs with strength $g_{s}$ and mixes with the Higgs boson through a
small angle $\theta$,

$\mathcal{L}\ni-g_{s}s\bar{N}N-\frac{s_{\theta}\,m_{f}}{v}s\bar{f}f\,,$ (5)

where $s_{\theta}=\sin\theta$, $v=174\,$GeV is the complex Higgs vacuum
expectation value (VEV), and $f$ represents SM fermions with mass $m_{f}$. The
final states listed in Fig. 3 are such that decays of the singlet scalar into
fermions ($s\rightarrow f\bar{f}$) will only affect the signals for the first
two event candidates, i.e., $\nu\ell\bar{\ell}$ and $\nu q\bar{q}$ (see figure
3). However, the total decay width of the HNL will increase and so will the
probability for the HNL to decay inside the detector. Consequently, we replace
Eq. (4) with the fraction that includes the additional events from $s\to
f\bar{f}$ decays,

$\displaystyle f_{w+s}$ $\displaystyle=$ $\displaystyle\text{BR}(Z\rightarrow
N)\Big{[}\text{BR}_{w}\times(1-e^{-D_{L}/L_{N}})$ (6) $\displaystyle+$
$\displaystyle\text{BR}_{s}\times(1-e^{-D_{L}/(L_{N}+L_{s})})\Big{]}\,,$

where

$\displaystyle\text{BR}_{w}$ $\displaystyle=$
$\displaystyle\frac{\Gamma_{N}(\text{weak decays})}{\Gamma_{N}(\text{weak
decays})+\Gamma(N\rightarrow s\nu)}\,,$ $\displaystyle\text{BR}_{s}$
$\displaystyle=$ $\displaystyle 1-\text{BR}_{w}\,,$ (7)

and $L_{s}$, the mean decay length of the gauge singlet scalar $s$, is
calculated following Refs. Fradette:2017sdd ; Winkler:2018qyg . The
corresponding decay length $L_{N}$ of the HNL is given in the Appendix. We
note that HNLs can also be produced in $Z$ boson 3-body decays $e^{+}e^{-}\to
Z\to sZ^{*}\to sN\bar{\nu}$ acciarri1996search , but these are suppressed by
$s_{\theta}^{2}$ and kinematical factors in comparison with $Z\to N\bar{\nu}$
and $Z\to N\nu$, so this production channel is not shown in Fig. 3 and,
consequently, $\text{BR}(Z\rightarrow N)$ remains the same as in Eq. (4). With
these modifications, we obtain the upper bounds for HNL-active neutrino mixing
shown in Fig. 1 (left).

The modified limits on $|U|^{2}$ can be understood as the result of the
competition between HNL weak (3-body) and scalar (2-body) decays, through the
factors BRw and BRs, and the interplay between the altered decay length of the
HNL $L_{N}(\text{weak}+s)$ and that of $s$. For example, the bottom right plot
of Fig. 1 shows the weakened limits starting at the kinematic threshold
$m_{N}>m_{s}=8$ GeV. As the coupling $g_{s}$ between $N$ and $s$ increases,
$N\to s\nu$ decays become more prevalent. These are invisible decays at small
mixing $\theta=10^{-7}$, decreasing the number of signal events and weakening
the limit on $|U|^{2}$. On the other hand, at larger singlet-Higgs mixing
$\theta\geq 5\times 10^{-5}$, the singlet decays to $f\bar{f}$ with a short
enough decay length for the final state particles to be observed, as though
they were coming from weak decays. We assume that experimental sensitivity to
these events is similar to that for the weak decays, resulting in an
unmodified limit relative to the published result.

As we move to the left in Fig. 1, looking at the columns corresponding to
smaller values of $m_{s}$, the kinematic threshold discontinuity for $N\to
s\nu$ also moves to the left, until the first column where it is no longer
visible since $m_{s}<m_{N}$ for the range of $m_{N}$ considered. The pattern
described for the right-most column is similar, except that the decay length
$L_{s}$ is additionally increased by the small $m_{s}$ suppressing $s\to
f\bar{f}$ decays, leading to more invisible decays and generally weaker
limits.

Exceptionally, there are several regions where the constraint on $|U|^{2}$ is
strengthened, most notably in the plot where $m_{s}=350\,$MeV,
$\theta=10^{-2}$, $m_{s}<m_{N}\lesssim 4\,$GeV. It can be understood through
the increased signal from $N\to s\nu$ followed by $s\to f\bar{f}$ compared to
weak decays. This excluded region eventually merges back to the pure weak
decay limit as $m_{N}$ increases, since the weak decay rate scales as
$m_{N}^{5}$, while the two-body rate scales as $m_{N}^{3}$.

Figure 4: Heavy neutral lepton production from $B$ meson decays. Figure 5:
Belle efficiencies $\varepsilon(R)$ for different HNL production modes. From
left to right $m_{N}=$ 1,2,3 GeV, respectively. Solid curves: our fits for the
efficiencies. To interpolate between different values of masses, we use the
mass-dependent efficiency curves shown in Fig. 2 of Ref. liventsev . Figure 6:
Belle: Light blue curves show our reconstruction of the Belle limits on weak
decays of HNLs. Other colors show the relaxed constraints for different values
of the coupling $g_{s}$ due to $N\to s\nu$ decays.

###  II.2 Belle

The Belle experiment searched for direct HNL decays
($N\rightarrow\ell^{\pm}\pi^{\mp}$, $\ell=e,\mu$) at the KEKB $e^{+}e^{-}$
collider, where $B\bar{B}$ pairs created at the $\Upsilon(4S)$ resonance in
$e^{+}e^{-}$ collisions. The HNL production mechanism from $B$ decays is
illustrated in figure 4. The number of detected heavy neutral leptons
($N_{N}$) is given by liventsev

$\displaystyle N_{N}$ $\displaystyle=$ $\displaystyle
2N_{B\bar{B}}\,\text{BR}(B\rightarrow N)\,\text{BR}(N\rightarrow\ell\pi)$ (8)
$\displaystyle\times$
$\displaystyle\int\frac{m_{N}\Gamma}{p_{N}}\exp\left(-\frac{m_{N}\Gamma
R}{p_{N}}\right)\,\varepsilon(R)\,dR\,,$

where $p_{N}$ is the momentum of the HNL, $\Gamma$ is its total decay width,
and $\varepsilon(R)$ is its reconstruction efficiency as a function of the
distance $R$ from the interaction point. $\varepsilon(R)$ is shown for three
values of $m_{N}$ and the most relevant $B$ decay channels in Fig. 5.

The most favorable mass range in which to look for HNLs at Belle is
$M_{K}<m_{N}<M_{B}$. For this reason, the total branching fraction for HNL
production

$\text{BR}(B\rightarrow N)=\sum_{X}\text{BR}(B\rightarrow X\ell_{1}N),$ (9)

includes $X=\pi,\eta,\rho,\omega,\eta^{\prime},\phi,D,D^{*}$ for semileptonic
decays and

$\text{BR}(B\rightarrow N)=\text{BR}(B\rightarrow\ell_{1}N),$ (10)

for purely leptonic decays liventsev . Assuming HNLs to be Majorana,333For
Dirac HNLs, only the processes where the signal fermion $\ell_{2}$ (produced
by the HNL) has the opposite electric charge with respect to the production
fermion $\ell_{1}$ (coming from the decay of the $B$ meson) should be
considered.

$\displaystyle\text{BR}(N\rightarrow\ell_{2}\pi)$ $\displaystyle=$
$\displaystyle\text{BR}(N\rightarrow\ell_{2}^{-}\pi^{+})$ (11)
$\displaystyle+$ $\displaystyle\text{BR}(N\rightarrow\ell_{2}^{+}\pi^{-})\,,$

where
$\text{BR}(N\rightarrow\ell_{2}^{+}\pi^{-})=\text{BR}(N\rightarrow\ell_{2}^{-}\pi^{+})$.
Therefore the signal events of the form $\ell_{1}\ell_{2}\pi$ are
$e^{+}e^{+}\pi^{-},$ $e^{+}e^{-}\pi^{+},$ $e^{-}e^{+}\pi^{-},$ and
$e^{-}e^{-}\pi^{+}$. For the calculations of the HNL production and decay
products we followed Ref. Bondarenko:2018ptm ,444We used
$M_{\text{pole}}=\infty$ FlavourLatticeAveragingGroup:2019iem ;
DeVries:2020jbs  instead of $5.65$ GeV Bondarenko:2018ptm  for
$f_{0}^{B\rightarrow\pi}$ in the calculations of the $B$ meson form factors.
which is updated relative to the values used in the Belle analysis, and leads
to some differences in our determination of the standard HNL constraint
compared to the published version.

We start by considering the standard assumption of weak decays only. At large
masses, the decay length $L_{N}$ is so short that the HNLs decay close to the
interaction point, where the reconstruction efficiencies at short distances
$R\cong 0$ are negligible (see figure 5), and no constraint on $|U|^{2}$
arises. At lower masses, and for sufficiently small mixing, HNLs decay outside
of the detector, again leading to no constraint. On the other hand, if
$|U|^{2}$ is sufficiently large, the mean HNL decay length can be so small
that their decays do not meet experimental selection criteria, similarly to
the case of heavy HNLs. Therefore constraints arise only for relatively small
$m_{N}$ and an intermediate range of $|U|^{2}$, as shown in Fig. 2 (center).
The difference between our reconstructed limit (solid) and the original one
(dashed) is due to the updated branching ratios mentioned above.

Like for the DELPHI search, the $N$-$s$ interaction introduces competition
between three-body weak decays and two-body scalar decays of the HNL. But in
contrast to DELPHI, at Belle the $N\rightarrow s\nu$ decays cannot contribute
to signal events, which are taken to be pions and charged leptons in the final
state. Although $s$ could decay to $\ell^{+}\ell^{-}$ or $\pi\pi$, it can
never produce the combination $\pi\ell$ which is required by the search.
Therefore $N\to s\nu$ is invisible in the Belle analysis, and these decays can
only weaken the limit on $|U|^{2}$, independently of the size of the singlet-
Higgs mixing angle $\theta$.

Fig. 6 shows our results for three choices of $m_{s}$ and a range of values
for $g_{s}$, displaying again the kinematic threshold discontinuity whenever
$m_{N}=m_{s}$. As stated, the effect of the invisible decays is only to weaken
the bounds. A borderline case is the coupling $g_{s}=10^{-5}$ (green lines),
where at high HNL masses there is no change relative to the purely weak decay
bounds, since the weak decays dominate at large $m_{N}$: their rate scales as
$m_{N}^{5}$, while the rate for $N\to s\nu$ goes as $m_{N}^{3}$. At lower
$m_{N}$, provided that $m_{s}<m_{N}$ (left column), one observes a weakening
of the limits. On the other hand, for large enough $g_{s}$, the limits can
disappear entirely, again provided that $m_{s}<m_{N}$.

###  II.3 CHARM

Figure 7: Recasted CHARM limits on $|U_{e}|^{2}$ versus HNL mass $m_{N}$, on a
grid of scalar mixing versus mass.  Figure 8: Like Fig. 8, but constraining
$|U_{\mu}|^{2}$ from the CHARM experiment.  Figure 9: Recasted CHARM limits on
$|U_{\tau}|^{2}$.

Heavy neutral leptons can be produced in the semileptonic decays of $D^{\pm}$
and $D^{0}$, and in the leptonic decays of $D^{\pm}$. The CHARM Collaboration
searched for HNLs in the mass range $0.5-1.8$ GeV charm . In subsequent
analyses Gronau:1984ct ; Boiarska:2021yho , this range was expanded up to 2
GeV by considering also the production of $D_{s}$, which we emulate here.

In the CHARM search, 400 GeV protons were stopped by a copper beam dump,
producing $D$ mesons. These can decay into HNLs via mixing, whose subsequent
decay produces one or two separate electromagnetic showers ($N\rightarrow
e^{+}e^{-}\nu_{e}$), two tracks ($N\rightarrow\mu^{+}\mu^{-}\nu_{\mu}$), or
one track and one electromagnetic shower ($N\rightarrow e^{+}\mu^{-}\nu_{e}$
or $N\rightarrow\mu^{+}e^{-}\nu_{\mu}$), in an empty decay region of length
$D_{L}=35$ m and cross-sectional area $3\times 3\,{\rm m}^{2}.$

The expected number of events is given by

$\displaystyle N$ $\displaystyle=$ $\displaystyle
N_{D}\,\text{BR}(D\rightarrow
N)\,\text{BR}(N\rightarrow\ell^{\prime}\ell\nu_{\ell})\,\mathcal{A}$ (12)
$\displaystyle\times$ $\displaystyle
e^{-d/L_{N}}\left(1-e^{-D_{L}/L_{N}}\right)\varepsilon$

where $N_{D}$ is the number of $D$ mesons produced by protons in the dump,
$\mathcal{A}$ is the acceptance factor (fraction of HNLs that enter the decay
region), $d=480$ m is the distance from the interaction point to the beginning
of the decay region, $D_{L}$ is the length of the decay region, $L_{N}$ is the
mean decay length of the HNL, and $\ell,\ell^{\prime}=e,\mu$. This formula is
similar to Eq. (8) in the case where the reconstruction efficiency
$\varepsilon$ is a constant, and the integration limits correspond to the
boundaries of the detection region.555In the original CHARM analysis charm ,
the distance from the interaction point to the beginning of the detector was
ignored (taking $d=0$) incorrectly leading to exclusion of arbitrarily large
mixing. The efficiency is $\varepsilon\sim 0.6$ for HNLs of mass $m_{N}\sim
1\,$GeV charm .

To account for the singlet scalar decay channel, we modify Eq. (12) similarly
to the recasting of DELPHI, making the replacement

$\displaystyle\left(1-e^{-D_{L}/L_{N}}\right)$ $\displaystyle\to$
$\displaystyle{\rm BR}_{w}\left(1-e^{-D_{L}/L_{N}}\right)$ $\displaystyle+$
$\displaystyle{\rm BR^{\prime}}_{s}\left(1-e^{-D_{L}/(L_{N}+L_{s})}\right)\,,$

where

${\rm BR^{\prime}}_{s}={\rm BR}_{s}\times\frac{{\rm
BR}(s\to\ell\bar{\ell})}{{\rm BR}(N\to\ell^{\prime}\ell\nu)},$ (14)

because the singlet $s$ not always decays into light lepton pairs. This is
similar to the case of Belle, where BR${}^{\prime}_{s}=0$ since the singlet
cannot decay into $\pi\ell$.

The CHARM limits are sensitive to lepton flavor, so we present respective
constraints on $U_{e}$, $U_{\mu}$, $U_{\tau}$ for the cases of HNL coupling to
a single family, in Figures 7-9. Although the original CHARM analysis did not
include $U_{\tau}$ constraints, Ref. Boiarska:2021yho  extended their results
to do so by including neutral current contributions to the HNL decays, and we
have done likewise.

Like the case of DELPHI, not only can constraints be weakened by the singlet
decay channel, but in some regions of parameter space the signal can be
enhanced by singlet decays into $f\bar{f}$, leading to new excluded regions
when singlet mixing times coupling ($\theta\,g_{s}$) is large enough. For
example in the upper left plot of Fig. 7, the singlet decays within the
detector for most values of $g_{s}$, even when $|U_{e}|^{2}=1$, allowing
exclusion of large mixing angles. At larger $m_{s}$, the singlet starts to
decay before reaching the experimental decay region, and the upper boundaries
on the excluded regions reappear. The bottom right graph has a disconnected
excluded region at large $|U_{e}|^{2}$, due to the singlet decay length
starting to fall within the detection region. The constraints on
$|U_{\mu}|^{2}$ shown in Fig. 8 are quite similar. Those on $|U_{\tau}|^{2}$
in Fig. 9 are qualitatively distinct, but display similar general features. We
used the results of Ref. Boiarska:2021yho  to calibrate the sensitivity curves
at low $|U_{\tau}|$.

##  III Relation to singlet scalar bounds

The Higgs mixing versus mass parameter space of the singlet, in which we have
displayed our recasted results for DELPHI and CHARM, is independently
constrained by a variety of experiments or astrophysical considerations. In
Fig. 10 we have shown how the regions considered in our previous results
compare with the previously constrained parameter space. There it can be seen
that the largest mixing angle $s_{\theta}=0.01$ we considered is ruled out by
beam dump experiments or LHCb, except in the case of heavy singlets,
$m_{s}\gtrsim 4\,$GeV. Moreover for lighter singlets $m_{s}\lesssim 250\,$MeV,
the region of small mixing angles $s_{\theta}\sim 10^{-7}$ that we considered
is excluded by the effects of singlet decays on supernova 1987A or by BBN. We
have nevertheless included these regions in our analysis to give a complete
picture of the qualitative trends. It can be seen that our results overlap
with a significant region of singlet parameters that is currently still open.

Just as the interactions of the singlet scalar with HNLs can alter the
constraints on the HNL-active neutrino mixing angle,666even if the singlet
does not decay into the signal that a given experiment (e.g., Belle) is
looking for they can also affect the constraints on the singlet scalar-Higgs
boson mixing. These changes in the limits come from modifications of the
production and decays of the HNL and the singlet scalar.

Figure 10: Current bounds on singlet scalar-Higgs boson mixing from LEP
acciarri1996search , CMS cms2012search , LHCb LHCb:2018cjc , CHARM charm ,
E949 BNL-E949:2009dza , BaBar PhysRevLett.107.221803 ; PhysRevD.88.071102 ;
PhysRevLett.114.171801 , KTeV KTeV , BBN Fradette:2017sdd , and SN1987a
Winkler:2018qyg . Black and blue boxes show the $s_{\theta}-m_{s}$ parameter
space used to get the results presented in Fig. 1 and Figs. 7-9, respectively.
The results for Belle (see Fig. 4) are independent from the singlet scalar-
Higgs boson mixing $\theta$. Limits adopted fromWinkler:2018qyg .

The HNL width can be increased by the new decay channel $N\to s\nu$ or by
$N\to f\bar{f}\nu$ mediated by virtual $s$ exchange, competing with the weak
HNL 3-body decays. The singlet width can be increased by $s\to N\bar{\nu}$ or
the analogous process with off-shell $N$ which decays weakly. The on-shell
decays $N\to s\nu$ and $s\to N\bar{\nu}$ open when $m_{N}>m_{s}$ and
$m_{s}>m_{N}$, respectively. Because these two regions of the full parameter
space are mutually exclusive, our previous analysis is not affected by the new
decay channels of the singlet. Off-shell contributions would not modify our
results either, since these mediate processes like
$s\to\bar{\nu}N^{*}\to\bar{\nu}\ell W$, which is kinematically forbidden for
the light scalars studied in this work.

##  IV Conclusions

In this work we have estimated the changes to heavy neutral lepton mixing
constraints due to their possible decays into a light singlet scalar and an
active neutrino, for $m_{N}$ between 0.5 and $\sim 5$ GeV. One motivation for
focusing on this mass range is the possibility that one generation of such
HNLs could be the dark matter of the Universe if they have sufficiently small
mixing Cline:2020mdt , while the other generations would be subject to the
constraints investigated here.

It is possible that the limits derived here could be adapted to other
qualitatively similar models. For example, if HNLs couple to a light vector
$Z^{\prime}$, which kinetically mixes with the standard model hypercharge, it
would give rise to similar effects as we have studied, with $g_{s}$
representing the new gauge coupling and $\theta y_{f}$ mapping onto $\epsilon
e$, where $y_{f}$ is the fermion Yukawa coupling (typically for the muon) and
$\epsilon$ the kinetic mixing parameter.

Beyond the specific limits presented here, it may be that the general method
described could be useful for recasting other experimental constraints,
especially in the case of older experiments where access to original data is
not available, or Monte Carlo simulations would be difficult to carry out.

Acknowledgments. We are grateful to Dmitri Liventsev for extensive help in
understanding and reproducing the Belle HNL limits, and to Gordan Krnjaic,
Maksym Ovchynnikov and Jonathan Rosner for very helpful correspondence. This
work was supported by NSERC (Natural Sciences and Engineering Research
Council, Canada). GG acknowledges support from CNPq grant No.141699/2016-7
(Brasil), McGill Space Institute, and McGill Graduate $\&$ Postdoctoral
Studies.

##  Appendix A Methodology

The first step in reproducing the observed limit of a given experiment is to
compute the fractional number $f$ of signal events, relative to decaying
parent particles, on a grid in the
$\log_{10}|U|^{2}-\log_{10}(m_{N}/\text{GeV})$ plane. This is done without
necessarily knowing the overall normalization for the efficiency, or setting
it to unity if the efficiency can be approximated as constant as a function of
decay distance. The normalized efficiency is inferred by plotting contours of
$f$, and choosing that value of $f$ that best reproduces the published limit.
In addition, there may be other parameters that can be tuned in order to
optimize the fit, namely the size $D_{L}$ of the decay region observed by the
experiment.

The fraction $f$ is generally determined by the product of three
probabilities,

$f={N_{E}\over
N_{P}}=\sum_{P,X,Y}\mathcal{P}_{1}(P,X)\,\mathcal{P}_{2}(Y)\,\mathcal{P}_{3}(P,X)\,,$
(15)

where $N_{E}$ is the number of observed events and $N_{P}$ is the number of
parent particles $P$ whose decays could produce HNLs. $N_{P}$ might be given
or it might be computable from, for example, the number of protons on target
and the production fractions Graverini:2133817 ; SHiP:2018xqw . The three
probability factors are specified as follows.

$\mathcal{P}_{1}=\text{BR}(P\rightarrow NX)$ is the probability of HNL
production for a given decay mode of $P$, where the HNL is accompanied by
particles $X$. For example, for CHARM in the case of pure mixings with
electron neutrinos, $X=e^{+},\,e^{+}K^{0},\,e^{+}K^{*0},\,e^{+}\pi^{0}$ with
$P=D^{+}$. $X$ could be used to trigger for event candidates (Belle) or not
(DELPHI, CHARM).

$\mathcal{P}_{2}=\text{BR}(N\rightarrow Y)$ is the probability for the HNL to
decay into the signal being searched for. For DELPHI, $\mathcal{P}_{2}=1$
since all HNL decays compete with the decays of the Z bosons. In the case of
Belle, $\mathcal{P}_{2}=\text{BR}(N\rightarrow e\pi)$ for HNLs that mix only
with $\nu_{e}$. (Recall that the signal events are then $ee\pi$ where the
second lepton and the pion have opposite electric charge.) In contrast with
DELPHI and Belle, the CHARM detector is far from the interaction point, so one
must account for the fact that not all decay products of the HNLs travel
toward the detector. The acceptance factor $\mathcal{A}$ quantifies this
effect charm . $\mathcal{A}$ generally depends on the mass of the HNL and the
geometry of the experiment. In our analysis, it is taken as a free parameter
to be fit by reproducing the sensitivity of the experiment, as further
described below.

$\mathcal{P}_{3}$ is the probability of reconstructing the HNL from its decays
when it was inside the detector. Its general form is

$\mathcal{P}_{3}=\int_{d}^{d+D_{L}}\frac{e^{-R/L_{N}}}{L_{N}}\,\varepsilon(R,P,X)\,dR,$
(16)

where $\varepsilon(R,P,X)$ is the reconstruction efficiency, which depends on
the mass of the HNL, its production mode ($P,X$), and the distance $R$ from
the interaction point to where it decays. The decay length $L_{N}$ depends on
the HNL momentum, due to time dilation, and this also depends on the
production modes, but it is more sensitive to the type of decay, i.e., 2-body
or 3-body, as we describe below.

In Eq. (16), $d$ is the distance from the interaction point to the beginning
of the decay region inside the detector. For DELPHI where HNLs are created
from the decays of Z bosons at rest, $d=0$. For Belle, since the background is
higher near the interaction point, the experiment is insensitive to small-$R$
events, which are rejected by selection criteria such that $\varepsilon\to 0$
at $R=0$. Since the exact behavior of $\varepsilon$ is uncertain near $R=0$,
we take $d$ to be an undetermined small cutoff to be fit by matching Belle
constraints. For CHARM the value $d=480$ m is specified in their paper.

In order to calculate the mean decay length of the HNLs,
$L_{N}=p_{N}/(m_{N}\Gamma_{N})$, we calculated the momenta in the rest frame
of the parent particles and boosted them to the laboratory frame, neglecting
departures from the axis of the parent mesons. For 3-body decays $P\rightarrow
N\ell x$ the maximum value of the momentum of the HNL is

$\displaystyle|\vec{p}_{N}^{\,(\text{max})}|$ $\displaystyle=$
$\displaystyle\frac{1}{2m_{P}}\bigg{[}\bigg{(}m_{P}^{2}-(m_{\ell}+m_{x}+m_{N})^{2}\bigg{)}$
(17) $\displaystyle\times$
$\displaystyle\bigg{(}m_{P}^{2}-(m_{\ell}+m_{x}-m_{N})^{2}\bigg{)}\bigg{]}^{1/2}\\!\\!.$

We used $|\vec{p}_{N}|=|\vec{p}_{N}^{\,(\text{max})}|/2$ before applying
Lorentz transformations to the laboratory reference frame, assuming
$|\vec{p}_{P}|\approx$ 67 GeV, for $P=B^{\pm},D^{\pm},D^{0},$ and $D_{s}$
Gronau:1984ct .

Figure 11: Contours of detectable fraction $f$ of $Z$ bosons decaying into
HNLs at DELPHI. Black curves illustrate the dependence on $f$, while red
curves show the dependence on $D_{L}$, the assumed size of the decay detection
length, Eq. (4). Light blue curve is the best fit to the published limit
(dashed curve).

Further details of this general procedure are next presented for the
experiments of interest. Fig. 11 shows contours of the fraction $f$ of signal
events for DELPHI. The best-fit value is $\log_{10}f=-6.05$, indicating that
the experiment was sensitive to one part in $\sim 10^{6}$ $Z$ bosons decaying
into HNLs. This is the correct order of magnitude, since DELPHI produced $\sim
10^{6}$ $Z$ bosons and observed one event, compared to an expected background
of $0.8$ events. The shape of the exclusion curve can be further tuned by
varying the length of the decay region $D_{L}$. Choosing $D_{L}=200$ cm
provides the optimal fit, as illustrated by the red curves in Fig. 11.

For Belle, one must consider the dependence on the reconstruction efficiencies
$\varepsilon(R)$ (see figure 5). In this case, we have considered the limits
integrating $R$ in Eq. (8) to be from $d=1.75$ cm to $d+D_{L}=$ 60 cm. The
curve that best reproduces the original constraint on $|U|^{2}$ is obtained
with $\log_{10}f=-10.1$, regardless of the HNL flavor. With this choice, the
total number of events is $N_{E}=2\,N_{BB}\times 10^{-10.1}=0.12$, where
$N_{BB}=772\times 10^{6}$ is the number of $B\bar{B}$ pairs at Belle,
indicating the limits plotted Fig. 2 (center) are in accordance with the null
results in the search for HNLs at Belle. For these calculations we have used
updated formulas for $B$ meson branching ratios and HNL decay widths from Ref.
Bondarenko:2018ptm , after reproducing Belle’s original limits based on
superseded branching ratios Gorbunov:2007ak . The main difference between the
updated and original constraints is seen in the region near $m_{N}\gtrsim$ 2
GeV because of revisions in the branching ratios for the $B\to\rho N\ell$ and
$B\to\pi N\ell$ production modes.

For CHARM, we integrate from $d=480$ m to $d+D_{L}=515\,$m Gronau:1984ct . The
excluded region in the $|U|^{2}-m_{N}$ plane is located to the left of the
light blue curve in Fig. 2 (Right). In order to interpret $f$ as the fraction
of events for a null search for HNLs in the CHARM experiment, we take the
efficiencies to be $60\%$ as reported by the CHARM collaboration for $m_{N}=$1
GeV charm , and the acceptance factor is $\mathcal{A}\cong 10^{-3}$ (see Fig.
7 of Ref. Boiarska:2021yho ). We then fix $\log_{10}f=-15.6$ to match the
original limit, giving the number of events $N_{E}=N_{D}\times f=2.24<2.3$ at
90$\%$ C.L.. The number of $D$ mesons was determined using
$N_{D_{i}}=N_{POT}\times\chi_{c\bar{c}}\times f_{c\rightarrow D_{i}}$,
$N_{POT}=2.4\times 10^{18}$, $\chi_{c\bar{c}}\approx 4\times 10^{-3}$ for a
400 GeV proton beam CERN-SHiP-NOTE-2015-009 , $f_{c\rightarrow{D^{+}}}=0.207$,
$f_{c\rightarrow{D^{0}}}=0.632$, and $f_{c\rightarrow{D_{s}}}=0.088$
Graverini:2133817 ; SHiP:2018xqw . (“POT” denotes “protons on target.”)

## References

  * (1) G. Cvetič, C. S. Kim, and J. Zamora-Saá, “CP violations in $\pi^{\pm}$ Meson Decay,” J. Phys. G 41 (2014) 075004, arXiv:1311.7554 [hep-ph].
  * (2) G. Cvetič, C. S. Kim, and J. Zamora-Saá, “CP violation in lepton number violating semihadronic decays of $K,D,D_{s},B,B_{c}$,” Phys. Rev. D 89 no. 9, (2014) 093012, arXiv:1403.2555 [hep-ph].
  * (3) G. Cvetič, C. S. Kim, S. Mendizabal, and J. Zamora-Saá, “Exploring CP-violation, via heavy neutrino oscillations, in rare B meson decays at Belle II,” Eur. Phys. J. C 80 no. 11, (2020) 1052, arXiv:2007.04115 [hep-ph].
  * (4) M. Dittmar, A. Santamaria, M. C. Gonzalez-Garcia, and J. W. F. Valle, “Production Mechanisms and Signatures of Isosinglet Neutral Heavy Leptons in $Z^{0}$ Decays,” Nucl. Phys. B 332 (1990) 1–19.
  * (5) L3 Collaboration, O. Adriani et al., “Search for isosinglet neutral heavy leptons in $Z^{0}$ decays,” Phys. Lett. B 295 (1992) 371–382.
  * (6) DELPHI Collaboration, P. Abreu et al., “Search for neutral heavy leptons produced in Z decays,” Z. Phys. C 74 (1997) 57–71. [Erratum: Z.Phys.C 75, 580 (1997)].
  * (7) L3 Collaboration, P. Achard et al., “Search for heavy neutral and charged leptons in $e^{+}e^{-}$ annihilation at LEP,” Phys. Lett. B 517 (2001) 75–85, arXiv:hep-ex/0107015.
  * (8) L3 Collaboration, P. Achard et al., “Search for heavy isosinglet neutrino in $e^{+}e^{-}$ annihilation at LEP,” Phys. Lett. B 517 (2001) 67–74, arXiv:hep-ex/0107014.
  * (9) A. Pais and S. B. Treiman, “Neutral Heavy Leptons as a Source for Dimuon Events: A Criterion,” Phys. Rev. Lett. 35 (1975) 1206.
  * (10) NuTeV, E815 Collaboration, A. Vaitaitis et al., “Search for neutral heavy leptons in a high-energy neutrino beam,” Phys. Rev. Lett. 83 (1999) 4943–4946, arXiv:hep-ex/9908011.
  * (11) T. Asaka, S. Eijima, and A. Watanabe, “Heavy neutrino search in accelerator-based experiments,” JHEP 03 (2013) 125, arXiv:1212.1062 [hep-ph].
  * (12) MicroBooNE Collaboration, P. Abratenko et al., “Search for Heavy Neutral Leptons Decaying into Muon-Pion Pairs in the MicroBooNE Detector,” Phys. Rev. D 101 no. 5, (2020) 052001, arXiv:1911.10545 [hep-ex].
  * (13) M. Breitbach, L. Buonocore, C. Frugiuele, J. Kopp, and L. Mittnacht, “Searching for physics beyond the Standard Model in an off-axis DUNE near detector,” JHEP 01 (2022) 048, arXiv:2102.03383 [hep-ph].
  * (14) M. Gronau, C. N. Leung, and J. L. Rosner, “Extending Limits on Neutral Heavy Leptons,” Phys. Rev. D 29 (1984) 2539.
  * (15) S. N. Gninenko, D. S. Gorbunov, and M. E. Shaposhnikov, “Search for GeV-scale sterile neutrinos responsible for active neutrino oscillations and baryon asymmetry of the Universe,” Adv. High Energy Phys. 2012 (2012) 718259, arXiv:1301.5516 [hep-ph].
  * (16) W. Bonivento et al., “Proposal to Search for Heavy Neutral Leptons at the SPS,” arXiv:1310.1762 [hep-ex].
  * (17) M. Drewes, J. Hajer, J. Klaric, and G. Lanfranchi, “NA62 sensitivity to heavy neutral leptons in the low scale seesaw model,” JHEP 07 (2018) 105, arXiv:1801.04207 [hep-ph].
  * (18) 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.
  * (19) N. Sabti, A. Magalich, and A. Filimonova, “An Extended Analysis of Heavy Neutral Leptons during Big Bang Nucleosynthesis,” JCAP 11 (2020) 056, arXiv:2006.07387 [hep-ph].
  * (20) A. Boyarsky, M. Ovchynnikov, O. Ruchayskiy, and V. Syvolap, “Improved big bang nucleosynthesis constraints on heavy neutral leptons,” Phys. Rev. D 104 no. 2, (2021) 023517, arXiv:2008.00749 [hep-ph].
  * (21) K. Bondarenko, A. Boyarsky, J. Klaric, O. Mikulenko, O. Ruchayskiy, V. Syvolap, and I. Timiryasov, “An allowed window for heavy neutral leptons below the kaon mass,” JHEP 07 (2021) 193, arXiv:2101.09255 [hep-ph].
  * (22) G. Cvetič, A. Das, and J. Zamora-Saá, “Probing heavy neutrino oscillations in rare $W$ boson decays,” J. Phys. G 46 (2019) 075002, arXiv:1805.00070 [hep-ph].
  * (23) S. Tapia and J. Zamora-Saá, “Exploring CP-Violating heavy neutrino oscillations in rare tau decays at Belle II,” Nucl. Phys. B 952 (2020) 114936, arXiv:1906.09470 [hep-ph].
  * (24) S. Tapia, M. Vidal-Bravo, and J. Zamora-Saá, “Discovering heavy neutrino oscillations in rare Bc$\pm$ meson decays at HL-LHCb,” Phys. Rev. D 105 no. 3, (2022) 035003, arXiv:2109.06027 [hep-ph].
  * (25) G. Cvetič, C. Dib, S. K. Kang, and C. S. Kim, “Probing Majorana neutrinos in rare $K$ and $D,D_{s},B,B_{c}$ meson decays,” Phys. Rev. D 82 (2010) 053010, arXiv:1005.4282 [hep-ph].
  * (26) G. Cvetič, C. Dib, C. S. Kim, and J. Zamora-Saá, “Probing the Majorana neutrinos and their CP violation in decays of charged scalar mesons $\pi,K,D,D_{s},B,B_{c}$,” Symmetry 7 (2015) 726–773, arXiv:1503.01358 [hep-ph].
  * (27) A. Blondel, A. de Gouvêa, and B. Kayser, “Z-boson decays into Majorana or Dirac heavy neutrinos,” Phys. Rev. D 104 no. 5, (2021) 055027, arXiv:2105.06576 [hep-ph].
  * (28) T. Asaka, S. Blanchet, and M. Shaposhnikov, “The nuMSM, dark matter and neutrino masses,” Phys. Lett. B 631 (2005) 151–156, arXiv:hep-ph/0503065.
  * (29) T. Asaka and M. Shaposhnikov, “The $\nu$MSM, dark matter and baryon asymmetry of the universe,” Phys. Lett. B 620 (2005) 17–26, arXiv:hep-ph/0505013.
  * (30) 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].
  * (31) J. M. Cline, M. Puel, and T. Toma, “A little theory of everything, with heavy neutral leptons,” JHEP 05 (2020) 039, arXiv:2001.11505 [hep-ph].
  * (32) CMS Collaboration, A. M. Sirunyan et al., “Search for heavy neutral leptons in events with three charged leptons in proton-proton collisions at $\sqrt{s}=$ 13 TeV,” Phys. Rev. Lett. 120 no. 22, (2018) 221801, arXiv:1802.02965 [hep-ex].
  * (33) ATLAS Collaboration, G. Aad et al., “Search for heavy neutral leptons in decays of $W$ bosons produced in 13 TeV $pp$ collisions using prompt and displaced signatures with the ATLAS detector,” JHEP 10 (2019) 265, arXiv:1905.09787 [hep-ex].
  * (34) L3 Collaboration, M. Acciarri et al., “Search for neutral Higgs boson production through the process $e^{+}e^{-}\rightarrow Z^{*}H^{0}$,” Physics Letters B 385 no. 1-4, (1996) 454–470.
  * (35) CMS Collaboration, S. Chatrchyan et al., “Search for a light pseudoscalar Higgs boson in the dimuon decay channel in pp collisions at $\sqrt{s}$= 7 TeV,” Phys. Rev. Lett. 109 (2012) 121801, arXiv:1206.6326.
  * (36) LHCb Collaboration, R. Aaij et al., “Search for a dimuon resonance in the $\Upsilon$ mass region,” JHEP 09 (2018) 147, arXiv:1805.09820 [hep-ex].
  * (37) A. Falkowski, C. Gross, and O. Lebedev, “A second Higgs from the Higgs portal,” JHEP 05 (2015) 057, arXiv:1502.01361 [hep-ph].
  * (38) F. Bojarski, G. Chalons, D. Lopez-Val, and T. Robens, “Heavy to light Higgs boson decays at NLO in the Singlet Extension of the Standard Model,” JHEP 02 (2016) 147, arXiv:1511.08120 [hep-ph].
  * (39) L. Leitao and A. Megevand, “Gravitational waves from a very strong electroweak phase transition,” JCAP 05 (2016) 037, arXiv:1512.08962 [astro-ph.CO].
  * (40) P. Huang, A. J. Long, and L.-T. Wang, “Probing the Electroweak Phase Transition with Higgs Factories and Gravitational Waves,” Phys. Rev. D 94 no. 7, (2016) 075008, arXiv:1608.06619 [hep-ph].
  * (41) K. Hashino, M. Kakizaki, S. Kanemura, P. Ko, and T. Matsui, “Gravitational waves and Higgs boson couplings for exploring first order phase transition in the model with a singlet scalar field,” Phys. Lett. B 766 (2017) 49–54, arXiv:1609.00297 [hep-ph].
  * (42) V. Vaskonen, “Electroweak baryogenesis and gravitational waves from a real scalar singlet,” Phys. Rev. D 95 no. 12, (2017) 123515, arXiv:1611.02073 [hep-ph].
  * (43) J. M. Cline, A. Friedlander, D.-M. He, K. Kainulainen, B. Laurent, and D. Tucker-Smith, “Baryogenesis and gravity waves from a UV-completed electroweak phase transition,” Phys. Rev. D 103 no. 12, (2021) 123529, arXiv:2102.12490 [hep-ph].
  * (44) G. W. Anderson and L. J. Hall, “The Electroweak phase transition and baryogenesis,” Phys. Rev. D 45 (1992) 2685–2698.
  * (45) J. McDonald, “Electroweak baryogenesis and dark matter via a gauge singlet scalar,” Phys. Lett. B 323 (1994) 339–346.
  * (46) J. Choi and R. R. Volkas, “Real Higgs singlet and the electroweak phase transition in the Standard Model,” Phys. Lett. B 317 (1993) 385–391, arXiv:hep-ph/9308234.
  * (47) J. R. Espinosa, T. Konstandin, and F. Riva, “Strong Electroweak Phase Transitions in the Standard Model with a Singlet,” Nucl. Phys. B 854 (2012) 592–630, arXiv:1107.5441 [hep-ph].
  * (48) J. M. Cline and K. Kainulainen, “Electroweak baryogenesis and dark matter from a singlet Higgs,” JCAP 01 (2013) 012, arXiv:1210.4196 [hep-ph].
  * (49) M. W. Winkler, “Decay and detection of a light scalar boson mixing with the Higgs boson,” Phys. Rev. D 99 no. 1, (2019) 015018, arXiv:1809.01876 [hep-ph].
  * (50) M. Hostert and M. Pospelov, “Novel multilepton signatures of dark sectors in light meson decays,” Phys. Rev. D 105 no. 1, (2022) 015017, arXiv:2012.02142 [hep-ph].
  * (51) 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, arXiv:1911.01447 [hep-ph].
  * (52) B. L. Sánchez-Vega, J. C. Montero, and E. R. Schmitz, “Complex Scalar DM in a B-L Model,” Phys. Rev. D 90 no. 5, (2014) 055022, arXiv:1404.5973 [hep-ph].
  * (53) C. E. Alvarez-Salazar and O. L. G. Peres, “Constraining the $3-3-1$ model with heavy neutral leptons using $(g-2)_{\mu}$ and dark matter observables,” Phys. Rev. D 103 no. 3, (2021) 035029, arXiv:1906.06444 [hep-ph].
  * (54) CHARM Collaboration, J. Dorenbosch et al., “A search for decays of heavy neutrinos in the mass range 0.5 - 2.8 GeV,” Phys. Lett. B 166 (Nov, 1985) 473–478. 15 p. https://cds.cern.ch/record/164101.
  * (55) Belle Collaboration, D. Liventsev et al., “Search for heavy neutrinos at Belle,” Phys. Rev. D 87 no. 7, (2013) 071102, arXiv:1301.1105 [hep-ex]. [Erratum: Phys.Rev.D 95, 099903 (2017)].
  * (56) C. A. Argüelles, N. Foppiani, and M. Hostert, “Heavy neutral leptons below the kaon mass at hodoscopic neutrino detectors,” Phys. Rev. D 105 no. 9, (2022) 095006, arXiv:2109.03831 [hep-ph].
  * (57) A. Das and N. Okada, “Improved bounds on the heavy neutrino productions at the LHC,” Phys. Rev. D 93 no. 3, (2016) 033003, arXiv:1510.04790 [hep-ph].
  * (58) A. Das, P. Konar, and S. Majhi, “Production of Heavy neutrino in next-to-leading order QCD at the LHC and beyond,” JHEP 06 (2016) 019, arXiv:1604.00608 [hep-ph].
  * (59) E. J. Chun, A. Das, S. Mandal, M. Mitra, and N. Sinha, “Sensitivity of Lepton Number Violating Meson Decays in Different Experiments,” Phys. Rev. D 100 no. 9, (2019) 095022, arXiv:1908.09562 [hep-ph].
  * (60) K. Bondarenko, A. Boyarsky, D. Gorbunov, and O. Ruchayskiy, “Phenomenology of GeV-scale Heavy Neutral Leptons,” JHEP 11 (2018) 032, arXiv:1805.08567 [hep-ph]. arXiv:1805.08567.
  * (61) I. Boiarska, A. Boyarsky, O. Mikulenko, and M. Ovchynnikov, “Constraints from the CHARM experiment on heavy neutral leptons with tau mixing,” Phys. Rev. D 104 no. 9, (2021) 095019, arXiv:2107.14685 [hep-ph].
  * (62) A. Fradette and M. Pospelov, “BBN for the LHC: constraints on lifetimes of the Higgs portal scalars,” Phys. Rev. D 96 no. 7, (2017) 075033, arXiv:1706.01920 [hep-ph].
  * (63) Flavour Lattice Averaging Group Collaboration, S. Aoki et al., “FLAG Review 2019: Flavour Lattice Averaging Group (FLAG),” Eur. Phys. J. C 80 no. 2, (2020) 113, arXiv:1902.08191 [hep-lat].
  * (64) J. De Vries, H. K. Dreiner, J. Y. Günther, Z. S. Wang, and G. Zhou, “Long-lived Sterile Neutrinos at the LHC in Effective Field Theory,” JHEP 03 (2021) 148, arXiv:2010.07305 [hep-ph].
  * (65) BNL-E949 Collaboration, A. V. Artamonov et al., “Study of the decay $K^{+}\to\pi^{+}\nu\bar{\nu}$ in the momentum region $140<P_{\pi}<199$ MeV/c,” Phys. Rev. D 79 (2009) 092004, arXiv:0903.0030 [hep-ex].
  * (66) BaBar Collaboration, J. P. Lees et al., “Search for hadronic decays of a light higgs boson in the radiative decay $\Upsilon\rightarrow\gamma{A}^{0}$,” Phys. Rev. Lett. 107 (Nov, 2011) 221803. https://link.aps.org/doi/10.1103/PhysRevLett.107.221803.
  * (67) BaBar Collaboration, J. P. Lees et al., “Search for a low-mass scalar higgs boson decaying to a tau pair in single-photon decays of $\Upsilon(1s)$,” Phys. Rev. D 88 (Oct, 2013) 071102. https://link.aps.org/doi/10.1103/PhysRevD.88.071102.
  * (68) BaBar Collaboration, J. P. Lees et al., “Search for long-lived particles in ${e}^{+}{e}^{-}$ collisions,” Phys. Rev. Lett. 114 (Apr, 2015) 171801. https://link.aps.org/doi/10.1103/PhysRevLett.114.171801.
  * (69) A. Alavi-Harati et al., “Search for the decay ${K}_{L}\rightarrow{\mathit{\pi}}^{0}{\mathit{\mu}}^{+}{\mathit{\mu}}^{-}$,” Phys. Rev. Lett. 84 (Jun, 2000) 5279–5282. https://link.aps.org/doi/10.1103/PhysRevLett.84.5279.
  * (70) SHiP Collaboration, E. Graverini, E. Van Herwijnen, and T. Ruf, “Mass dependence of branching ratios into HNL for FairShip,”. https://cds.cern.ch/record/2133817.
  * (71) SHiP Collaboration, C. Ahdida et al., “Sensitivity of the SHiP experiment to Heavy Neutral Leptons,” JHEP 04 (2019) 077, arXiv:1811.00930 [hep-ph].
  * (72) D. Gorbunov and M. Shaposhnikov, “How to find neutral leptons of the $\nu$MSM?,” JHEP 10 (2007) 015, arXiv:0705.1729 [hep-ph]. [Erratum: JHEP 11, 101 (2013)].
  * (73) SHiP Collaboration, “Heavy Flavour Cascade Production in a Beam Dump,”. https://cds.cern.ch/record/2115534.