# Discovering the Origin of Neutrino Masses at SHiP

Seyda Ipek<EMAIL_ADDRESS>Douglas Tuckler
<EMAIL_ADDRESS>Department of Physics, Carleton University,
Ottawa, ON K1S 5B6, Canada

###### Abstract

In $U(1)_{R}$ extensions of supersymmetric models, the bino and its Dirac
partner, the singlino, can play the role of right-handed neutrinos. The bino
and the singlino form a pseudo-dirac pair, dubbed the ‘bi$\nu$o’, which can
generate Standard Model neutrino masses via the inverse seesaw mechanism. We
investigate the prospects for detecting long-lived bi$\nu$os at SHiP, where
GeV scale bi$\nu$os can be copiously produced in the decays of mesons. We show
that SHiP can probe new regions of parameter space that are complementary to
searches for the lepton flavor-violating decay $\mu\to e\gamma$. This scenario
provides a well-motivated benchmark for future experiments of a right-handed
neutrino that mixes with all Standard Model neutrinos, and is directly related
to the generation of neutrino masses.

##  I Introduction

The observation of neutrino oscillations indicates that at least two of the
Standard Model (SM) neutrinos have non-zero masses. Measurements of
atmospheric, reactor, solar, and accelerator neutrinos have determined the
neutrino mass differences and mixing angles to be Gonzalez-Garcia:2021dve

$\displaystyle\Delta m_{21}^{2}\simeq 7.4\times 10^{-5}\ \text{eV}^{2},\ \
|\Delta m_{31}^{2}|\simeq 2.5\times 10^{-3}\ \text{eV}^{2},$
$\displaystyle\sin^{2}\theta_{12}\simeq 0.3,\ \sin^{2}\theta_{23}\simeq 0.45,\
\sin^{2}\theta_{13}\simeq 0.022\ .$ (1)

In the SM, neutrinos are massless and an explanation of non-zero neutrino
masses requires beyond-the-SM (BSM) physics. A simple way to generate neutrino
masses is to introduce Majorana fermions that are SM gauge singlets called
right-handed neutrinos which lead to a suppression of SM neutrino masses via
the seesaw mechanism Minkowski:1977sc ; Gell-Mann:1979vob ; Mohapatra:1979ia ;
Yanagida:1979as . In this mechanism the light neutrino masses are inversely
proportional to the Majorana mass and hence, a very heavy scale, $M\sim
10^{16}$ GeV, is needed to explain the smallness of SM neutrino masses. In
comparison, in the Inverse Seesaw (ISS) mechanism right-handed neutrinos are
pseudo-Dirac fermions, with both Majorana and Dirac masses Mohapatra:1986aw ;
Mohapatra:1986bd ; Branco:1988ex . In this case, light neutrino masses are
_proportional_ to the (small) Majorana mass.

It has been shown that the ISS mechanism can be realized in a
$U(1)_{R}$-symmetric minimal supersymmetric SM (MSSM) Coloma:2016vod . In
$U(1)_{R}$-symmetric MSSM gauginos are necessarily pseudo-Dirac fermions. The
Dirac gaugino masses are produced via super-soft terms while the Majorana
masses, proportional to the gravitino mass $m_{3/2}$, can be produced via
anomaly mediation. The pseudo-Dirac bino can be considered as a pseudo-Dirac
right-handed neutrino, generating the neutrino masses. (The model is described
in more detail in the next section.) Such a bino is called “bi$\nu$o” to
stress the neutrino connection. Light neutrino masses in this model are
proportional to the ratio of the gravitino mass over the messenger scale
$\Lambda_{M}$. Low energy observables like BR$(\mu\to e\gamma)$ give a
constraint of $\Lambda_{M}\gtrsim 35$ TeV. For $\Lambda_{M}\sim O(100~{}{\rm
TeV})$, generating the correct neutrino mixing parameters requires
$m_{3/2}\sim O(10~{}{\rm keV})$.

The bi$\nu$o model also provides rich collider phenomenology Gehrlein:2019gtk
; Gehrlein:2021hsk . Since the bi$\nu$o mixes with neutrinos, it decays to a
mixture of leptons, quarks and missing energy. If the bi$\nu$o is heavier than
the weak scale, $M_{\tilde{B}}\gtrsim 90~{}$GeV, and the messenger scale is
not too high, $\Lambda_{M}<10^{8}$ TeV, the bi$\nu$o decays promptly. However,
if the bi$\nu$o is light, it needs to decay via off-shell $W/Z$ or Higgs to a
three-body final state, making it a long-lived particle even for
$\Lambda_{M}=100~{}$TeV. In Gehrlein:2021hsk  a long-lived bi$\nu$o signal at
MATHUSLA, FASER and CODEX-b was studied for $M_{\tilde{B}}>1$ GeV.

These earlier collider studies relied on a certain choice of the
supersymmetric mass spectrum: the bi$\nu$o is the next-to-lightest
supersymmetric particle (NLSP) and all other gauginos are heavier than the
squarks. In this scenario, the bi$\nu$o is produced via squark decays. Hence,
the constraints are cast on a combination of the squark mass, the bi$\nu$o
mass and the messenger scale. However, if the squark masses are beyond the
reach of high-energy colliders like the LHC, producing the bi$\nu$o in large
quantities could be challenging.

Because the bi$\nu$o mixes with SM neutrinos, it can be produced in any
process where a SM neutrino is produced. In particular, if they are light
enough, bi$\nu$os can be copiously produced in the decays of mesons at high
energy beam dump experiments.111Neutralino decays in R-parity violating (RPV)
MSSM can also be probed at experiments like SHiP, see deVries:2015mfw . The
model we describe here is distinctly different than generic RPV. In addition,
bi$\nu$os can become long-lived even for relatively low messenger scales. In
this paper, we investigate the prospects for GeV-scale bi$\nu$os at beam dump
experiments such as SHiP. We show that the large production rate of bi$\nu$os
from meson decays enables the probe of new regions of parameter space that are
not excluded by $\mu\to e\gamma$ and Big Bang Nucleosynthesis (BBN)
constraints. Currently, long-lived particle experiments like SHiP can set the
leading exclusion limits for bi$\nu$o masses in the $\sim 1-5$ GeV range. Our
results are summarized in Fig. 2.

This paper is organized as follows. In Sec. II we introduce the model and fix
our notation. In Sec. III we briefly describe the SHiP experiment and discuss
meson production in the proton beam dump. The production and decay
phenomenology of the bi$\nu$o is described in Sec IV. The sensitivity of SHiP
is discussed in Sec. V. We conclude in Sec. VII.

##  II Model

In this section we summarize the relevant parts of the model we study. The
details can be found in Coloma:2016vod ; Gehrlein:2019gtk .

In $U(1)_{R}$-symmetric MSSM, a global $U(1)_{R}$ is imposed on the
supersymmetric sector. The SM particles are not charged under this symmetry
while the supersymmetric partners carry $+1$ $U(1)_{R}$ charges. We work with
a modified version of this model where the global symmetry is instead
$U(1)_{R-L}$, with $L$ being the lepton number. Due to this global symmetry,
gauginos cannot be Majorana fermions. In order to give Dirac masses to
gauginos, adjoint partners with $U(1)_{R-L}$ charges of $-1$ are introduced.
For example, for the bino, $\tilde{B}$, there will be a superfield $\Phi_{S}$
whose fermionic component $S$, _singlino_ , is a SM singlet with $-1$
$U(1)_{R-L}$ charge. (Similarly, there is a tripletino and an octino as Dirac
partners to the weakinos and the gluino respectively.) We only focus on the
bino, lepton and Higgs superfields here as they are the relevant particle
content for generating the neutrino masses.

We assume supersymmetry is broken in a hidden sector, via both $F-$ and
$D-$terms and is mediated to the visible sector at a messenger scale
$\Lambda_{M}$. Furthermore, as any global symmetry, $U(1)_{R-L}$ must be
broken due to gravity. At the end, bino will have both a Dirac Fox:2002bu  and
a Majorana mass Randall:1998uk ; Giudice:1998xp :

$\displaystyle M_{\tilde{B}}=c_{i}\frac{D}{\Lambda_{M}}\,,\quad
m_{\tilde{B}}=\frac{\beta(g_{Y})}{g_{Y}}m_{3/2}\,,$ (2)

where $c_{i}$ is an $O(1)$ coefficient and
$m_{3/2}=\sum(F_{i}^{2}+D_{i}^{2}/2)/\sqrt{3}M_{\rm Pl}^{2}$ is the gravitino
mass.

In Coloma:2016vod  it was shown that the following dimension-5 and dimension-6
operators

$\frac{f_{i}}{\Lambda_{M}^{2}}\int d^{2}\theta
W^{\prime}_{\alpha}W^{\alpha}_{\tilde{B}}H_{u}L_{i}\ ,\
\frac{d_{i}}{\Lambda_{M}}\int d^{4}\theta\phi^{\dagger}\Phi_{S}H_{u}L_{i}\,,$
(3)

can generate two non-zero neutrino masses via the inverse seesaw mechanism.
(Here $\phi=1+\theta^{2}m_{3/2}$ is the conformal compensator.) Together with
the mass terms for the bino, these interactions lead to the following
Lagrangian

$\mathcal{L}\supset
M_{\tilde{B}}\tilde{B}S+M_{\tilde{B}}\tilde{B}\tilde{B}+f_{i}\frac{M_{\tilde{B}},}{\Lambda_{M}}\ell_{i}h_{u}\tilde{B}+d_{i}\frac{m_{3/2}}{\Lambda_{M}}\ell_{i}h_{u}S\,,$
(4)

where $f_{i}$ and $d_{i}$ are determined by the neutrino mass differences as

$\displaystyle f_{i}\simeq\begin{pmatrix}0.35\\\ 0.85\\\
0.35\end{pmatrix}\,,\quad d_{i}\simeq\begin{pmatrix}-0.06\\\ 0.44\\\
0.89\end{pmatrix}.$ (5)

After electroweak symmetry breaking (EWSB) the light neutrino masses are given
by

$m_{1}=0,\ m_{2}=\frac{m_{3/2}v^{2}}{\Lambda_{M}^{2}}(1-\rho),\
m_{3}=\frac{m_{3/2}v^{2}}{\Lambda_{M}^{2}}(1+\rho),$ (6)

where $\rho\simeq 0.7$ is determined by the neutrino mass splittings in Eq.
(I). Note, that the coupling of $S$ to SM particles induced after electroweak
symmetry breaking are proportional to the small gravitino mass $m_{3/2}\sim\
\mathcal{O}$(1-10 keV), and will not play a phenomenological role. Therefore,
we only focus on $\tilde{B}$ which we will call the “bi$\nu$o” for the rest of
this paper.

##  III The SHiP Experiment

The Search for Hidden Particles (SHiP) experiment is a proposed proton beam
dump experiment that uses the 400 GeV proton beam at the CERN Super Proton
Synchrotron (SPS) accelerator to search for long-lived particles
Alekhin:2015byh ; SHiP:2015vad ; SHiP:2021nfo . The proton beam extracted from
the SPS will be dumped onto a high density target and provides 2$\times
10^{20}$ protons-on-target in 5 years of operation. The entire SHiP experiment
consists of a high density target followed by a hadron absorber, a muon
shield, and a neutrino detector which all together have a length of
$\ell_{\text{sh}}=$ 64 m. Immediately following the neutrino detector is a
decay volume with a length $\ell_{\text{decay}}=50$ m, and a Hidden Sector
Decay Spectrometer (HSDS) to detect the decay products of long-lived hidden
sector particles.

Long-lived particles with MeV-GeV masses can be produced from the decays of
mesons. Because of the high energy of the proton beam and high density of the
target, a large number of mesons will be produced in the beam dump when the
incoming proton beam collides with the target. The production of kaons and
heavy flavor mesons at SHiP has been previously studied in Alekhin:2015byh ;
SHiP:2018xqw ; Gorbunov:2020rjx .

To determine the number of charged kaons produced at SHiP, we use the kaon
production fractions from Gorbunov:2020rjx . With a GEANT4 simulation that
takes into account production and propagation of kaons, it was found that
$\sim 8$ $K^{+}$s and $\sim 3.5$ $K^{-}$s are produced per proton-on-target.
About half of these are absorbed in the target, and a large fraction of the
remaining kaons decay at rest. The resulting final states are isotropic and
only a small number of the decay products will have trajectories that are in
the direction of the HSDS. Therefore, kaons that are stopped and decay at rest
are neglected. The usable kaons are those that decay in-flight and it was
found that 0.29 $K^{+}$ and 0.07 $K^{-}$ are produced per proton-on-target
which can be used for the production of BSM particles.

The number of heavy flavor mesons produced at SHiP is given by Alekhin:2015byh
; SHiP:2018xqw

$N_{M}=N_{\text{POT}}\times(2\times X_{\bar{q}q}\times
f^{q}_{\text{cascade}})\times f(q\to M)\,,$ (7)

where $N_{\text{POT}}=2\times 10^{20}$ is the number of protons-on-target,
$X_{\bar{q}q}$ is the quark-antiquark production fraction, i.e. the
probability of producing a $q\bar{q}$ pair in $pp$ collisions,
$f^{q}_{\text{cascade}}~{}(q=c,b)$ is a cascade enhancement factor that takes
into account secondary meson production in the target, and $f(q\to M)$ is the
meson production fraction – the probability that a quark $q$ will hadronize
into a meson $M$. Values for these parameters are given in Tab. 1 and detailed
discussions can be found in CERN-SHiP-NOTE-2015-009 ; SHiP:2018xqw . Note that
we only show values for $D^{\pm}_{(s)}$ and $B^{\pm}_{(c)}$ since these have
the largest branching ratios to bi$\nu$os above the kaon mass.

The total number of charged mesons produced at SHiP with
$N_{\text{POT}}=2\times 10^{20}$ are

$\displaystyle\begin{split}N_{K^{+}}=&5.8\times
10^{19},~{}~{}~{}N_{K^{-}}=1.4\times 10^{19}\,,\\\ N_{D^{\pm}}=&3.2\times
10^{17},~{}~{}~{}N_{D_{s}^{\pm}}=1.4\times 10^{17}\,,\\\
N_{B^{\pm}}=&4.5\times 10^{13},~{}~{}~{}N_{B_{c}^{\pm}}=2.8\times
10^{11}\,,\end{split}$ (8)

where we have saturated the upper bound of $f(b\to B_{c}^{\pm})$ to calculate
$N_{B_{c}^{\pm}}$. Because of the large number of mesons produced, SHiP can
have remarkable sensitivity to MeV-GeV-scale BSM particles.

$f(c\to D^{\pm})$ | 0.207 | $X_{\bar{c}c}$ | $1.7\times 10^{-3}$
---|---|---|---
$f(c\to D_{s}^{\pm})$ | 0.088 | $X_{\bar{b}b}$ | $1.6\times 10^{-7}$
$f(b\to B^{\pm})$ | 0.417 | $f^{c}_{\text{cascade}}$ | 2.3
$f(b\to B_{c}^{\pm})$ | $\leq 2.6\times 10^{-3}$ | $f^{b}_{\text{cascade}}$ | 1.7
Table 1:  Meson production fractions Graverini:2133817 , $q\bar{q}$ production
fractions HERA-B:2007rfd ; Lourenco:2006vw , and cascade enhancement factors
CERN-SHiP-NOTE-2015-009  for the heavy flavor mesons that are most relevant
for bi$\nu$o production.

Figure 1:  Left: Number of bi$\nu$os produced at SHiP assuming 5 years of
operation and $\Lambda_{M}=1$ TeV. (This messenger scale is taken as a
benchmark for comparison to other relevant models.) Right: Branching ratios of
the bi$\nu$o into different final states as a function of the bi$\nu$o mass
$M_{\tilde{B}}$.

##  IV Bi$\nu$o Phenomenology at SHiP

The bi$\nu$o mixes with the SM neutrinos via interactions given in (4). After
EWSB, this mixing will induce couplings of the bi$\nu$o to the EW gauge bosons
$W^{\pm},Z$, given by

$\mathcal{L}\supset\frac{g_{2}}{\sqrt{2}}f_{i}\frac{M_{\tilde{B}}}{\Lambda_{M}}W^{+}_{\mu}\ell_{i}\bar{\sigma}^{\mu}\tilde{B}+\frac{g_{2}}{2\cos\theta_{W}}f_{i}\frac{M_{\tilde{B}}}{\Lambda_{M}}Z_{\mu}\nu_{i}\bar{\sigma}^{\mu}\tilde{B}\,,$
(9)

where $g_{2}$ is the $SU(2)_{L}$ coupling constant and $\theta_{W}$ is the
Weinberg angle. These interactions will allow the bi$\nu$o to be produced in
any process where a SM neutrino is produced, and to decay directly to EW gauge
bosons if $M_{\tilde{B}}>M_{W,Z}$ or to SM fermions via off-shell gauge bosons
when $M_{\tilde{B}}<M_{W,Z}$. Since the mixing parameters $f_{i}$ are fixed by
the observed neutrino mass differences, unlike minimal heavy neutral lepton
(HNL) scenarios, the bi$\nu$o phenomenology is completely determined once the
bi$\nu$o mass $M_{\tilde{B}}$ and the messenger scale $\Lambda_{M}$ are fixed.

###  IV.1 Bi$\nu$o Production

In proton beam dump experiments, bi$\nu$os will be copiously produced in two-
or three-body decays of mesons in addition to secondary production from the
decays of $\tau^{\pm}$ produced in the decay of $D^{\pm}_{s}$ meson. The
number of bi$\nu$os produced is given by

$N_{\tilde{B}}=N_{i}\text{BR}(i\to\tilde{B}+X)\,,$ (10)

where $N_{i}$ is the number of mesons, $N_{M}$, or $\tau$ leptons, $N_{\tau}$,
produced and BR($i\to\tilde{B}+X$) is the branching ratio for the meson or
$\tau^{\pm}$ decay to a bi$\nu$o and final states $X$. The full expressions
for the branching ratios of meson and $\tau$ decays can be found in
Bondarenko:2018ptm ; Coloma:2020lgy .

In the left plot of Fig. 1 we show the number of bi$\nu$os produced at SHiP
assuming 5 years of operation and a messenger scale $\Lambda_{M}=1$ TeV. Note
that, even though it is excluded by low energy observables, this value of
$\Lambda_{M}$ is chosen for illustration purposes and for ease of translating
into other new physics scales. We see that the production rate of the bi$\nu$o
in association with muons (dashed curves) is typically larger than that of
electrons (solid curves) and tau leptons (dot-dashed curves) since
$f_{\mu}>f_{e,\tau}$.

###  IV.2 Bi$\nu$o Decays

Once produced, bi$\nu$os will decay via the weak interactions into two-body
final states, $\tilde{B}\to\ell^{\pm}M$, $\tilde{B}\to\nu M$, where $M$ is a
meson, and three-body final states, $\tilde{B}\to ff^{\prime}\nu_{i}$ for
final state fermions $f,f^{\prime}$. The partial widths for the various decays
are given in Gorbunov:2007ak ; Atre:2009rg ; Bondarenko:2018ptm ;
Berryman:2017twh ; Coloma:2020lgy . Note that above $M_{\tilde{B}}\sim 1$ GeV,
the decay to hadrons is more appropriately described by quark production in
the final state. Thus, for $M_{\tilde{B}}<1$ GeV we determine the total
hadronic decay rate by summing partial widths for exclusive decays to mesons.
Above 1 GeV we switch to the inclusive decay rates to quarks $\tilde{B}\to
qq^{\prime}\nu$, following the approach of Bondarenko:2018ptm .

In the right plot of Fig. 1 we show the branching ratios (BRs) as a function
of the bi$\nu$o mass $M_{\tilde{B}}$. The most promising channels to search
for are those involving charged particles in the final state. The BRs are
independent of the ratio $M_{\tilde{B}}/\Lambda_{M}$ and depend only on the
mixing parameters $f_{i}$, which are fixed by the neutrino mixing observables.
When $M_{\tilde{B}}\lesssim m_{\pi}$ the most important decay is $\tilde{B}\to
ee\nu$, while in the region $m_{\pi}\lesssim M_{\tilde{B}}\lesssim 1$ GeV the
decays $\tilde{B}\to e^{\pm}\pi^{\mp}$ and $\tilde{B}\to\mu^{\pm}\pi^{\mp}$
become dominant. The decay $\tilde{B}\to\mu^{\pm}\rho^{\mp}$ becomes important
between 1-2 GeV, while for $M_{\tilde{B}}\gtrsim 2$ GeV the decay
$\tilde{B}\to e\mu\nu$ is the most relevant.

It is important to emphasize that the bi$\nu$o is directly related to neutrino
mass generation, unlike minimal HNL scenarios that can be searched for at
SHiP.222By minimal HNLs we mean the existence of a single HNL that mixes with
only one SM neutrino flavor, which is an important benchmark for future
experimental searches Beacham:2019nyx . As a result, it is possible to gain
insights into the mechanism responsible for neutrino mass generation if a
positive signal involving both electrons and muons, for example, are observed.
Because the BRs are completely determined by the neutrino mixing parameters,
ratios of BRs could tell us if the new particle is involved in neutrino mass
generation and would be a smoking gun signal for the model described in Sec.
II. We also emphasize that, unlike RPV SUSY models, the relavant bi$\nu$o
phenomenology described here does not depend on sfermion masses. In this way,
SHiP will able to probe this model even if the sfermions are much heavier than
the LHC capabilities.

##  V SHiP Sensitivity to bi$\nu$o

Figure 2:  SHiP sensitivity to bi$\nu$os assuming 5 years of operation using
our conservative method (solid black curve) and the HNL@SHiP method (dashed
black curve). The greater sensitivity with the conservative method below
$M_{\tilde{B}}\sim 400$ MeV is due to the bi$\nu$o production from kaons,
which is not included in the HNL@SHiP package. Existing limits from $\mu\to
e\gamma$ MEG:2016leq  are depicted by the gray shaded region below the dashed
gray line labelled “$\mu\to e\gamma$”. We also show the projected exclusion
limit from the proposed Mu2e experiment Mu2e:2014fns , which will look for
$\mu\to e$ conversion in nuclei. BBN constraints are shown by the orange
shaded regions. See text for further details.

Given the SHiP experimental setup described in Sec. III, we can determine the
geometrical acceptance and efficiency for detecting the bi$\nu$o decay
products via Monte Carlo simulation. To determine the geometric acceptance, we
require that the visible final states reach the detector after the bi$\nu$o
has decayed anywhere in the decay volume. For a given position $z$ along the
decay volume, we require that the final states enter the detector, which has
$x$ and $y$ dimensions of 5 m and 10 m, respectively SHiP:2021nfo . These
conditions can be quantified by the following inequalities Giffin:2022rei :

$\displaystyle\begin{split}x_{f}=&\bigg{|}\frac{p^{\tilde{B}}_{x}}{p^{\tilde{B}}_{z}}z+\frac{p^{f}_{x}}{p^{f}_{z}}(\ell_{\text{sh}}+\ell_{\text{decay}}-z)\bigg{|}<2.5\text{m}\,,\\\
y_{f}=&\bigg{|}\frac{p^{\tilde{B}}_{y}}{p^{\tilde{B}}_{z}}z+\frac{p^{f}_{y}}{p^{f}_{z}}(\ell_{\text{sh}}+\ell_{\text{decay}}-z)\bigg{|}<5\text{m}\,,\end{split}$
(11)

where $p^{\tilde{B}}_{x,y,z}$ are the bi$\nu$o momentum components and
$p^{f}_{x,y,z}$ are the final state fermion momentum components. Signal events
are chosen to be those in which the bi$\nu$o decays into two charged particles
in the fiducial decay volume enclosed by $z_{\text{min}}=\ell_{\text{sh}}=64$
m and $z_{\text{max}}=\ell_{\text{sh}}+\ell_{\text{decay}}=114$ m, and satisfy
the conditions in (11).

The geometric acceptance depends on the location where the bi$\nu$o decays.
Therefore, the total efficiency is an integral over the length of the decay
volume given by Giffin:2022rei

$\text{eff}=M_{\tilde{B}}\Gamma\int_{z_{\text{min}}}^{z_{\text{max}}}dz\sum_{\text{events}~{}\in~{}\text{geom.}}\frac{e^{-zM_{\tilde{B}}\Gamma/p_{z}}}{N_{\text{MC}}p_{z}}\,,$
(12)

where $M_{\tilde{B}}$, $\Gamma$, and $p_{z}$ are the mass, decay width, and
$z-$component of the momentum of the bi$\nu$o, respectively. The sum is over
the events that fall within the geometric acceptance of the detector, and
$N_{\text{MC}}$ is the total number of simulated events. Efficiency plots for
minimal HNLs can be found in Batell:2020vqn ; Giffin:2022rei , for example.

Given the production, decay, and efficiency information discussed above, the
total number of signal events is given by

$N_{\text{sig}}=N_{\tilde{B}}\times\text{BR}(\tilde{B}\to
i)\times\text{eff}_{i}\,,$ (13)

where $N_{\tilde{B}}$ is the number of bi$\nu$os produced in a particular
production mode, BR($\tilde{B}\to i$) is the branching ratio for $\tilde{B}$
decaying to a final states $i$, and effi is the efficiency for detecting the
final states. Here we assume that the detector efficiency for reconstructing
the charged final states is 100%.

To determine the reach of SHiP to GeV-scale bi$\nu$os we use the median
expected exclusion significance Cowan:2010js ; Kumar:2015tna :

$Z_{\text{excl}}=\sqrt{2[s-b\ln(1+s/b)]}\,,$ (14)

where $s$ and $b$ are the number of signal and background events,
respectively. For 5 years of operation $b=0.1$ background events are expected
at SHiP Ahdida:2654870 . To set 95% C.L. exclusion limits we require
$Z_{\text{excl}}>1.645$, which corresponds to $s\approx 2$ signal events.

We use two methods to determine the SHiP sensitivity to GeV-scale long-lived
bi$\nu$os. (i) In the first method, we perform a Monte Carlo simulation to
simulate the production and decay of the bi$\nu$o. For the production, we use
Pythia8 Sjostrand:2006za ; Sjostrand:2014zea  to simulate a 400 GeV proton
beam striking a proton at rest, and extract the four momenta of mesons that
are produced. We use this output to determine the lab frame momentum of the
bi$\nu$o.

For two-body decays of the bi$\nu$o, we can analytically solve for the final
state momenta of its decay products in the rest frame of the bi$\nu$o. For
three-body decays, we use the the publicly available code muBHNL which uses
the differential decay distributions of HNLs to generate a weighted sample of
final state momenta in the rest-frame of the HNL Kelly:muBHNL ;
deGouvea:2021ual ; Kelly:2021xbv . Note that, for this purpose, a bi$\nu$o
with GeV-scale mass is a type of HNL with generic mixing with SM neutrinos,
and the kinematic output of muBHNL is directly applicable.

For both two- and three-body decays, we then use the lab-frame momentum of the
bi$\nu$o to boost the final state momenta into the lab frame. With the lab-
frame momenta of the bi$\nu$o and its decay products, we can determine the
geometric acceptance and efficiency using (12), calculate the signal rate
using (13), and find the 95% C.L. exclusion limits from (14).

For each bi$\nu$o mass, we consider only the production and decay mode which
leads to the best sensitivity, rather than a combination of all possible
production and decay modes. These exclusion limits can thus be considered to
be relatively conservative.

(ii) The second method we employ to obtain exclusion limits is by using the
Mathematica based package HNL@SHiP HNLatSHIP  that determines the SHiP
sensitivity to HNLs with arbitrary mixings with SM neutrinos; see SHiP:2018xqw
for more details. The HNL@SHiP package combines all production and decay modes
for a given bi$\nu$o mass, and the exclusion limits with this method can be
considered to be much more aggressive compared to the first method described
above.

Another important difference between the two methods is the omission of kaons
in the HNL@SHiP package. Initially, it was expected that the interaction
length of kaons in the SHiP beam dump is much shorter than their decay length
and that the kaons would be absorbed before decaying. However, it was shown in
Gorbunov:2020rjx  that the kaon production rate can be substantial. We’ve
included bi$\nu$o production from kaons in method (i), which sets the leading
SHiP constraints below $M_{\tilde{B}}\simeq 400$ MeV.

##  VI Results and Discussion

The exclusion limits derived using these two analyses are depicted in Fig. 2.
The solid curve corresponds to the more conservative analysis, where we only
consider the channels that lead to the best sensitivity rather than a
combination of all production and decay channels. The dashed black curve
depicts the exclusion limits derived using the HNL@SHiP code.

Additionally, we show current and future constraints from $\mu\to e\gamma$
MEG:2016leq  and $\mu\to e$ conversion in nuclei Mu2e:2014fns  by the
horizontal dashed gray lines labeled “$\mu\to e\gamma$” and “Mu2e Projected”,
respectively. The current bound of BR($\mu\to e\gamma)<4.2\times 10^{-13}$
results in a lower bound on the messenger scale $\Lambda_{M}\gtrsim 35$ TeV,
while future projections of the sensitivity of the Mu2e experiment to $\mu\to
e$ conversion near nuclei is expected to probe upto $\Lambda_{M}\lesssim 100$
TeV. These constraints are independent of the bi$\nu$o mass.

Finally, the shaded orange regions are constraints from BBN. To obtain these
bounds, we assume that the bi$\nu$o mixes with a single neutrino flavor
$\nu_{i}$ with a mixing parameter $U_{i}\sim f_{i}M_{\tilde{B}}/\Lambda_{M}$,
with $f_{i}$ given by (5). We then use the results of Sabti:2020yrt ;
Boyarsky:2020dzc  for BBN constraints on HNLs that mix with a single neutrino
flavor to constrain the bi$\nu$o parameter space. This is depicted by the
three orange shaded regions. The regions enclosed by the solid, dashed, and
dot-dashed orange curves are BBN constraints for electron-, muon-, and tau-
mixed HNLs, respectively. For the bi$\nu$o, which mixes with all three
neutrino flavors, BBN constraints will generally lie in between the dashed and
the dot-dashed curves enclosing the upper-left and lower-right orange shaded
region.

We observe that already with our conservative analysis (the solid black curved
in Fig. 2) SHiP can begin to probe the $M_{\tilde{B}}=1-2$ GeV window that is
not excluded by the current $\mu\to e\gamma$ and BBN limits, with messenger
scales up to $\Lambda_{M}\sim 60$ TeV. Combining all possible production and
decay modes is expected to lead to better exclusion limits allowing SHiP to
probe the mass region above the $D^{\pm}_{s}$ meson mass for the messenger
scales not excluded by current limits. We also see that SHiP can even probe
messenger scales beyond the projected sensitivity of the Mu2e experiment.

These results show that SHiP is highly complementary to experiments looking
for charged lepton flavor violation, as well as constraints from BBN for the
bi$\nu$o model. SHiP is able to probe new parameter space for bi$\nu$o masses
in the $M_{\tilde{B}}\sim 1-5$ GeV range, corresponding to a messenger scale
in the $\Lambda_{M}\sim 60-200$ TeV range.

##  VII Conclusion

In this paper, we investigated the sensitivity of the SHiP experiment to the
bino mass and the SUSY messenger scale in a $U(1)_{R-L}$ extension of the MSSM
which explains the origin of neutrino masses. (In this scenario, the bino and
its Dirac partner singlino act like right-handed neutrinos and are responsible
for generating the neutrino masses.) Previous studies have considered the
production of bi$\nu$os at the LHC from the decays of squarks, and constrained
a combination of the squark mass, the bi$\nu$o mass, and the messenger scale
$\lambda_{M}$. However, if the squarks are too heavy to be produced at the
LHC, weak interactions might be too weak to produce bi$\nu$os at the LHC.

On the other hand, if the bi$\nu$o mass is in the MeV-GeV range, a large
production rate can be achieved in high energy beam dump experiments. Mixing
with SM neutrinos allows the bi$\nu$o to be produced from meson decays, which
have large production rates at SHiP. We showed that SHiP can probe messenger
scales up to 200 TeV and probe a parameter region complementary to experiments
looking for charged-lepton flavor violation.

The model we investigated is directly related to neutrino mass generation that
can be discovered at SHiP. If final states involving both muons and electrons
are observed, their relative widths are fully predicted by neutrino-mass-
mixing measurements. This would provide an incredible opportunity to search
for the explanation of one of the most pressing problems with the SM.

## Acknowledgements

This work is supported in part by the Natural Sciences and Engineering
Research Council (NSERC) of Canada. The work of DT is supported by the Arthur
B. McDonald Canadian Astroparticle Physics Research Institute.

## References

  * (1) M. C. Gonzalez-Garcia, M. Maltoni, and T. Schwetz, “NuFIT: Three-Flavour Global Analyses of Neutrino Oscillation Experiments,” Universe 7 (2021) no. 12, 459, arXiv:2111.03086 [hep-ph].
  * (2) P. Minkowski, “$\mu\to e\gamma$ at a Rate of One Out of $10^{9}$ Muon Decays?,” Phys. Lett. B 67 (1977) 421–428.
  * (3) M. Gell-Mann, P. Ramond, and R. Slansky, “Complex Spinors and Unified Theories,” Conf. Proc. C 790927 (1979) 315–321, arXiv:1306.4669 [hep-th].
  * (4) R. N. Mohapatra and G. Senjanovic, “Neutrino Mass and Spontaneous Parity Nonconservation,” Phys. Rev. Lett. 44 (1980) 912.
  * (5) T. Yanagida, “Horizontal gauge symmetry and masses of neutrinos,” Conf. Proc. C 7902131 (1979) 95–99.
  * (6) R. N. Mohapatra, “Mechanism for Understanding Small Neutrino Mass in Superstring Theories,” Phys. Rev. Lett. 56 (1986) 561–563.
  * (7) R. N. Mohapatra and J. W. F. Valle, “Neutrino Mass and Baryon Number Nonconservation in Superstring Models,” Phys. Rev. D 34 (1986) 1642.
  * (8) G. C. Branco, W. Grimus, and L. Lavoura, “The Seesaw Mechanism in the Presence of a Conserved Lepton Number,” Nucl. Phys. B 312 (1989) 492–508.
  * (9) P. Coloma and S. Ipek, “Neutrino masses from a pseudo-Dirac Bino,” Phys. Rev. Lett. 117 (2016) no. 11, 111803, arXiv:1606.06372 [hep-ph].
  * (10) J. Gehrlein, S. Ipek, and P. J. Fox, “Bi$\boldsymbol{\nu}$o Phenomenology at the LHC,” JHEP 03 (2019) 073, arXiv:1901.09284 [hep-ph].
  * (11) J. Gehrlein and S. Ipek, “Long-lived bi$\boldsymbol{\nu}$o at the LHC,” JHEP 05 (2021) 020, arXiv:2103.01251 [hep-ph].
  * (12) J. de Vries, H. K. Dreiner, and D. Schmeier, “R-Parity Violation and Light Neutralinos at SHiP and the LHC,” Phys. Rev. D 94 (2016) no. 3, 035006, arXiv:1511.07436 [hep-ph].
  * (13) P. J. Fox, A. E. Nelson, and N. Weiner, “Dirac gaugino masses and supersoft supersymmetry breaking,” JHEP 08 (2002) 035, arXiv:hep-ph/0206096.
  * (14) L. Randall and R. Sundrum, “Out of this world supersymmetry breaking,” Nucl. Phys. B 557 (1999) 79–118, arXiv:hep-th/9810155.
  * (15) G. F. Giudice, M. A. Luty, H. Murayama, and R. Rattazzi, “Gaugino mass without singlets,” JHEP 12 (1998) 027, arXiv:hep-ph/9810442.
  * (16) S. Alekhin et al., “A facility to Search for Hidden Particles at the CERN SPS: the SHiP physics case,” Rept. Prog. Phys. 79 (2016) no. 12, 124201, arXiv:1504.04855 [hep-ph].
  * (17) SHiP Collaboration, M. Anelli et al., “A facility to Search for Hidden Particles (SHiP) at the CERN SPS,” arXiv:1504.04956 [physics.ins-det].
  * (18) SHiP Collaboration, C. Ahdida et al., “The SHiP experiment at the proposed CERN SPS Beam Dump Facility,” Eur. Phys. J. C 82 (2022) no. 5, 486, arXiv:2112.01487 [physics.ins-det].
  * (19) SHiP Collaboration, C. Ahdida et al., “Sensitivity of the SHiP experiment to Heavy Neutral Leptons,” JHEP 04 (2019) 077, arXiv:1811.00930 [hep-ph].
  * (20) D. Gorbunov, I. Krasnov, Y. Kudenko, and S. Suvorov, “Heavy Neutral Leptons from kaon decays in the SHiP experiment,” Phys. Lett. B 810 (2020) 135817, arXiv:2004.07974 [hep-ph].
  * (21) SHiP Collaboration, “Heavy Flavour Cascade Production in a Beam Dump,”. https://cds.cern.ch/record/2115534.
  * (22) 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.
  * (23) HERA-B Collaboration, I. Abt et al., “Measurement of D0, D+, D+(s) and D*+ Production in Fixed Target 920-GeV Proton-Nucleus Collisions,” Eur. Phys. J. C 52 (2007) 531–542, arXiv:0708.1443 [hep-ex].
  * (24) C. Lourenco and H. K. Wohri, “Heavy flavour hadro-production from fixed-target to collider energies,” Phys. Rept. 433 (2006) 127–180, arXiv:hep-ph/0609101.
  * (25) 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].
  * (26) P. Coloma, E. Fernández-Martínez, M. González-López, J. Hernández-García, and Z. Pavlovic, “GeV-scale neutrinos: interactions with mesons and DUNE sensitivity,” Eur. Phys. J. C 81 (2021) no. 1, 78, arXiv:2007.03701 [hep-ph].
  * (27) 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)].
  * (28) A. Atre, T. Han, S. Pascoli, and B. Zhang, “The Search for Heavy Majorana Neutrinos,” JHEP 05 (2009) 030, arXiv:0901.3589 [hep-ph].
  * (29) J. M. Berryman, A. de Gouvêa, K. J. Kelly, and Y. Zhang, “Dark Matter and Neutrino Mass from the Smallest Non-Abelian Chiral Dark Sector,” Phys. Rev. D 96 (2017) no. 7, 075010, arXiv:1706.02722 [hep-ph].
  * (30) J. Beacham et al., “Physics Beyond Colliders at CERN: Beyond the Standard Model Working Group Report,” J. Phys. G 47 (2020) no. 1, 010501, arXiv:1901.09966 [hep-ex].
  * (31) MEG Collaboration, A. M. Baldini et al., “Search for the lepton flavour violating decay $\mu^{+}\rightarrow\mathrm{e}^{+}\gamma$ with the full dataset of the MEG experiment,” Eur. Phys. J. C 76 (2016) no. 8, 434, arXiv:1605.05081 [hep-ex].
  * (32) Mu2e Collaboration, L. Bartoszek et al., “Mu2e Technical Design Report,” arXiv:1501.05241 [physics.ins-det].
  * (33) P. Giffin, S. Gori, Y.-D. Tsai, and D. Tuckler, “Heavy Neutral Leptons at Beam Dump Experiments of Future Lepton Colliders,” arXiv:2206.13745 [hep-ph].
  * (34) B. Batell, J. A. Evans, S. Gori, and M. Rai, “Dark Scalars and Heavy Neutral Leptons at DarkQuest,” JHEP 05 (2021) 049, arXiv:2008.08108 [hep-ph].
  * (35) G. Cowan, K. Cranmer, E. Gross, and O. Vitells, “Asymptotic formulae for likelihood-based tests of new physics,” Eur. Phys. J. C 71 (2011) 1554, arXiv:1007.1727 [physics.data-an]. [Erratum: Eur.Phys.J.C 73, 2501 (2013)].
  * (36) N. Kumar and S. P. Martin, “Vectorlike Leptons at the Large Hadron Collider,” Phys. Rev. D 92 (2015) no. 11, 115018, arXiv:1510.03456 [hep-ph].
  * (37) SHiP Collaboration, C. Ahdida et al., “SHiP Experiment - Progress Report,” tech. rep., CERN, Geneva, 2019.  https://cds.cern.ch/record/2654870.
  * (38) T. Sjostrand, S. Mrenna, and P. Z. Skands, “PYTHIA 6.4 Physics and Manual,” JHEP 05 (2006) 026, arXiv:hep-ph/0603175 [hep-ph].
  * (39) T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C. O. Rasmussen, and P. Z. Skands, “An Introduction to PYTHIA 8.2,” Comput. Phys. Commun. 191 (2015) 159–177, arXiv:1410.3012 [hep-ph].
  * (40) K. J. Kelly and P. A. N. Machado, “muBHNL.” https://github.com/kjkellyphys/muBHNL.
  * (41) A. de Gouvêa, P. J. Fox, B. J. Kayser, and K. J. Kelly, “Three-body decays of heavy Dirac and Majorana fermions,” Phys. Rev. D 104 (2021) no. 1, 015038, arXiv:2104.05719 [hep-ph].
  * (42) K. J. Kelly and P. A. N. Machado, “MicroBooNE experiment, NuMI absorber, and heavy neutral leptons,” Phys. Rev. D 104 (2021) no. 5, 055015, arXiv:2106.06548 [hep-ph].
  * (43) SHiP Collaboration, “HNL sensitivity of SHiP experiment.” https://doi.org/10.5281/zenodo.1472071.
  * (44) 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].
  * (45) A. Boyarsky, M. Ovchynnikov, O. Ruchayskiy, and V. Syvolap, “Improved big bang nucleosynthesis constraints on heavy neutral leptons,” Phys. Rev. D 104 (2021) no. 2, 023517, arXiv:2008.00749 [hep-ph].