aainstitutetext: Discovery Center, Niels Bohr Institute, University of
Copenhagen,
Blegdamsvej 17, DK-2010, Copenhagen, Denmarkbbinstitutetext: Institute of
Physics, Laboratory for Particle Physics and Cosmology,
École Polytechnique Fédérale de Lausanne, CH-1015 Lausanne, Switzerland

# Dirac vs. Majorana HNLs (and their oscillations) at SHiP

J.-L. Tastet111To whom correspondence should be addressed. b I. Timiryasov
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

SHiP is a proposed high-intensity beam dump experiment set to operate at the
CERN SPS. It is expected to have an unprecedented sensitivity to a variety of
models containing feebly interacting particles, such as Heavy Neutral Leptons
(HNLs). Two HNLs or more could successfully explain the observed neutrino
masses through the seesaw mechanism. If, in addition, they are quasi-
degenerate, they could be responsible for the baryon asymmetry of the
Universe. Depending on their mass splitting, HNLs can have very different
phenomenologies: they can behave as Majorana fermions—with lepton number
violating (LNV) signatures, such as same-sign dilepton decays—or as Dirac
fermions with only lepton number conserving (LNC) signatures. In this work, we
quantitatively demonstrate that LNV processes can be distinguished from LNC
ones at SHiP, using only the angular distribution of the HNL decay products.
Accounting for spin correlations in the simulation and using boosted decision
trees for discrimination, we show that SHiP will be able to distinguish
Majorana-like and Dirac-like HNLs in a significant fraction of the currently
unconstrained parameter space. If the mass splitting is of order
$10^{-6}\;$\mathrm{e}\mathrm{V}$$, SHiP could even be capable of resolving HNL
oscillations, thus providing a direct measurement of the mass splitting. This
analysis highlights the potential of SHiP to not only search for feebly
interacting particles, but also perform model selection.

##  1 Introduction

The experimentally observed non-vanishing neutrino mass differences are among
a few firmly established deviations from the Standard Model (SM) predictions.
An economic way of generating the light neutrino masses is to introduce heavy
singlet fermions with Majorana mass terms into the model Minkowski:1977sc ;
GellMann:1980vs ; Mohapatra:1979ia ; Yanagida:1980xy ; Schechter:1980gr ;
Schechter:1981cv . The masses of the active neutrinos in this extension of the
SM are determined by the type-I seesaw formula and at least two singlet
fermions are needed to accommodate the two observed mass differences of light
neutrinos. A consequence of this mechanism is the presence of heavy Majorana
fermions which mix with active neutrinos. The mass scale of these Majorana
fermions—Heavy Neutral Leptons (HNLs)—is not fixed. It can be below the
electroweak scale,222An argument in favour of the low-scale seesaw comes from
the measured values of the Higgs and top masses. HNLs with masses below the
electroweak scale are not destabilising the Higgs mass Vissani:1997ys ;
Bezrukov:2012sa .  like in the $\mathrm{\nu MSM}$ Asaka:2005an ; Asaka:2005pn
, where two HNLs are responsible for the light neutrino masses and generating
the Baryon Asymmetry of the Universe (BAU) via $CP$-violating oscillations
during their production.

From the FIP (feebly interacting particles) search point of view pbc_2019 ,
HNLs with masses below that of a $B$ meson are the most accessible in the
foreseeable future. There is a vast program to search for HNLs at _intensity
frontier_ experiments, either LHC-based, such as MATHUSLA chou_new_2017 ;
curtin_long-lived_2018 , FASER faser_2018 , CODEX-b gligorov_searching_2018
and AL3X gligorov_leveraging_2018 , or at beam-dump facilities, such as NA62++
drewes_na62_2018  and SHiP ship_collaboration_facility_2015 ;
alekhin_facility_2016 ; ship_collaboration_sensitivity_2018 . The latter two
experiments will be sensitive to both the mass and mixing angles of HNLs.

SHiP is a proposed beam-dump experiment (represented in fig. 1) set to operate
at the CERN SPS. It will use an intense,
$400\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ proton beam from the SPS, dumped
on a thick target in order to produce a large number of heavy hadrons, which
subsequently decay into Standard Model (SM) or feebly-interacting particles.
SHiP is designed to provide a background-free environment to look for the
decays of these heavy FIPs. To this end, a hadron absorber located right after
the target absorbs most SM particles. It is followed by an active muon shield
which deflects the muons away from the experimental cavern. The main detector
consists of a decay volume—evacuated in order to reduce the neutrino
background, and surrounded by vetos—with a tracker and a calorimeter located
at its far end, enabling it to reconstruct the decay event.

Figure 1: Sketch of the SHiP experiment, with the decay chain $H\to
h^{\prime}l_{\alpha}(N\to l_{\beta}h^{\prime\prime})$.

In order to generate the light neutrino masses via the seesaw mechanism, HNLs
must be Majorana fermions, which violate the total lepton number. However, in
some theories—such as the $\mathrm{\nu MSM}$—they can pair to form a coherent
superposition of two quasi-degenerate Majorana fermions, which behaves almost
like a Dirac fermion. Such a combination is dubbed “ _quasi-Dirac pair_ ”. In
this case, the mixing angles can exceed the naive seesaw limit $U^{2}\approx
m_{\nu}/M_{N}$ Shaposhnikov:2006nn ; Kersten:2007vk ; Moffat:2017feq , where
$m_{\nu}$ and $M_{N}$ are respectively the mass scales of light neutrinos and
HNLs. This is possible because a quasi-Dirac fermion approximately conserves
the total lepton number, hence protecting the light neutrino masses. Finally,
the near degeneracy of the HNL masses leads to coherent HNL oscillations. In
the $\mathrm{\nu MSM}$, these oscillations in the early Universe are
responsible for baryogenesis.

For sufficiently light ($\lesssim$10\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$)
HNLs like the ones accessible at SHiP, LNV may be experimentally observable
even when they form a quasi-Dirac pair anamiati_quasi-dirac_2016 ;
drewes_lepton_2019 . We can distinguish three cases,333To be generic, we have
included the more exotic cases of a single Dirac or Majorana HNL. The limits
presented below are for a quasi-Dirac pair, which only differs from those in
the number of events produced.  depending on the scale of the oscillation
phase $\delta M\tau$, where $\delta M$ is the mass splitting of the quasi-
Dirac pair and $\tau$ the typical proper time probed:

  1. 1.

Dirac-like HNL: One Dirac HNL or a quasi-Dirac pair with an oscillation period
exceeding the HNL lifetime or detector size ($\delta M\tau\ll 2\pi$).444As
pointed out in drewes_lepton_2019 , for most experiments, this possibility
might be technically unnatural due to the very small mass splitting needed to
satisfy the inequality. Only LNC processes can be observed.

  2. 2.

Majorana-like HNL: One Majorana HNL or a quasi-Dirac pair with a lifetime and
detector size exceeding the oscillation period ($\delta M\tau\gg 2\pi$). Both
LNC and LNV processes can be observed, with equal integrated rates (see
section 2.2).

  3. 3.

Manifestly quasi-Dirac HNLs: An interesting case occurs when the oscillation
period is comparable to the HNL lifetime or to the size of the
detector555Interestingly, the mass difference needed to generate DM in the
$\mathrm{\nu MSM}$, as found in ref. Canetti:2012kh , is exactly in this
borderline range.  ($\delta M\tau\sim 2\pi$): the experiment may then be
sensitive to the coherent oscillations of HNLs.

If HNLs were to be observed at SHiP, the detection or non-observation of
lepton number violation and HNL oscillations would allow constraining models
and their parameters. The most relevant LNV process at SHiP is the well-
studied _same-sign dilepton decay_ : $H\to[h^{\prime}]l_{\alpha}^{+}(N\to
h^{\prime\prime}l_{\beta}^{+})$, where $H$, $h^{\prime}$ and
$h^{\prime\prime}$ are hadrons (with $h^{\prime}$ possibly missing), and
$l_{\alpha}^{+}$, $l_{\beta}^{+}$, $\alpha,\beta=e,\mu,\tau$ are charged
leptons of potentially different generations. Due to suppressed background,
this type of signature is a smoking gun for HNLs in accelerator searches.
However, at beam-dump experiments, the heavy hadron decay which produces the
HNL takes place inside the target, and therefore the charge of the primary
lepton $l_{\alpha}$ cannot be observed. Naively, it seems that the information
about the HNL production is lost, since the charge of the secondary lepton
$l_{\beta}$, by itself, is not enough to tell apart LNC and LNV processes. As
we shall see in this paper, it turns out that the HNL decay products
nevertheless carry important information. Namely, their distribution is
different for LNC and LNV processes. Not only does this allow distinguishing
Majorana-like from Dirac-like HNLs given sufficiently many events, but the
knowledge of these distributions can also be used to _resolve_ HNL
oscillations and directly measure the mass splitting.

Estimating these two distributions is complicated by the presence of a variety
of two- and three-body production channels. In addition, the parent hadrons
are produced with a finite spectrum. As we shall see in section 3.3, this
smears the distributions, making them look more similar. Therefore, in order
to assess whether SHiP will be able to discriminate between Majorana- and
Dirac-like HNLs, an accurate treatment of all production channels, including
spin correlations, is required. This is accomplished using a Monte-Carlo
simulation.

The angular distribution of HNL decay products has been studied in a collider
setting for decays which are not fully reconstructible arbelaez_probing_2018 ;
dib_signatures_2017 ; 1810.07210  (such as trilepton decays), as well as for
beam-dump experiments cvetic_probing_2012 ; balantekin_addressing_2019 . Our
analysis improves on the latter by not relying on HNLs being produced as
helicity eigenstates, by handling a larger class of production channels, by
considering the full phase-space distribution of the HNL decay products
(instead of just their energy) and by producing a concrete sensitivity
estimate using a realistic geometry and heavy meson spectrum for SHiP.

This paper is organized as follows. In section 2, we review the Standard Model
extended with HNLs, and discuss lepton number violation and coherent HNL
oscillations. In section 3, we analyze the different signatures of LNC and LNV
processes at the SHiP experiment. In section 4, we propose a strategy to
detect LNV and reconstruct HNL oscillations. Finally, in section 5, we present
the sensitivity of SHiP to LNV achieved through this method, as well as a
possible signature of HNL oscillations. Technical details about the simulation
and the statistical analysis are respectively provided in appendices A and B.

##  2 Model

###  2.1 Heavy Neutral Leptons

We consider the Standard Model extended with $\mathcal{N}$ HNLs $N_{I}$, which
are spin-$\frac{1}{2}$ SM singlets with Majorana masses $M_{I}$, and new
Yukawa couplings $Y_{\alpha I}^{\nu}$, with $\alpha=e,\mu,\tau$ the lepton
flavor index. Using the conventions from dreiner_two-component_2010 :

$\mathcal{L}=\mathcal{L}_{\mathrm{SM}}+\frac{i}{2}N_{I}^{\dagger}(\bar{\sigma}\cdot\partial)N_{I}-(Y_{\alpha
I}^{\nu})^{*}(\phi\cdot
L_{\alpha})N_{I}-\frac{M_{I}}{2}N_{I}N_{I}+\mathrm{h.c.}$ (1)

After electroweak symmetry breaking, the Yukawa interaction generates a Dirac
mass term $(m_{D})_{\alpha I}=\frac{v}{\sqrt{2}}(Y_{\alpha I}^{\nu})^{*}$,
resulting in a non-diagonal, symmetric Dirac-Majorana mass term for neutrinos
Giunti:2007ry :

$\\!M_{\mathrm{DM}}=-\frac{1}{2}\begin{pmatrix}\nu^{T}&N^{T}\end{pmatrix}\begin{pmatrix}0&m_{D}^{T}\\\
m_{D}&M_{M}\end{pmatrix}\begin{pmatrix}\nu\\\ N\end{pmatrix}+\mathrm{h.c.}$
(2)

where $M_{M}=\mathrm{diag}\left(M_{I}\dots\right)$. Using a unitary
transformation of the fields (Takagi factorization takagi_1927 ), the mass
matrix can be brought to a diagonal form:

$\displaystyle\nu_{\alpha}=U_{\alpha i}n_{i}\quad\text{ and }\quad
N_{I}=U_{Ii}n_{i}$ (3) $\displaystyle
M_{\mathrm{DM}}=-\frac{m_{i}}{2}(n_{i}n_{i}+n_{i}^{\dagger}n_{i}^{\dagger})$
(4)

In the limit $|M_{M}|\gg|m_{D}|$, we can use an approximate block
factorization, leading to the mass eigenstates $n_{i}\cong\nu_{i},N_{I}$
mixing with the flavor fields as:

$\displaystyle\nu_{\alpha}$ $\displaystyle\cong U_{\alpha
i}^{\mathrm{PMNS}}\nu_{i}+\Theta_{\alpha I}N_{I}$ (5)
$\displaystyle\Theta_{\alpha I}$ $\displaystyle\cong M_{I}^{-1}(m_{D})_{\alpha
I}$ (6)

and the following mass sub-matrices:

$\displaystyle m_{\alpha\beta}$
$\displaystyle\cong-\sum_{I}\frac{(m_{D})_{\alpha I}(m_{D})_{\beta
I}}{M_{I}}\cong-\sum_{I}M_{I}\Theta_{\alpha I}\Theta_{\beta I}$ (7)
$\displaystyle m_{IJ}$ $\displaystyle\cong M_{I}\delta_{IJ}$ (8)

The choice of the mass scale $M_{M}$ and Yukawa couplings $Y_{\alpha I}^{\nu}$
is not uniquely dictated by low-energy neutrino observables, and should be
fixed otherwise.

The Standard Model features an accidental symmetry—lepton number—which, at
tree level, is conserved for massless or Dirac neutrinos, but is violated by
the Majorana mass term of HNLs. Charged leptons and neutrinos have lepton
number $+1$, while charged anti-leptons and anti-neutrinos have lepton number
$-1$. If lepton number is conserved (LNC), then the only allowed Feynman
diagrams are those with a conserved flow of lepton number (represented by the
arrow on the fermion lines of leptons), like the opposite-sign dilepton decay
of a heavy hadron shown in fig. 2(a). On the other hand, in the presence of
lepton number violating (LNV) operators, processes like the same-sign dilepton
decay shown in fig. 2(b) become possible. Lepton number violation can also
manifest itself in neutral-current processes or in neutrinoless double-$\beta$
decay. Whether such LNV transitions actually happen depends on the specific
model.

(a) LNC

(b) LNV

Figure 2: Lepton number conserving and violating decay chains for
$H\rightarrow h^{\prime}l_{\alpha}(N\rightarrow l_{\beta}h^{\prime\prime})$.

In the past decade, a class of low-scale seesaw models have risen in
popularity, such as the $\nu$MSM Asaka:2005pn , not least because of their
falsifiability at existing or proposed experiments. In these models, $M_{M}$
is postulated to be below the electroweak scale. The seesaw formula (7)
requires at least $2$ HNLs to explain the two observed mass differences. If
their parameters are arbitrary, then the smallness of the light neutrino
masses is achieved through small Yukawa couplings of order
$Y^{\nu}\sim\frac{1}{v}\sqrt{|m_{\nu}||M_{M}|}$, leading to squared mixing
angles $|\Theta|^{2}\sim|m_{\nu}|/|M_{M}|$. For a typical HNL with
$M_{M}\sim$1\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$, this gives
$|\Theta|^{2}\sim 10^{-11}$, a number that is too small to be probed at any
current or proposed experiment.

However, multiple HNLs can have mixing angles well above the seesaw limit, yet
at the same time produce the correct neutrino masses in a technically natural
way, if a certain symmetry is imposed on their Yukawa couplings. If we
consider for simplicity $\mathcal{N}=2$ nearly degenerate HNLs $N_{1,2}$,
their mixing angles should be related by $\Theta_{\alpha 2}\approx\pm
i\Theta_{\alpha 1}$ Shaposhnikov:2006nn ; Kersten:2007vk . Such HNLs form a
quasi-Dirac fermion, which approximately conserves the total lepton number.
This implies that the usual searches for naive LNV effects (e.g. same-sign
dilepton decays), may return null results even if HNLs _are there_.

Below we discuss an important consequence of the approximate nature of this
lepton number conservation: HNL oscillations, and how quasi-Dirac HNLs can
phenomenologically behave either as Majorana or Dirac HNLs depending on their
mass splitting $\delta M$ and the length scale probed at the experiment.

###  2.2 Coherent oscillations of Heavy Neutral Leptons

The SHiP experiment is only sensitive to
$\mathrm{G}\mathrm{e}\mathrm{V}$-scale HNLs, with mixing angles significantly
above the seesaw limit ship_collaboration_sensitivity_2018 . Therefore it can
only probe the quasi-Dirac regime described above. Apart from a small mass
splitting $\delta M\ll M$, the two HNLs are otherwise identical. Since these
two HNLs cannot be distinguished in any realistic experiment, they both
mediate the same processes and each contribute to the total transition
amplitude, resulting in interference. Only the initial and final-state
particles, which strongly interact with the environment, are measured in the
quantum mechanical sense. In order to accurately describe processes involving
multiple HNLs, it is therefore necessary to consider them as intermediate
particles within a larger process consisting of the HNL production,
propagation and decay, and only square the overall transition amplitude
between the observed, external particles. This can be formulated rigorously
within the framework of the external wave packet model sachs_interference_1963
; giunti_treatment_1993  (see also beuthe_oscillations_2003 ;
akhmedov_neutrino_2010 ; akhmedov_neutrino_2011 ; akhmedov_quantum_2019  and
references therein for recent reviews). Let us note in passing that this
description automatically takes care of spin correlations between the
particles taking part in the HNL production and decay.

In what follows, we consider a typical reconstructible decay chain at SHiP, as
depicted in fig. 2. We will postpone the detailed discussion of this process
to section 3. A heavy hadron $H$ produced in the target decays at space-time
coordinates $x_{P}$ into an HNL $N_{I}$, a charged lepton $l_{\alpha}$ (the
_primary_ lepton), and an optional hadron $h^{\prime}$. If the HNL is
sufficiently long-lived, it can propagate a macroscopic distance before
decaying at $x_{D}$ into a charged lepton $l_{\beta}$ (the _secondary_ lepton)
and a hadron $h^{\prime\prime}$.

The slightly different masses of the HNLs mediating the process lead to
different dispersion relations $q_{I}^{2}=M_{I}^{2}$. As a consequence, the
space-time-dependent phase $e^{-iq_{I}\cdot(x_{D}-x_{P})}$ acquired by the HNL
between its production and decay will differ slightly for each mass
eigenstate. When squaring the amplitude in order to obtain the differential
decay rate, the interference terms between the partial amplitudes coming from
different mass eigenstates will therefore feature a space-time-dependent
modulation: HNL oscillations. The external wave packet model allows one to
unambiguously establish the expression for the oscillation phase and check
that the entire process remains coherent in all experimentally relevant
situations.

The present paper does not aim to be a detailed study of HNL oscillations,
which have already been covered in various settings and limits in the
literature Asaka:2005pn ; eijima_parameter_2018 ; 1409.4265 ; 1505.04749 ;
anamiati_quasi-dirac_2016 ; 1709.03797 ; 1805.00070 ; 1810.07210 . Therefore,
we will only quote the main result. Let
$\mathrm{d}\hat{\Gamma}_{\alpha\beta}^{\pm\pm}$ be the differential decay rate
for the above-described process
$H\rightarrow[h^{\prime}]l_{\alpha}^{\pm}(N\rightarrow
l_{\beta}^{\pm}h^{\prime\prime})$ with a single HNL, in the (unphysical) limit
of unit mixing angles between the HNL and all three active flavors, and
without the absorptive part. The coherent differential width
$\smash{\mathrm{d}\Gamma_{\alpha\beta}^{\pm\pm}(\tau)}$ for $\mathcal{N}$
HNLs, as a function of the proper time $\tau=\sqrt{(x_{D}-x_{P})^{2}}$ between
the HNL production and decay vertex, is then:

$\mathrm{d}\Gamma_{\alpha\beta}^{\pm\pm}(\tau)=\left|\sum_{I=1}^{\mathcal{N}}\Theta_{\alpha
I}^{\pm}\Theta_{\beta
I}^{\pm}e^{-iM_{I}\tau-\frac{\Gamma_{I}}{2}\tau}\right|^{2}\mathrm{d}\hat{\Gamma}_{\alpha\beta}^{\pm\pm}$
(9)

where $M_{I}$ is the mass of the $I$-th heavy mass eigenstate, $\Gamma_{I}$
its total width, and we have used the shorthand notation
$\Theta^{+}\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny def}}}{=}}\Theta^{*}$
and $\Theta^{-}\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny
def}}}{=}}\Theta$.

In the case of $\mathcal{N}=2$ HNLs forming a quasi-Dirac pair, i.e.
$M_{1}=M-\frac{\delta M}{2}$, $M_{2}=M+\frac{\delta M}{2}$, $\Theta_{\alpha
2}\cong\pm i\Theta_{\alpha 1}$ and
$\Gamma_{1}\cong\Gamma_{2}\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny
def}}}{=}}\Gamma$, the coherent differential width reduces to:

$\mathrm{d}\Gamma_{\alpha\beta}^{\pm\pm}(\tau)\cong 2\left|\Theta_{\alpha
1}\right|^{2}\left|\Theta_{\beta 1}\right|^{2}\left(1\pm\cos\left(\delta
M\tau\right)\right)e^{-\Gamma\tau}\mathrm{d}\hat{\Gamma}_{\alpha\beta}^{\pm\pm}$
(10)

where the $+$ sign is for lepton number conserving processes
($\mathrm{d}\Gamma_{\alpha\beta}^{+-}$ and
$\mathrm{d}\Gamma_{\alpha\beta}^{-+}$), and the $-$ sign for lepton number
violating ones ($\mathrm{d}\Gamma_{\alpha\beta}^{++}$ and
$\mathrm{d}\Gamma_{\alpha\beta}^{--}$). Notice how in the quasi-Dirac limit,
the oscillation pattern does not explicitly depend on the lepton flavors
$\alpha$ and $\beta$, but only on whether the process is LNC or LNV. If
$\delta M$ vanishes exactly, HNLs form a Dirac fermion and LNV effects are
completely absent. Recently, $CP$-violating HNL oscillations have attracted
some interest 1904.04787 ; 1905.03097 ; 1906.09470 ; abada_interference_2019 .
However, here we can see that $CP$-violation is suppressed in the quasi-Dirac
limit.

Throughout this paper, we will focus on the case where $\Gamma\tau\ll 1$,
which is the most relevant for SHiP, and drop the exponentially decaying
factor. Analysing formula (10), we see that there are three regimes of
interest, depending on the mass splitting $\delta M$ and proper time scale
$\tau$ probed at the experiment:

  * •

If $\delta M\tau\ll 2\pi$, the HNL pair is observed before the onset of
oscillations, and it behaves like a single Dirac HNL, i.e. we cannot observe
lepton-number violation.

  * •

If $\delta M\tau\gg 2\pi$, fast oscillations are averaged out, and the HNL
pair behaves like a single Majorana HNL, with equal integrated decay rates for
LNC and LNV channels.666 In the rest frame of a single on-shell, Majorana HNL,
the only “memory” of the production process is the HNL spin. To perform the
phase-space integration for the HNL decay, one can always choose a frame where
the HNL is at rest and with a fixed spin projection, hence resulting in the
same integrated rates for LNC and LNV processes.

  * •

If $\delta M\tau\sim 2\pi$, oscillations must be accounted for. If it is
possible to experimentally reconstruct, for each selected event, the proper
time $\tau$ between the production and decay vertex of the HNL, then
oscillations can be resolved, i.e. the $\tau$-differential event rates for LNC
/ LNV will show a periodic modulation according to eq. 10.

At SHiP, the proper time scale $\tau$ is about $2\text{\,}\mathrm{m}$ for
sufficiently long-lived HNLs. It corresponds to the average time between the
production and decay of an _observed_ HNL, in its rest frame. Therefore, the
critical mass splitting separating the three regimes—near which oscillations
are resolvable—is about $10^{-6}\;$\mathrm{e}\mathrm{V}$$.

##  3 Probing lepton number violation at SHiP

Many collider searches for Majorana HNLs lhcb_collaboration_search_2014 ;
atlas_collaboration_search_2015 ; cms_collaboration_search_2015 ;
cms_collaboration_search_2016  are sensitive to lepton number violation
through the charges of the leptons produced at the HNL production and decay
vertex. Indeed, due to the chiral nature of the weak interaction, they
unambiguously tell the chiral projection through which the HNL interacts at a
given vertex. In theory, a same-sign dilepton decay (either prompt or
displaced) would thus provide clear evidence for lepton number violation
(although, in practice, significant standard model backgrounds exist for
prompt decays).

At SHiP, similar numbers of mesons and anti-mesons are expected to be
produced.777Unless cascade production significantly alters the results from
lebc_d-meson_1987 . The charm spectrum will be measured at SHiP prior to data
taking ship_collaboration_measurement_2017 . Asymmetries, if present, can only
improve the classification accuracy, since the secondary lepton charge would
then carry _some_ information.  This leads to similar numbers of HNLs being
produced along with positively and negatively charged primary leptons.
Consequently, the secondary lepton charge contains very little information as
to whether the process is LNC or LNV. To lift this degeneracy, it becomes
necessary to look at new observables.

Luckily, the HNL lepton number is not the only quantum number conserved by the
weak interaction. The HNL also carries spin $\frac{1}{2}$, and the total
angular momentum is always conserved. When the HNL is produced, its spin is
correlated (opposite if $H$ and $h^{\prime}$ are pseudoscalar) with that of
the primary lepton. Due to chiral suppression, the spin of the primary lepton
is itself correlated with its lepton number (see for example the left part of
fig. 3). This suggests that by looking at the angular distribution of the
secondary particles—which may be observable—we should be able to obtain
information about the primary interaction, and thus whether the process was
LNC or LNV (see the right part of fig. 3). This realization was the starting
point of the present work. More generally, we expect LNC and LNV decay chains
to have different kinematics due to their different Lorentz structures,
potentially allowing us to distinguish them without directly observing the
primary decay.

In section 3.1, we describe the relevant HNL production and decay channels at
SHiP; in section 3.2, we quantitatively compare the angular distributions for
LNC and LNV processes, and in section 3.3 we discuss how this affects the
observable momenta in a beam-dump setting.

Figure 3: This sketch explains the origin of the different angular
correlations for LNC and LNV processes. For simplicity, here we consider two-
body primary and secondary decays involving only pseudoscalar mesons, and the
masses of the charged leptons and of $h^{\prime\prime}$ are neglected. For
definiteness, the charge of the primary lepton—which is produced inside the
target and thus inaccessible—is also fixed to $+$. Since the HNL is a Majorana
fermion, the secondary lepton $l_{\beta}$ can have either charge. However, due
to angular momentum conservation, the lepton $l^{+}_{\alpha}$ and the HNL $N$
are produced with opposite spin projections in the rest frame of the heavy
meson $H$. Because of chiral suppression (which is more effective for light
fermions), the charge of the primary lepton is correlated with its spin (e.g.
in the massless limit, $l_{\alpha}^{+}$ has helicity $+\frac{1}{2}$) and hence
with the HNL spin. For the same reason, the angular distribution of the decay
products of the resulting HNL spin eigenstate (which is unaffected by a boost
along the quantization axis) will therefore depend on the secondary lepton
charge. The very same formula for the probability $\mathcal{P}$ also holds for
$CP$-conjugated channels, with the $+$ sign for LNC and the $-$ sign for LNV.
The general case (massive, with two- or three-body primary decay) is discussed
in section 3.2.

###  3.1 HNL production and decay at SHiP

At SHiP, most HNLs are produced in heavy meson decays through flavor-changing
charged currents, as discussed in bondarenko_phenomenology_2018 . In addition,
for the present analysis, we will only consider fully reconstructible HNL
decays such as $N\rightarrow l_{\beta}^{\mp}\pi^{\pm}$, producing only charged
particles which are sufficiently long-lived to be detected by the tracking
station located at the end of the decay vessel. Those are also mediated by the
charged-current interaction.

Without losing generality, we can therefore consider the generic lepton number
conserving and violating processes
$H\rightarrow[h^{\prime}]l_{\alpha}(N\rightarrow l_{\beta}h^{\prime\prime})$
represented in figs. 2(a) and 2(b), respectively, as well as their
$CP$-conjugates. $H$ denotes a heavy hadron (typically a $D_{[s]}$ or
$B_{[c]}$ meson at SHiP), $h^{\prime}$ and $h^{\prime\prime}$ are hadrons
(with $h^{\prime}$ missing for two-body primary decays), and
$l_{\alpha}^{\pm}$ and $l_{\beta}^{\pm}$ are respectively the primary and
secondary leptons.

Since the heavy hadron $H$ is typically short-lived, the primary decay takes
place inside the target and cannot be observed. If the HNL is sufficiently
long-lived (we will assume this to be the case throughout this paper), it can
propagate a macroscopic distance before decaying, and leave a very displaced
vertex inside the SHiP decay vessel. For the selected decay channels
$N\rightarrow l_{\beta}^{\mp}\pi^{\pm}$, this secondary vertex can be fully
reconstructed.

In the present study, we will restrict ourselves to HNL masses between the $K$
and $D_{s}$ thresholds. Masses below the $K$ threshold have already been
heavily constrained pbc_2019 , while above the $D_{s}$ mass, HNLs are mainly
produced in $B$ meson decays, whose spectrum cannot be directly measured at
the beam dump, making our analysis more sensitive to modeling errors.

###  3.2 Angular correlations in LNC and LNV decay chains

In order to study the angular correlations between all final-state particles,
spin correlations between the primary and secondary decay must be accounted
for. Those result from the non-observation of the HNL spin, which leads to
interference between the two spin eigenstates $N_{s}$, $s=\pm\frac{1}{2}$
(similarly to how the non-observation of its precise mass allows for flavor
oscillations). To compute the overall transition amplitude, we can therefore
use the same trick as for oscillations, i.e. treat the primary and secondary
decays as a single process.

To simplify the calculations, in this section we will focus on the case of a
single Majorana HNL, which mediates both LNC and LNV decay chains with equal
rates, and we will omit the absorptive part of the amplitude (i.e. we will
study $\mathrm{d}\hat{\Gamma}_{\alpha\beta}^{\pm\pm}$ instead of
$\mathrm{d}\Gamma_{\alpha\beta}^{\pm\pm}(\tau)$). We do not lose generality in
doing so, because the effect of multiple nearly degenerate HNLs and their
finite lifetime can be factored out, and subsequently recovered, using eqs. 9
and 10. To keep the notation light, we will from now on drop the HNL index
$I=1$.

Since we are only concerned with long-lived HNLs, which are produced on their
mass shell and have well separated, localized production and decay vertices,
the momentum $q$ of the HNL is practically fixed, which allows factorizing the
transition amplitude as:

$\mathcal{A}\left(H\rightarrow
h^{\prime}l_{\alpha}l_{\beta}h^{\prime\prime}\right)\Bigr{|}_{N\text{ long-
lived}}\propto\sum_{s=\pm\frac{1}{2}}\mathcal{A}\left(H\rightarrow
h^{\prime}l_{\alpha}N_{s}(q)\right)\mathcal{A}\left(N_{s}(q)\rightarrow
l_{\beta}h^{\prime\prime}\right)$ (11)

where we have omitted the complex phase $e^{-iq\cdot(x_{D}-x_{P})}$ resulting
from the HNL propagation, which is unimportant in the case of one HNL. The
sub-amplitudes for the primary and secondary polarized decays are then
straightforward to compute using the usual Feynman rules with two-component
spinors dreiner_two-component_2010 .

Consider now the LNC and LNV processes
$H\rightarrow[h^{\prime}]l_{\alpha}(N\rightarrow l_{\beta}h^{\prime\prime})$
where $H,h^{\prime},h^{\prime\prime}$ are pseudoscalar mesons and $h^{\prime}$
may be missing. They are respectively represented in figs. 2(a) and 2(b), with
the arrows denoting the flow of lepton number. Their $CP$-conjugates have been
omitted, since in the absence of oscillations (as is the case for the
incoherent width), $CP$ is conserved. As can be seen in fig. 5, the primary
decays $H\to[h^{\prime}]l_{\alpha}N$ with $h^{\prime}$ a pseudoscalar meson or
missing indeed produce the majority of HNLs with masses
$\gtrsim$0.7\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$ and below the $D_{s}$
mass.888Below $M_{N}\approx$0.7\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$, a
non-negligible fraction of HNLs is produced along with a _vector_ meson. In
this case, we expect the angular correlations to reverse compared to the
pseudoscalar case.  Let $J_{W\mu}^{h}$ be the hadronic charge-lowering
current,
$j_{1\mu}^{-}=\smash{\big{\langle}h^{\prime}|J_{W\mu}^{h}|H\big{\rangle}}$ and
$j_{2\mu}^{\mp}=\smash{\big{\langle}h^{\prime\prime}|J_{W\mu}^{h(\dagger)}|0\big{\rangle}}$
the hadronic matrix elements, $p_{\alpha,\beta}$ the charged lepton momenta,
and $q$ the HNL momentum. If the primary decay is purely leptonic, then
$|h^{\prime}\rangle=|0\rangle$. Since SHiP cannot directly measure the spin or
helicity of the particles detected, we sum incoherently over all possible spin
configurations of final state particles. The spin-summed, squared amplitudes
are then, in the Fermi approximation:

$\displaystyle\overline{\big{|}\mathcal{A}_{\mathrm{LNC}}(H\rightarrow
h^{\prime}l_{\alpha}^{+}l_{\beta}^{-}h^{\prime\prime})\big{|}^{2}}$
$\displaystyle=\frac{\left|\Theta_{\alpha}\right|^{2}\left|\Theta_{\beta}\right|^{2}}{v^{8}}\,\mathrm{tr}\left(P_{R}\not{p}_{\alpha}\not{j}_{1}^{*}\not{q}\not{j}_{2}^{*}\not{p}_{\beta}\not{j}_{2}\not{q}\not{j}_{1}\right)$
(12) $\displaystyle\overline{\big{|}\mathcal{A}_{\mathrm{LNV}}(H\rightarrow
h^{\prime}l_{\alpha}^{+}l_{\beta}^{+}h^{\prime\prime})\big{|}^{2}}$
$\displaystyle=\frac{\left|\Theta_{\alpha}\right|^{2}\left|\Theta_{\beta}\right|^{2}}{v^{8}}M_{N}^{2}\,\mathrm{tr}\left(P_{R}\not{p}_{\alpha}\not{j}_{1}^{*}\not{j}_{2}^{*}\not{p}_{\beta}\not{j}_{2}\not{j}_{1}\right)$
(13)

where we have omitted the ± for brevity if they can be inferred from context,
$\Theta_{\alpha,\beta}$ are the mixing angles, and
$v=\langle\left|\phi\right|\rangle{}\approx$246\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$
is the vacuum expectation value of the Higgs field. These results are
consistent with the polarized decay rates from ballett_heavy_2019 , but
generalize to the case where the primary decay produces a superposition of HNL
helicity eigenstates. The above two expressions differ in the trace, therefore
we generically expect them to produce different momentum distributions for LNC
and LNV processes. However, in their current form, this difference is not
manifest. To understand it, it is interesting to consider the special case
where the production process is a two-body decay. As can be seen in figs. 5
and 5, it is actually the main production channel for HNLs with masses
$\gtrsim$1\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$ and below the $D_{s}$
mass.

Figure 4: Number of HNLs produced at SHiP as a function of the primary decay
multiplicity, for a coupling to one flavor.

Figure 5: Fraction of HNLs produced at SHiP as a function of the primary decay
multiplicity and spin of the outgoing meson, for a coupling to one flavor.

When both the production and decay process are two-body decays, the hadronic
matrix elements are $j_{1}^{\mu}=-iV_{UD}f_{H}p_{H}^{\mu}$ and
$j_{2}^{\mu}=+iV_{U^{\prime}D^{\prime}}f_{h^{\prime\prime}}p_{h^{\prime\prime}}^{\mu}$,
where $V_{UD}$ denotes the relevant CKM matrix element and $f_{h}$ is the
meson decay constant. Neglecting the masses of the final state particles,
which give
$\mathcal{O}\left(\nicefrac{{m_{\alpha,\beta,h^{\prime\prime}}^{2}}}{{M_{H,N}^{2}}}\right)$
corrections, the traces from eqs. 12 and 13, respectively for LNC and LNV
processes, simplify to:

$\displaystyle\mathrm{tr}\left(P_{R}\not{p}_{\alpha}\not{j}_{1}^{*}\not{q}\not{j}_{2}^{*}\not{p}_{\beta}\not{j}_{2}\not{q}\not{j}_{1}\right)$
$\displaystyle\cong\left|V_{UD}\right|^{2}\left|V_{U^{\prime}D^{\prime}}\right|^{2}f_{H}^{2}f_{h^{\prime\prime}}^{2}\cdot
M_{N}^{6}\left(M_{H}^{2}-M_{N}^{2}-s_{ll}\right)$ (14) $\displaystyle
M_{N}^{2}\mathrm{tr}\left(P_{R}\not{p}_{\alpha}\not{j}_{1}^{*}\not{j}_{2}^{*}\not{p}_{\beta}\not{j}_{2}\not{j}_{1}\right)$
$\displaystyle\cong\left|V_{UD}\right|^{2}\left|V_{U^{\prime}D^{\prime}}\right|^{2}f_{H}^{2}f_{h^{\prime\prime}}^{2}\cdot
M_{N}^{6}s_{ll}$ (15)

where $s_{ll}\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny
def}}}{=}}(p_{\alpha}+p_{\beta})^{2}$ is the invariant dilepton mass. Note the
linear and opposite dependences of the LNC and LNV spin-summed squared
amplitudes on $s_{ll}$. To understand their origin, it is enlightening to
reexpress $s_{ll}$ in the rest frame of the HNL, in terms of the angle
$\theta_{ll}^{\mathrm{CM}}=\angle(\mathbf{p}_{\alpha}^{\mathrm{CM}},\mathbf{p}_{\beta}^{\mathrm{CM}})$
between the two lepton momenta. Still in the massless limit, we find:

$s_{ll}=\frac{M_{H}^{2}-M_{N}^{2}}{2}\Big{(}1-\cos\big{(}\theta_{ll}^{\mathrm{CM}}\big{)}\Big{)}$
(16)

Therefore,

$\displaystyle\overline{\big{|}\mathcal{A}_{\mathrm{LNC}}\big{|}^{2}}$
$\displaystyle\propto 1+\cos\big{(}\theta_{ll}^{\mathrm{CM}}\big{)}$ (17)
$\displaystyle\overline{\big{|}\mathcal{A}_{\mathrm{LNV}}\big{|}^{2}}$
$\displaystyle\propto 1-\cos\big{(}\theta_{ll}^{\mathrm{CM}}\big{)}$ (18)

We observe that opposite-sign leptons (LNC) tend to be produced in the same
direction, and same-sign leptons (LNV) in opposite directions. As explained in
fig. 3, this is a consequence of the chirality of the weak interaction and the
conservation of the total angular momentum. In the absence of any other
dynamics, spin projections lead to the characteristic angular dependence in
$\cos\left(\nicefrac{{\theta_{ll}^{\mathrm{CM}}}}{{2}}\right)$ and
$\sin\left(\nicefrac{{\theta_{ll}^{\mathrm{CM}}}}{{2}}\right)$ of the
transition amplitude, respectively for LNC and LNV. Equations 17 and 18 then
directly follow from squaring the amplitude.

In the massive case, the finite masses of the decay products can results in
helicity flips, and in the three-body case, the QCD matrix elements result in
non-trivial correlations between the momenta of the primary decay products.
These effects complicate the correlations between the various momenta.
Nevertheless, they can be accounted for when _sampling_ events. To this end,
we have implemented the full matrix elements from eqs. 12 and 13 in our Monte-
Carlo simulation, as discussed in section A.4.

###  3.3 Angular distribution in the laboratory frame

At SHiP, the invariant mass $s_{ll}$ (or angle $\theta_{ll}^{\mathrm{CM}}$)
cannot be reconstructed. This is because neither the heavy hadron momentum nor
the momenta of its decay products (other than the HNL) can be determined.
Indeed, the heavy hadrons producing the HNLs do not have a monochromatic
spectrum, and the primary decay cannot be observed since it takes place inside
the target. One can then reasonably wonder if some difference between the LNC
and LNV distributions subsists when looking only at the (observable) secondary
decay products, in the laboratory frame, or if it is washed out.

Figure 6: This sketch shows how the different distributions of $l_{\beta}$ in
the HNL rest frame for LNC vs. LNV processes affect the corresponding
distributions in the rest frame of the heavy hadron $H$ and in the laboratory
frame. The various momenta shown for $l_{\beta}$ represent multiple
realizations of the decay. In the $H$ frame, LNV processes typically result in
larger momenta for $l_{\beta}$ than LNC ones. In the laboratory frame, this
effect partly survives the averaging over the heavy hadron spectrum and
manifests itself as a broadening of the distribution of the secondary lepton
momentum $p_{\beta}$.

To start answering this question, it is instructive to go back to the
simplified case discussed in section 3.2, where the HNL is produced and decays
through two-body processes involving pseudoscalar mesons. In the HNL rest
frame, we obtained the following correlation: for LNV processes, the direction
of the secondary lepton momentum is positively correlated with the boost
direction (denoted by $z$ on figs. 3 and 6) from the heavy meson rest frame to
the HNL rest frame; while for LNC processes it is anti-correlated. This is
depicted in the left panel of fig. 6. Furthermore, in two-body decays, the
magnitudes of all momenta in the rest frame of the parent particle are fixed
by four-momentum conservation, and depend only on the particle masses.
Consequently, in the heavy meson rest frame, the magnitude of the secondary
lepton momentum will on average be larger for LNV processes compared to LNC
ones. This argument is still valid for three-body decays involving
pseudoscalar mesons. A non-trivial asymmetry thus subsists in the heavy meson
rest frame (see the middle panel of fig. 6).

As a final step, the momenta must be boosted back to the laboratory frame.
Since the heavy hadron momentum is not fixed, this has the potential to wash
out the correlations. At SHiP, heavy mesons have a large momentum spread along
the beam axis ($\mathcal{O}($10\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$)$,
much larger than the yield of the meson decay), and a significantly smaller
one ($\mathcal{O}($1\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$)$) in the
transverse direction (see section A.3). The asymmetry between the LNC and LNV
distributions is therefore more likely to be visible in the transverse plane
than along the beam axis. For it to be significant, the HNL kinetic energy in
the heavy hadron rest frame should be similar to or exceed the transverse
momentum spread of the hadron spectrum. As a result, we expect the $p_{T}$
spectrum of the secondary lepton $l_{\beta}$ to be broader for LNV processes
than for LNC ones (see the right panel of fig. 6), provided that both of them
are broader than the irreducible $p_{T}$ spread of the heavy meson spectrum.

Alternatively, one could try to approximate the angle
$\theta_{ll}^{\mathrm{CM}}$ in the HNL rest frame. If the heavy hadron
momentum is fixed, this can be done exactly, and results in the maximal
classification accuracy allowed by spin projections (e.g. $a=3/4$ in the two-
body, massless case). It is then equivalent to measuring the (observable)
momentum $p^{\mathrm{CM}}$ of the secondary lepton $l_{\beta}$ in the HNL rest
frame. However, when the heavy hadron has a finite spectrum, the boost
direction from its rest frame to the HNL rest frame is not fixed any more.
This partially decorrelates $\theta_{ll}^{\mathrm{CM}}$ and $p^{\mathrm{CM}}$,
hence reducing the discriminating power of the latter.

As we shall see in section 4.2, the features discussed above can indeed be
used to discriminate between LNC and LNV processes (see for example fig. 7).
More generally, any difference—in the laboratory frame—between the
distributions of the visible decay products of LNC and LNV processes opens up
the possibility of measuring their relative rates, given sufficiently many
events. Although discriminating between these two classes of events would be
very challenging analytically, this problem is well suited to multivariate
analysis.

Further complications arise, however, due to HNLs being produced from a mix of
various two- and three-body decays, and because of the geometrical acceptance
of the experiment, which alters the distribution of visible particles.
Generating a training set which faithfully reproduces the angular correlations
discussed above while including these effects is therefore best done using a
Monte-Carlo simulation. In the next section, we discuss the simulation used to
generate the training set (section 4.1), then how we use it to train a binary
classifier (section 4.2), and finally how we use the classifier output in
order to perform model selection (section 4.3) and reconstruct HNL
oscillations (section 4.4).

##  4 Simulation and analysis

###  4.1 Simulation

In order to accurately estimate the distribution of the momenta of the HNL
decay products, we have devised a simple Monte-Carlo simulation, which
generates the primary and secondary decays at once, using the matrix elements
presented in section 3.2. The first step is to generate $D$ mesons with a
realistic spectrum. Generating these spectra from simulation would be a
difficult undertaking, so instead we chose to use experimental data collected
by the LEBC-EHS collaboration lebc_d-meson_1987 , at the CERN SPS running at
$400\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ with a hydrogen target. We then
randomly select a production and decay channel according to the relative
abundances of charmed mesons from alekhin_facility_2016  and the branching
fractions from bondarenko_phenomenology_2018 . Finally, we generate the
momenta of both the primary and secondary decay products at once. This is done
by first sampling all the momenta according to phase-space, independently for
each decay, and finally performing rejection sampling on these momenta using
the matrix element for the combined process. As a last step, we simulate the
geometrical acceptance by requiring the HNL to decay within the hidden sector
decay vessel, into two long-lived, charged particles which both intersect the
tracking station. In order to account for the (small) probability of the HNL
decaying inside the fiducial volume, each event is weighted by
$P_{\mathrm{decay}}(\tau)=\Gamma e^{-\Gamma\tau}$, where $\tau$ is the proper
time between the HNL production and decay. Throughout this paper, we assume
the particle identification to be perfectly efficient, which should be a
reasonably good approximation at SHiP hosseini_particle_2017 . The simulation
is described in details in appendix A.

###  4.2 LNC / LNV classification

For a given choice of relative squared mixing angles
$\left|\Theta_{\alpha}\right|^{2}$ (which are supposed to be known by the time
LNV is studied at SHiP), we generate a dataset for a range of HNL masses
between the $K$ and $D_{s}$ thresholds. For each HNL mass, we sample $3\cdot
10^{6}$ events with uniform weights, and keep only those passing the
acceptance cuts. The HNL is simulated as a single Majorana particle, which
ensures that the dataset contains equal numbers of LNC and LNV events, and is
also balanced with respect to the primary and secondary lepton charges.

Each event is labelled with a boolean flag set to false for LNC and true for
LNV, using the MC truth. The only observable quantities come from the HNL
decay in the vacuum vessel. They are: the momenta and charges of the lepton
$l_{\beta}^{\pm}$ and pion $\pi^{\mp}$, and the decay vertex $x_{D}$. Of these
quantities, we record a total of $19$ primary or derived features. Their
definitions can be found in table 1, and some typical distributions are
presented, as an example, in fig. 7, for both LNC and LNV processes. Finally,
from each dataset, we set aside $30\%$ of events for testing and $20\%$ for
validation, leaving us with $50\%$ of events for training the classifier.

Figure 7: Corner plot showing the correlations between five selected features, for a $1\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ HNL coupling to the muon. See table 1 for a description of the features. Each subplot shows, on the same scale, the marginal distributions of LNC and LNV events as a function of either one (on-diagonal plots) or two (off-diagonal plots) features. $1\mathrm{d}$ distributions are represented as histograms, and $2\mathrm{d}$ distributions as contour plots of the probability density.  Feature(s) | Description
---|---
Ql2 | Charge of the secondary lepton $l_{\beta}$
E1, p1x, p1y, p1z | Reconstructed HNL momentum $p_{N}=p_{l_{\beta}}+p_{\pi}$ (lab frame)
E2, p2x, p2y, p2z | Secondary lepton momentum $p_{l_{\beta}}$ (lab frame)
E3, p3x, p3y, p3z | Secondary pion momentum $p_{\pi}$ (lab frame)
pCMx, pCMy, pCMz | Secondary lepton momentum $p_{\mathrm{CM}}$ (HNL frame)
xD, yD, zD | Decay vertex (lab frame)
Table 1: The $19$ features recorded for each event.

For each dataset, we train a binary classifier to discriminate between LNC and
LNV decay chains. For this study, we use the LightGBM ke_lightgbm_2017
decision tree boosting algorithm, through the Python interface to the
reference implementation lightgbm . In order to perform simple classification,
we choose the binary objective. The training is discussed in more details in
section B.2. The accuracy of the trained classifier (as evaluated on the test
set) is presented in fig. 8 as a function of the HNL mass for three scenarios,
corresponding to an HNL coupling to electrons, muons, or equally to both.

Figure 8: Classification accuracy as a function of the mass, for an HNL
coupling to $e$, $\mu$, or equally to both.

###  4.3 Model selection

Assuming the true event distribution to match (or be sufficiently close to)
the simulated one, we can then use our trained classifier to classify each
event as either LNC or LNV. As stated in section 1, our main goal is to
distinguish the following two hypotheses:

  * •

$\mathcal{H}_{1}$: HNLs are Dirac or quasi-Dirac with $\delta M\tau\ll 1$ (LNC
decays only).

  * •

$\mathcal{H}_{2}$: HNLs are Majorana or quasi-Dirac with $\delta M\tau\gg 1$
(as many LNC / LNV decays).

Since the classifier is not perfectly accurate, its decision cannot be used to
directly confirm the presence of LNV processes, or constrain their existence.
If we knew the full distribution in feature space $\rho(z)$ for each
hypothesis, we could obtain an optimal test statistics by constructing the
corresponding likelihood ratio neyman_1933 . However, accurately estimating
$\rho(z)$ is a non-trivial task and would be error-prone, so we elected to use
a less powerful but more reliable, simplified model. Knowing the
classification accuracy $a$ for a given binary classifier, we compute the
likelihood of classifying $k$ events out of $N$ as LNV, and $N-k$ events as
LNC (independently of their specific feature vectors $z$) assuming that the
true fraction of LNV events is $f$. We then compute the best-fit value for $f$
and use Wilk’s theorem wilks_large-sample_1938  in order to determine whether
it significantly deviates from either $f=0$ (corresponding to
$\mathcal{H}_{1}$) or $f=\frac{1}{2}$ (corresponding to $\mathcal{H}_{2}$).

In order to estimate the “model-selection” sensitivity of SHiP, we then
compute, under each hypothesis and as a function of the HNL mass $M_{N}$ and
squared mixing angles $\left|\Theta_{\alpha}\right|^{2}$, the median
confidence level at which we can exclude the other hypothesis assuming $5$
years of nominal operation (i.e. $2\cdot 10^{20}$ protons on target). For each
true hypothesis, we finally draw the sensitivity limit by plotting, for each
$M_{N}$, the smallest $\left|\Theta_{\alpha}\right|^{2}$ for which this median
confidence level is at least $0.9$. In other words, for mixing angles above
this limit, SHiP has a probability of at least $1/2$ of disfavouring one
hypothesis at $\mathrm{CL}=0.9$ if the other is realized. The construction of
these confidence limits is described in details in section B.3, and the
resulting sensitivity plots are presented in section 5.1.

###  4.4 Resolving HNL oscillations

So far we have only considered the two extreme cases ($\mathcal{H}_{1}$ and
$\mathcal{H}_{2}$), where the HNL(s) behave either as a single Dirac or
Majorana particle. However, as discussed in section 2.2, if two nearly
degenerate HNLs form a quasi-Dirac pair, both LNC and LNV decay chains will be
present, with a non-trivial ratio $\neq 0,1$, and the corresponding decay
rates will feature oscillations as a function of the proper time $\tau$
between the HNL production and decay events, with the characteristic
$1\pm\cos(\delta M\tau)$ dependence described by eq. 10, where ($+$)
corresponds to LNC and ($-$) to LNV.

For $\delta M\sim 10^{-6}\;$\mathrm{e}\mathrm{V}$$, $\delta M\tau$ will be of
order $2\pi$ at SHiP, leading to potentially resolvable oscillations, provided
we can accurately reconstruct the proper time $\tau$ between the HNL
production and decay. Expressing it as $\tau=\nicefrac{{L}}{{\beta\gamma}}$,
we see that this can be accomplished if we have sufficiently accurate
vertexing and energy reconstruction. At SHiP, the precision on $L$ will be
limited by the impossibility of reconstructing the primary vertex within the
target. The energy resolution, despite being sufficient for particle
identification, is not enough for reconstructing $\tau$ (see sections 4.7 and
4.10 in ship_collaboration_facility_2015 ). However, the momentum resolution,
combined with the dispersion relation (assuming the HNL mass to be known
already with sufficient accuracy) should allow reconstructing $\gamma$ much
more precisely. The high vertexing and momentum resolution permitted by the
SHiP tracker, together with our method for (statistically) distinguishing LNC
from LNV processes (described in section 4.3), should therefore make it
possible to resolve the oscillation pattern in part of the parameter space.

In order to search for HNL oscillations, we first classify the observed events
using a model trained (for _one_ HNL) at the corresponding mass. We thus
assume again that we have sufficiently many events that the HNL mass $M_{N}$
is well known. The events are then binned in proper time $\tau$, which is the
relevant variable for oscillations of massive, relativistic particles. Instead
of using the predicted class, here we implement the classifier decision as a
weight for the binned events, using the predicted probability
$p_{\mathrm{LNV}}$. This weight contains more information than the class does,
since it acts as a measure of uncertainty by taking values close to $1/2$ for
ambiguous events, and closer to $0$ or $1$ for unambiguous ones. However,
without applying further corrections, the sum of these probabilities would
average to $N\langle p_{\mathrm{LNV}}\rangle$ for the entire sample of $N$
events. If used directly as weights, they would therefore cause the
oscillatory pattern to be hidden among Poisson fluctuations. In order to
reveal this pattern, we instead weight the events by
$p_{\mathrm{LNV}}-\overline{p_{\mathrm{LNV}}}$, where
$\overline{p_{\mathrm{LNV}}}$ is the sample average of the estimated
$p_{\mathrm{LNV}}$. This weight averages to zero over the entire sample, which
limits the impact of Poisson fluctuations.

HNL oscillations are implemented in our simulation by first generating events
without taking interference into account then, in a second time, performing
rejection sampling based on the proper time $\tau$, following eq. 10. The
results obtained using this simulated data set are presented in section 5.2.

##  5 Results

###  5.1 Sensitivity to Lepton Number Violation

In order to easily compare our results to existing exclusion bounds or to the
sensitivities of future experiments, let us consider two simplified models
where a single HNL exclusively mixes with the electron or muon
neutrino.999Within the seesaw mechanism, it is impossible to generate the two
observed light neutrino mass differences with a single HNL, or if HNLs mix
with one generation only drewes_minimal_2019 . The two benchmarks presented in
figs. 9(a) and 9(b) are thus simplifications, used here because they are
consistent with the parametrization employed by the PBC working group.  As can
be seen in fig. 8, more generic mixing patterns with the $e$ and $\mu$ flavors
do not significantly degrade the classification accuracy; therefore they
should leave the limits presented below mostly unchanged. However, if a
significant fraction of HNLs is produced through mixing with the $\tau$
neutrino, then the present analysis would need to be modified to handle
secondary production of HNLs in $\tau$ decays, including spin correlation
effects.

As discussed in section 4.3, we define the sensitivity to lepton number
violation as the smallest mixing angles for which SHiP has a $1/2$ probability
of either rejecting or detecting LNV, if it is respectively absent or present
with the same rate as LNC. The results are presented in footnote 11, along
with various existing exclusion bounds and detection sensitivity101010The
usual sensitivity, by opposition to the sensitivity to lepton number violation
discussed here. limits for planned or proposed experiments, extracted from the
report of the Physics Beyond Colliders working group pbc_2019 . We only show
the sensitivities of experiments which can not only set exclusion bounds, but
also reconstruct the HNL mass, should it be observed. Note that in order to be
consistent with the SHiP detection sensitivity, which was computed for one
Majorana HNL, we present our results for one HNL as well. In the realistic
case of $\mathcal{N}\geq 2$ HNLs, both curves must be scaled down by a factor
of $\mathcal{N}^{1/2}$. Above the black dashed line, SHiP should be able to
distinguish Dirac-like ($\mathcal{H}_{1}$) and Majorana-like
($\mathcal{H}_{2}$) HNLs. We have discarded the HNL masses for which the early
stopping criterion returned the first iteration as the best, since it suggests
that the classifier has failed to learn anything about the data. Below
$0.7\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$, additional production channels
$H\to h_{V}^{\prime}l_{\alpha}N$ (where $h_{V}^{\prime}$ denotes a vector
meson) become significant, and have not been implemented with spin
correlations in our Monte-Carlo simulation. Therefore we also restrict the HNL
mass to $M_{N}\gtrsim$0.7\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$.
Additionally, since the sensitivity is almost identical for excluding
$\mathcal{H}_{1}$ or $\mathcal{H}_{2}$, we only plot one limit, which
corresponds to excluding $\mathcal{H}_{1}$ at $90\%~{}\mathrm{CL}$ if LNV is
actually present.

We can see that the larger number of accepted events at higher masses
initially compensates for the worse classification accuracy, but is not
sufficient any more as we approach the $D$ threshold. In practice, we expect
that systematic uncertainties about the $D$ spectrum and the simulation will
decrease the sensitivity at both ends of the mass range, where the
classification accuracy is already close to $1/2$. Comparing the results to
the SHiP detection sensitivity, we see that around
$1\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$, the model-selection sensitivity
limit is about one order of magnitude above the detection one, while remaining
well below the planned NA62++ limit as well as existing bounds.

This leads us to an interesting conclusion: there exists a non-trivial region
of parameter space, unconstrained by current or near-future experiments, where
SHiP would not only be able to detect HNLs, but also characterize them as
either Dirac-like or Majorana-like particles. As discussed in sections A.3 and
B.4, this conclusion is robust with respect to uncertainties on the heavy
meson spectrum.

(a) HNL mixing with $\nu_{e}$.

(b) HNL mixing with $\nu_{\mu}$.

Figure 9:  SHiP sensitivity to lepton number violation. The dashed curve is
the “model-selection” sensitivity computed in this work. Dotted curves are the
detection sensitivities for the proposed or planned experiments which can
reconstruct the HNL mass. Coloured, filled areas are regions of parameter
space which have been excluded by previous experiments. The grey filled area
denoted by BBN indicates the region which is incompatible with Big Bang
Nucleosynthesis. Below the seesaw limit111111The seesaw limit can only be
rigorously computed if the mixing angles are consistent with the seesaw
equation (7). This is not possible for HNLs mixing with only one generation,
nor for a single HNL. The limits presented here instead correspond to the
“naive” estimate $\sum m_{\nu}\leq
M_{N}\cdot\sum_{\alpha}|\Theta_{\alpha}|^{2}$, where we have assumed the
lightest neutrino to be massless. (hatched region), mixing angles are too
small to produce the observed neutrino masses.

###  5.2 Resolvable quasi-Dirac oscillations

Figure 10: Events binned by proper time $\tau$ and weighted by
$p_{\mathrm{LNV}}-\overline{p_{\mathrm{LNV}}}$, revealing the oscillatory
pattern, for $M_{N}=$1\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$,
$\left|\Theta_{\mu}\right|^{2}=2\cdot 10^{-8}$, $\Theta_{e}=\Theta_{\tau}=0$
and $\delta M=$4\text{\cdot}{10}^{-7}\text{\,}\mathrm{e}\mathrm{V}$$.

The result of the procedure described in section 4.4 is presented in fig. 10
for a new simulated dataset (independent from the training set), corresponding
to a quasi-Dirac pair of mass $1\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$, mass
splitting $4\text{\cdot}{10}^{-7}\text{\,}\mathrm{e}\mathrm{V}$, and mixing
with muon neutrinos only, with a squared mixing angle
$\left|\Theta_{\mu}\right|^{2}=2\cdot 10^{-8}$. The oscillatory pattern is
manifest at $\tau<$5\text{\,}\mathrm{m}$$, where most of the events fall. At
larger $\tau$ it is hidden in Poisson fluctuations. The uncertainty on $\tau$
at SHiP is dominated by the (boosted) length of the target
$\sim$0.1\text{\,}\mathrm{m}$$, which contains the unresolved primary vertex.
It could smear out fast oscillations, in which case an accurate treatment of
this uncertainty is needed in the simulation. However, for longer oscillation
periods like the one shown in fig. 10, its effect should be negligible.
Deriving precise sensitivity limits for HNL oscillations is beyond the scope
of this paper, since it is likely that no simple analytical expression exists
for them, due to the more complex test statistics required, compared to the
detection or model-selection limits. HNL oscillations might for instance be
amenable to methods such as maximum likelihood estimation, wavelets, or
matched filtering, for which the null distribution can be estimated
numerically using a (computationally expensive) bootstrapping procedure.

##  6 Conclusions

The SHiP experiment is set to have an unprecedented detection reach for a
variety of models containing feebly interacting particles, such as Heavy
Neutral Leptons (HNLs). A distinctive feature of SHiP among other intensity
frontier experiments is its decay spectrometer, which allows it to not only
place exclusion bounds, but also perform event reconstruction and measure the
HNL properties. The simplest consistent HNL model accessible at SHiP contains
two nearly-degenerate HNLs, which can undergo oscillations. Their mass
splitting $\delta M$ is of particular interest, since it greatly influences
their phenomenology as well as early-Universe cosmology (specifically,
baryogenesis and dark matter production).

In the present work, we have investigated to which extent SHiP may be able to
constrain or even measure $\delta M$. Depending on the scale of the
oscillation phase $\delta M\tau$ accessible at an experiment, HNLs may or may
not exhibit lepton number violation (LNV). The problem thus amounts to
distinguishing LNC from LNV decay chains (fig. 2) in a beam-dump setting (fig.
1), where the primary lepton cannot be observed. We have shown that the
angular distribution of the _visible_ secondary decay products provides a
partial solution to this problem, since, depending on the HNL mass, it can
significantly differ between LNC and LNV in the laboratory frame (fig. 7).
This result has been qualitatively understood in the simplified case of two-
body decays in the massless limit (figs. 3 and 6). In order to handle more
realistic cases, a Monte-Carlo simulation has been employed to generate
accurate data sets of LNC and LNV events, including spin correlations and
geometrical acceptance. The different distributions of the kinematic variables
thus allow discriminating between LNC and LNV events using multivariate
analysis; and with sufficiently many events, it becomes possible to
statistically detect or exclude lepton number violation.

In order to produce sufficiently accurate training sets, our simulation must
satisfy several requirements. It should be able to generate all the relevant
two- and three-body meson decays containing an HNL (fig. 5), as well as the
selected HNL decay channel $N\to\pi^{\mp}l_{\beta}^{\pm}$. It should be
accurate for $\mathrm{G}\mathrm{e}\mathrm{V}$-scale HNLs, and should account
for the spin correlations between the primary and secondary decays. Finally,
it should run sufficiently fast to allow producing large training sets for
various hypotheses and parameters. In order to meet all these requirements, we
have written our own Monte-Carlo simulation, the output of which is used to
train a binary classifier.

Knowing the accuracy of the classifier decision (fig. 8) for a given mass and
(relative) mixing angles, we can finally draw a “model-selection” sensitivity
limit in the $(M_{N},\left|\Theta\right|^{2})$ plane (shown in figs. 9(a) and
9(b)), above which SHiP should be able to either discover or rule out lepton
number violation from HNLs. Interestingly, this limit lies below the detection
sensitivity of near-future experiments such as NA62++. This leads to a
striking conclusion: SHiP might be able to not only _discover_ HNLs, but also
_characterize_ them as either “Dirac-like” or “Majorana-like” fermions
(depending on whether they feature LNV) even if previous experiments _see no
signal at all_. Better yet, if the mass splitting between the two HNLs is of
order $\delta M\sim 10^{-6}\;$\mathrm{e}\mathrm{V}$$, SHiP should be able to
_resolve the oscillations_ of HNLs (fig. 10), given sufficiently many events.
Intriguingly, this mass splitting falls within the range required for
producing dark matter in the $\mathrm{\nu MSM}$ Canetti:2012kh . Its
measurement—or constraining—would therefore be an important test of
cosmological models.

###### Acknowledgements.

The authors would like to thank Mikhail Shaposhnikov for his suggestion of
investigating HNL oscillations at SHiP, Oleg Ruchayskiy for his helpful
comments and proofreading of the present manuscript, Sonia Bouchiba and
Federico Leo Redi for stimulating discussions about helicity effects in HNL
decays, and Kyrylo Bondarenko and Elena Graverini for sharing with us their
code to compute HNL decay widths. We thank the SHiP collaboration for their
interest in our work. We are grateful to Annarita Buonaura, Richard Jacobsson,
and Nicola Serra for their helpful comments and clarifications about the
dimensions of the facility. JLT acknowledges support from the Carlsberg
foundation. The work of IT was supported in part by the ERC-AdG-2015 grant
694896.

##  Appendix A Simulation

###  A.1 Overview

It is not obvious whether the different angular correlations of LNC and LNV
events lead to an observable effect in a realistic beam-dump experiment. To
answer this question, we have devised a toy Monte-Carlo simulation, inspired
from the one used in ship_collaboration_sensitivity_2018 , to simulate the
production and decay of HNLs at the SHiP experiment
ship_collaboration_facility_2015 ; alekhin_facility_2016  (represented on fig.
1).

The simulation of rare BSM processes with spin correlations entails two main
requirements. First, we cannot afford to simulate all the possible processes,
since, due to the small HNL mixing angles, the decay chains mediated by an HNL
only represent a tiny fraction of all decays. Instead, we only simulate the
BSM processes, and use importance sampling (i.e. introduce weights) in order
to obtain the correct absolute number of events and expectation values
(section A.2).

Secondly, we cannot sample the primary and secondary decays separately, since
they are not independent. Instead, we construct all possible decay chains for
the production and decay processes of interest, and sample the entire chain at
once, with a probability proportional to its combined branching fraction. The
momenta of all the decay products are then sampled simultaneously, using the
matrix element for the entire chain (section A.4).

In addition, in order to accurately model the SHiP experiment, we need to
sample the heavy meson momenta from a realistic spectrum (section A.3) and
take into account the finite size of SHiP and its geometrical acceptance
(section A.5). Finally, since most machine learning algorithms take unweighted
data points as input, it is necessary to perform a last step of rejection
sampling in order to produce a training set consisting of events with equal
weights (section A.6).

###  A.2 Decay chains

As discussed in bondarenko_phenomenology_2018 , the dominant HNL production
process at SHiP is from weak decays of the lightest charmed or beauty mesons.
In the present study, we focus on HNL masses below the $D_{s}$ mass, and only
select the fully reconstructible secondary decays
$N\rightarrow\pi^{\pm}l_{\beta}^{\mp}$, By producing long-lived, charged
particles which can be measured by the decay spectrometer located at the end
of the decay vessel, they allow the HNL momentum to be reconstructed. The
efficiency of particle identification at SHiP is high enough
hosseini_particle_2017  that we can approximate it as one for the present
estimate. Therefore we do not need to simulate decay chains containing any
other secondary decays.

For the mixing angles of interest (i.e. below existing bounds), the fraction
of all decays which are mediated by an HNL is tiny. We therefore need to use
importance sampling in order to efficiently simulate only the processes of
interest. For every proton on target (POT), the probability of producing a
heavy hadron of species $H$ is:

$P(H)=\frac{\sigma_{cc}}{\sigma_{pN}}\cdot A_{H}$ (19)

where $\sigma_{cc}$ is the production cross-section for charmed hadrons,
$\sigma_{pN}$ the interaction cross-section for protons hitting the target
nuclei, and $A_{H}$ is the relative abundance of the charmed hadron species
$H$ (as given in appendix A of alekhin_facility_2016 ). The _nominal_ (i.e.
physical) probability of producing an HNL which mediates a given decay chain
$H\rightarrow[h^{\prime}]l_{\alpha}(N\rightarrow l_{\beta}h^{\prime\prime})$
(irrespective of whether the decay is observed in the detector) is then:

$P\left(H\rightarrow[h^{\prime}]l_{\alpha}(N\rightarrow
l_{\beta}h^{\prime\prime})\right)=P(H)\cdot P(h^{\prime}l_{\alpha}N|H)\cdot
P(l_{\beta}h^{\prime\prime}|h^{\prime}l_{\alpha}N)\\\
=\frac{\sigma_{cc}}{\sigma_{pN}}\cdot
A_{H}\cdot\mathrm{Br}_{\mathrm{prod}}(H\rightarrow[h^{\prime}]l_{\alpha}N)\cdot\mathrm{Br}_{\mathrm{decay}}(N\rightarrow
l_{\beta}h^{\prime\prime})$ (20)

where the last two terms are the production and decay branching ratios for
HNLs in the considered decay chain. The _importance_ distribution $P^{\prime}$
is defined as a uniform scaling for decay chains involving an HNL, and as zero
for all other outcomes:

$\begin{cases}P^{\prime}\left(H\rightarrow[h^{\prime}]l_{\alpha}(N\rightarrow
l_{\beta}h^{\prime\prime})\right)=\frac{1}{w_{\mathrm{prod}}}P\left(H\rightarrow[h^{\prime}]l_{\alpha}(N\rightarrow
l_{\beta}h^{\prime\prime})\right)\\\ P^{\prime}(\text{no HNL})=0\end{cases}$
(21)

where $w_{\mathrm{prod}}$ is the weight to be applied to all the chains
sampled from the importance distribution, and corresponds to the total
probability of producing an HNL according to the nominal distribution:

$w_{\mathrm{prod}}=\sum_{\mathrm{chains}}P\left(H\rightarrow[h^{\prime}]l_{\alpha}(N\rightarrow
l_{\beta}h^{\prime\prime})\right)$ (22)

When computing expected numbers of events over the entire duration of the SHiP
experiment, which represents an integrated $N_{\mathrm{POT}}=2\cdot 10^{20}$
protons on target for $5$ years of nominal operation, we must further multiply
by $N_{\mathrm{POT}}$ the expectation values obtained for one event. This is
most easily done by simply multiplying the total weights by
$N_{\mathrm{POT}}$.

###  A.3 Heavy meson spectrum

Once a chain is selected, we sample the momentum of the corresponding charmed
meson from the spectrum measured by the LEBC-EHS collaboration lebc_d-
meson_1987  at the CERN SPS running at
$400\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ with a hydrogen target. The
differential cross-section is parametrized as the product of a $\beta$
distribution in $x_{F}$ and an exponential distribution in $p_{T}^{2}$:

$\frac{\mathrm{d}^{2}\sigma}{\mathrm{d}x_{F}\mathrm{d}p_{T}^{2}}=\sigma\frac{(n+1)b}{2}(1-\left|x_{F}\right|)^{n}e^{-bp_{T}^{2}}$
(23)

with the best-fit values $n=4.9\pm 0.5$ and $b=(1.0\pm
0.1)\;$\mathrm{G}\mathrm{e}\mathrm{V}^{-2}$$. We thus implicitly assume the
spectrum to be separable. Due to their very similar mass, and to compensate
for the lack of data, we assume $D_{s}$ mesons to share the same spectrum as
$D$ mesons.

By using the spectrum for a hydrogen target, we effectively neglect cascade
production of heavy hadrons inside the target, leading us to underestimate the
number of hadrons produced at the low-energy end of the spectrum. This could
be problematic if their $p_{T}$ spectrum happens to be significantly different
from that of primary hadrons produced in $pp$ collisions. However, the lower
acceptance for these softer hadrons should help mitigate the issue. In fig.
11, we show how varying the width of the heavy meson $p_{T}$ spectrum affects
the final sensitivity. As expected, a larger $p_{T}$ spread reduces the
sensitivity, while a narrower spectrum improves it.

Figure 11: Effect of varying the width of the heavy meson $p_{T}$ spectrum on
the sensitivity to lepton number violation ($90\%$ CL), for an HNL coupling to
the muon. Black lines represent the model-selection sensitivity of SHiP for
various values of $\langle p_{T}^{2}\rangle$. The dashed line corresponds to
the best-fit value $\langle
p_{T}^{2}\rangle=$1\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}^{2}$$ from the
LEBC-EHS collaboration lebc_d-meson_1987 .

###  A.4 Decay product momenta

In order to preserve spin correlations between the HNL siblings and its decay
products, we simulate both the HNL production and decay processes at once. For
the masses and couplings of interest, the HNL is long-lived and can be assumed
to be on its mass shell. Therefore the phase-space sampling can be performed
independently for the primary and secondary decays. We use the $m$-generator
algorithm james_monte_1968  for that, as described in ilten_electroweak_2014 .
In order to sample events with a probability proportional to the squared
transition amplitude, we then perform rejection sampling, taking the phase-
space distribution as proposal distribution, and an acceptance probability
proportional to the spin-summed, squared matrix elements (12) and (13) for the
entire decay chain. Only the spin states of the external particles (which
interact with the detector and are thus “measured” in the quantum mechanical
sense) are summed over.

###  A.5 Geometry

In order to model the geometry of the SHiP experiment, we must account for the
finite size of the detector and its geometrical acceptance. In the current
SHiP design (represented on fig. 1), the fiducial volume consists of an
evacuated right pyramidal frustum of length $50\text{\,}\mathrm{m}$, located
at a distance of $50\text{\,}\mathrm{m}$ from the target, and with horizontal
and vertical sides $5\text{\,}\mathrm{m}$ and $10\text{\,}\mathrm{m}$
respectively at the far end. It is followed by a $10\text{\,}\mathrm{m}$ long
tracking station.

To estimate the probability of the HNL decaying within the fiducial volume and
passing the acceptance cuts, we use once again importance sampling for
sampling the decay vertex. This is required in order to overcome the
potentially very long HNL lifetime, which could cause most of them to decay
away from the experiment. We choose an importance distribution (approximately)
covering the fiducial volume, by sampling the decay vertex uniformly along the
HNL momentum, at a distance such that it falls inside the decay vessel. The
nominal decay probability density is, as a function of the proper time $\tau$
(or boost factor $\gamma$ and distance $L$) between the HNL production and
decay:

$P_{\mathrm{decay}}(\tau)=\Gamma e^{-\Gamma\tau}\quad\Longrightarrow\quad
P_{\mathrm{decay}}(L|\gamma)=\frac{\Gamma}{\beta\gamma}e^{-\frac{\Gamma
L}{\beta\gamma}}$ (24)

The partial weight resulting from this importance sampling step is therefore:

$w_{\mathrm{decay}}(L|\gamma)=\frac{\Gamma
L_{\mathrm{DV}}}{\beta\gamma\cos(\theta)}e^{-\frac{\Gamma L}{\beta\gamma}}$
(25)

where $L_{\mathrm{DV}}=$50\text{\,}\mathrm{m}$$ is the length of the decay
vessel and $\theta$ the angle between the HNL momentum and the beam axis. In
the linear regime, where $\Gamma\tau\ll 1$, this partial weight reduces to
$w_{\mathrm{decay}}(L|\gamma)\cong\nicefrac{{\Gamma
L_{\mathrm{DV}}}}{{\beta\gamma\cos(\theta)}}$.

We finally apply acceptance cuts by requiring the HNL to decay within the
decay vessel, and the trajectories of its two decay products
($l_{\beta}^{\mp}$ and $\pi^{\pm}$) to intersect the tracking station located
at its far end.

###  A.6 Unweighting

As a last step, we perform again rejection sampling on the weighted events in
order to obtain a set of events with equal weights, which are easier to
analyse and process with machine learning algorithms. This is done by
accepting events with a probability proportional to their weight, and can be
justified as follows.

Let $X$ denote a random variable representing the simulated event, and $x$ a
concrete realization of it. Let $f(x)=P(X=x)$ be the nominal (i.e. true)
distribution and $g(x)$ the importance distribution, such that $g(x)\neq 0$
for all outcomes $x$ in the domain of interest $\Omega$ (i.e. all relevant
observables must have their support in $\Omega$). If $x$ is sampled from the
importance distribution $g(x)$, its associated weight will be
$w(x)=f(x)/g(x)$. Let $M$ be an upper bound on $w(x)$, i.e. $M\geq
w(x),\forall x\in\Omega$. If we choose the acceptance probability to be
$a(x)\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny def}}}{=}}w(x)/M\leq 1$,
then it immediately follows that the accepted events, effectively drawn from
the new importance distribution $g(x)\cdot a(x)$, will have uniform weight
$M$.

It is therefore possible to perform rejection sampling a posteriori in order
to produce uniformly weighted events. However, storing all the generated
events, many of which will eventually be rejected, would be inefficient from a
memory perspective. A more economical solution, which we decided to use,
consists in performing rejection sampling directly as events are being
generated. This requires estimating an upper bound $M$ on the weights, during
an initial burn-in phase.

##  Appendix B LNC/LNV classification

At leading order in the light lepton and hadron masses, the matrix elements
for LNC and LNV decay chains have a straightforward analytical dependence on
the invariant mass $s_{ll}$ of the charged lepton pair. However, unlike in
collider experiments, this variable is not readily available in a beam-dump
setting, due to the primary lepton being unobservable. As we saw in section
3.2, the different angular correlations between the charged leptons can
nevertheless lead to residual correlations between the visible HNL decay
products. The absence of an obvious test statistics, along with the almost
background-free conditions and highly efficient PID at SHiP
hosseini_particle_2017 , makes the task of distinguishing LNC from LNV ideally
suited for multivariate analysis. In the following subsections, we describe
how we generate the training set (section B.1), the classifier used to
discriminate between LNC and LNV events (section B.2), how to produce a
sensitivity limit from its output (section B.3), and finally how sensitive is
the classification to systematic uncertainties on the heavy meson spectrum
(section B.4).

###  B.1 Dataset

As mentioned in section 4.2, we need to generate datasets for various HNL
masses $M_{N}$ and rays in $\left|\Theta_{\alpha}\right|^{2}$ space, where
$\alpha=e,\mu$ (the overall normalization does not matter). In practice, we
choose a mass range spanning the region between the $K$ and $D_{s}$
thresholds, and consider several benchmark models with fixed
$\left|\Theta_{e}\right|^{2}:\left|\Theta_{\mu}\right|^{2}$ ratios.121212We do
not consider HNL production through $\tau$ mixing in this work, since it would
have required to implement secondary production from $\tau$ decays. It is
negligible in the considered mass range unless the $\Theta_{\tau}$ mixing
angle is significantly larger than the others, as can be seen in fig. 5. In
addition, visible HNL decays through $\tau$ mixing are forbidden below the
$\tau$ threshold.  For each choice of physical parameters, we sample $3\cdot
10^{6}$ events with uniform weights. This is done by sampling sufficiently
many weighted events and, as they are being generated, “unweighting” them by
performing rejection sampling with an acceptance probability proportional to
their weight. Only events which pass the acceptance cuts are used for
training. In the simulation, the HNL is taken to be a single Majorana
particle, such that the dataset contains equal numbers of LNC and LNV events
and is balanced with respect to the primary and secondary lepton charges. We
select only the fully reconstructible HNL decays
$N\rightarrow\pi^{\mp}l_{\beta}^{\pm}$, which do not contain an unobservable
light neutrino, and produce long-lived charged particles which can be measured
by the decay spectrometer. For the sake of simplicity, we will assume the PID
to be perfectly efficient throughout this analysis. Non-trivial efficiencies
are expected to slightly reduce the final sensitivity reach. As explained in
section 4.2, each event is labelled as being either LNC or LNV, and we record
the $19$ observable features listed in table 1. The dataset is split into
training / validation / test sets with respective proportions $0.5:0.2:0.3$.

###  B.2 Classifier

We employ the LightGBM ke_lightgbm_2017  gradient boosting algorithm, accessed
through the Python interface to the reference implementation lightgbm . For
classification, we choose the binary objective. We use early stopping based on
the binary log-loss (binary_logloss) and the area-under-curve (auc) metrics,
with a $10$ round threshold. The hyperparameters num_leaves and learning_rate
are manually optimized by maximizing the above two metrics on the validation
set. The classification accuracy is presented in fig. 8 as a function of the
HNL mass $M_{N}$ for two orthogonal scenarios, corresponding to the HNL
coupling exclusively to electrons
($\left|\Theta_{e}\right|^{2}:\left|\Theta_{\mu}\right|^{2}:\left|\Theta_{\tau}\right|^{2}=1:0:0$)
or muons
($\left|\Theta_{e}\right|^{2}:\left|\Theta_{\mu}\right|^{2}:\left|\Theta_{\tau}\right|^{2}=0:1:0$),
and a third one where it couples equally to both
($\left|\Theta_{e}\right|^{2}:\left|\Theta_{\mu}\right|^{2}:\left|\Theta_{\tau}\right|^{2}=1:1:0$).

It is instructive to understand the origin of this dependence, if only to make
sure that it corresponds to a physical effect. LightGBM provides a way to
estimate the feature importance, by counting the number of times a feature is
used to split a tree. Those are listed in table 2 for a
$1\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ HNL coupling to muons (which
results in a classification accuracy of $63.5\%$). They reveal that the most
important features are the transverse components of the momenta of the HNL
decay products. Indeed, it is possible to successfully train a model using a
single feature such as the transverse momentum $p_{T,\mu}$ of the secondary
muon, while still obtaining a classification accuracy of $61.5\%$ for the same
dataset.

Inspecting the results more closely (see fig. 7) shows that LNV events have on
average a slightly larger transverse momentum than LNC ones. This is
consistent with our discussion from section 3.2, and allows us to understand
the mass dependence. At large HNL masses, as we approach the closing mass of
$D$ meson leptonic decays, the kinetic energy of the HNL in the heavy meson
rest frame decreases, until it becomes so small that the difference between
LNC and LNV becomes negligible compared to the transverse momentum spread of
the heavy meson spectrum. As the HNL mass decreases, 3-body semileptonic decay
channels open, and become dominant at lower masses. The additional meson takes
away part of the energy from the HNL, leaving it with insufficient kinetic
energy to “escape” the transverse momentum spread of the heavy meson spectrum.
Finally, the large boost of the heavy mesons along the beam axis washes out
most of the information contained in the longitudinal part of all laboratory
frame momenta, which explains their low importance.

Feature | p2y | p3y | p2x | p3x | pCMz | zD | xD | yD | p1x | pCMy
---|---|---|---|---|---|---|---|---|---|---
# splits | 302 | 282 | 243 | 238 | 141 | 114 | 105 | 97 | 91 | 85

Feature | pCMx | p1y | E1 | E2 | E3 | p3z | p2z | p1z | Ql2
---|---|---|---|---|---|---|---|---|---
# splits | 77 | 74 | 69 | 67 | 61 | 53 | 34 | 14 | 9

Table 2: Feature importance for a $1\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$
HNL coupling to $\mu$.

###  B.3 Sensitivity to lepton number violation

As stated in section 4.3, our main goal is to distinguish between the
following two hypotheses using exclusively the classifier decision (i.e. not
the underlying feature vector $z$):

  * •

$\mathcal{H}_{1}$: HNLs are Dirac or quasi-Dirac with $\delta M\tau\ll 1$ (LNC
decays only).

  * •

$\mathcal{H}_{2}$: HNLs are Majorana or quasi-Dirac with $\delta M\tau\gg 1$
(LNC and LNV decays).

Those can be expressed as special cases of a more general hypothesis
$\mathcal{H}(f)$, $f\in[0,1]$, parametrized by the relative frequency $f$ of
LNV events:

  * •

$\mathcal{H}(f)$: $(\text{LNV rate})=f\times(\text{total rate})$.

such that $\mathcal{H}_{1}=\mathcal{H}(f=0)$ and
$\mathcal{H}_{2}=\mathcal{H}(f=1/2)$.

We model the classifier decisions using a $2\times 2$ confusion matrix
$C_{ij}=P(i\text{ classified as }j)$, where $i,j=1,2$ correspond to the two
classes, respectively LNC and LNV. The confusion matrix can be expressed in
terms of the classification accuracies as:

$C=\begin{pmatrix}a_{1}&1-a_{1}\\\ 1-a_{2}&a_{2}\end{pmatrix}$ (26)

Suppose we observe $N$ events passing the selection cuts, $k$ of which are
classified as LNV. Then, under $\mathcal{H}(f)$, the likelihood of classifying
$N-k$ events in class 1 (LNC) and $k$ in class 2 (LNV) is given by the
following binomial distribution:

$\mathcal{L}(k;f)={N\choose
k}\big{(}a_{2}f+(1-a_{1})(1-f)\big{)}^{k}\big{(}a_{1}(1-f)+(1-a_{2})f\big{)}^{N-k}$
(27)

Under hypothesis $\mathcal{H}_{1}$, i.e. all events are LNC, this likelihood
reduces to:

$\mathcal{L}_{1}(k)=\mathcal{L}(k;f=0)={N\choose k}(1-a_{1})^{k}a_{1}^{N-k}$
(28)

while under hypothesis $\mathcal{H}_{2}$, i.e. events come from either class
with equal probability, it becomes:

$\mathcal{L}_{2}(k)=\mathcal{L}(k;f=1/2)={N\choose
k}\frac{(1+a_{2}-a_{1})^{k}(1+a_{1}-a_{2})^{N-k}}{2^{N}}$ (29)

For many models, including LightGBM (with a balanced training set),
$a_{1}\approx a_{2}\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny def}}}{=}}a$.
In this limit, $\mathcal{L}_{2}(k)$ simplifies to ${N\choose k}2^{-N}$.

Since $\mathcal{H}_{1,2}$ and $\mathcal{H}(f)$ are nested, then, assuming we
have sufficiently many events, we can use Wilk’s theorem131313A potential
issue in the case of $\mathcal{H}_{1}$ could be that the null value $f=0$ lies
on the boundary of the domain $[0,1]$ of $f$, while Wilk’s theorem requires
the true value to be in the interior of the parameter space. However,
$\ln(\mathcal{L}(k;f))$ has a well-behaved analytical continuation over a
domain larger than $[0,1]$. As long as the estimator $\hat{f}$ has a
sufficiently small variance, this boundary effect can therefore be ignored and
Wilk’s theorem still applies. See algeri_searching_2019  for a comprehensive
discussion of the validity conditions of Wilk’s theorem.  to try to exclude
$\mathcal{H}_{1,2}$ . To this end, we construct the two likelihood ratios
$\Lambda_{1,2}(k)$ as:

$\Lambda_{i}(k)=\frac{\mathcal{L}_{i}(k)}{\mathcal{L}(k;\hat{f})},\quad i=1,2$
(30)

where $\hat{f}$ is the maximum likelihood estimator for $f$:

$\hat{f}=\frac{1-a-k/N}{1-2a}$ (31)

Wilk’s theorem states that if $\mathcal{H}_{i}$ ($i=1\text{ or }2$) is
realized, then $-2\ln(\Lambda_{i}(k))$ follows a $\chi^{2}$ distribution with
one degree of freedom. Conversely, if we observe $-2\ln(\Lambda_{i}(k))>2.7$,
then $\mathcal{H}_{i}$ will be disfavoured at $90\%$ CL. If both hypotheses
$\mathcal{H}_{1,2}$ were disfavoured simultaneously, this would suggest
$\delta M\tau\sim 2\pi$ and potentially resolvable HNL oscillations.

If hypothesis $\mathcal{H}_{1}$ is actually realized, we expect $k$ to take a
value around the expected number of events misclassified as LNV: $(1-a)N$,
which, for large $N$, is approximately equal to the median. The median of the
log-likelihood-ratio when testing for $\mathcal{H}_{2}$ is therefore:

$\mathrm{med}_{1}\left(\ln(\Lambda_{2})\right)\approx-N\underbrace{\big{(}\ln(2)+a\ln(a)+(1-a)\ln(1-a)\big{)}}_{\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny
def}}}{=}}\;l_{1}(a)>0}$ (32)

If, instead, $\mathcal{H}_{2}$ is realized, then we expect $k$ to take a
median value of approximately $N/2$, such that:

$\mathrm{med}_{2}\left(\ln(\Lambda_{1})\right)\approx
N\underbrace{\left(\ln(2)+\frac{1}{2}\ln(a)+\frac{1}{2}\ln(1-a)\right)}_{\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny
def}}}{=}}\;l_{2}(a)<0}$ (33)

For a fixed confidence level, we can invert these two formulas to estimate,
for each true hypothesis $\mathcal{H}_{i}$, $i=1,2$, the median number of
events $N_{i}(a)$ required to exclude the other hypothesis:

$N_{i}(a)=\left|\frac{\ln(\Lambda_{\mathrm{cr}})}{l_{i}(a)}\right|$ (34)

with $-2\ln(\Lambda_{\mathrm{cr}})\approx 2.7$ for a $90\%~{}\mathrm{CL}$. The
higher the classification accuracy, the less events are required to reach the
target, while accuracies close to $1/2$ do not allow distinguishing the two
hypotheses, as $N_{i}(1/2)\rightarrow\infty$. So far we have only considered
the two extreme cases $f=0\text{ or }1/2$, i.e. $\delta M\tau\lessgtr 2\pi$.
We can generalize this analysis to the case where the true hypothesis or the
null hypothesis have a non-trivial LNV fraction $f$. A larger number of events
will then be required to reach the same confidence level. We will not discuss
these cases further in this paper, in order to avoid making the discussion
unnecessarily complicated.

As a final step, for each HNL mass $M$ and ratio
$\left|\Theta_{e}\right|^{2}:\left|\Theta_{\mu}\right|^{2}:\left|\Theta_{\tau}\right|^{2}$,
we compute the squared mixing angles $\left|\Theta_{\alpha}\right|_{i}^{2}(M)$
required to produce $N_{i}(a(M))$ events, thus producing for each true
hypothesis $\mathcal{H}_{i}$ a sensitivity limit, above which SHiP should be
able to exclude the other hypothesis with a probability of at least $1/2$. The
resulting sensitivity plots are presented in section 5.1.

###  B.4 Systematic uncertainties coming from the heavy meson spectrum

For a classifier to generalize well out of sample, i.e. on real-world data,
the distribution used for training should match the true, physical
distribution of features. This is in general not the case, since a simulation
never perfectly represents reality. We can, however, work around this
requirement by explicitly evaluating the classification accuracy over a set of
test distributions which is likely to encompass the true distribution. This
requires knowing and parametrizing the uncertainties coming from the
simulation. We can then obtain a conservative estimate for the classification
accuracy by varying the unknown parameters within their uncertainties, and
taking a lower bound. If this lower bound is high enough, we should still be
able to probe lepton number violation on real data.

At SHiP, the main uncertainty affecting the LNC / LNV classification accuracy
comes from the transverse momentum spread of the heavy meson spectrum, which
is only known with limited accuracy. In order to estimate the actual
sensitivity of SHiP to LNV for a realistic dataset, we therefore compute the
classification accuracy for a family of test sets generated using slightly
different $p_{T}$ spectra, and we take the lowest value as our estimate. The
change in the sensitivity resulting from varying $\langle p_{T}^{2}\rangle$ by
a factor of two up and down with respect to the best-fit value from LEBC-EHS
lebc_d-meson_1987  is shown in fig. 12. The planned charm spectrum
measurements at SHiP should be able to constrain $\langle p_{T}^{2}\rangle$ to
a much better accuracy than the range displayed in the figure.

Figure 12: Effect on the LNV sensitivity ($90\%$ CL) of computing the
classification accuracy on a test set generated with a different $p_{T}$
spectrum compared to the training set, for an HNL coupling to the muon. Black
lines represent the model-selection sensitivity of SHiP for various true
$\langle p_{T}^{2}\rangle$. Here, the training set is always generated with
$\langle p_{T}^{2}\rangle=$1\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$.

Interestingly, when comparing this result with fig. 11, we observe that the
classification accuracy seems to mostly depend on the $\langle
p_{T}^{2}\rangle$ of the test set, but not much on the one used for training.
This suggests that we might be able to safely use the best-fit spectrum for
training without worrying about biasing the results should the true spectrum
turn out to be different, provided that we use a conservative estimate for the
accuracy. In a more comprehensive study, one would likely want to vary
additional parameters related to the spectrum, geometry and simulation.

## References

  * (1) P. Minkowski, _$\mu\to e\gamma$ at a Rate of One Out of $10^{9}$ Muon Decays?_, _Phys. Lett._ 67B (1977) 421–428.
  * (2) M. Gell-Mann, P. Ramond and R. Slansky, _Complex Spinors and Unified Theories_ , _Conf. Proc._ C790927 (1979) 315–321, [1306.4669].
  * (3) R. N. Mohapatra and G. Senjanovic, _Neutrino Mass and Spontaneous Parity Nonconservation_ , _Phys. Rev. Lett._ 44 (1980) 912.
  * (4) T. Yanagida, _Horizontal Symmetry and Masses of Neutrinos_ , _Prog. Theor. Phys._ 64 (1980) 1103.
  * (5) J. Schechter and J. W. F. Valle, _Neutrino Masses in $SU(2)$ $\times$ U(1) Theories_, _Phys. Rev._ D22 (1980) 2227.
  * (6) J. Schechter and J. W. F. Valle, _Neutrino Decay and Spontaneous Violation of Lepton Number_ , _Phys. Rev._ D25 (1982) 774.
  * (7) F. Vissani, _Do experiments suggest a hierarchy problem?_ , _Phys. Rev._ D57 (1998) 7027–7030, [hep-ph/9709409].
  * (8) F. Bezrukov, M. Yu. Kalmykov, B. A. Kniehl and M. Shaposhnikov, _Higgs Boson Mass and New Physics_ , _JHEP_ 10 (2012) 140, [1205.2893].
  * (9) T. Asaka, S. Blanchet and M. Shaposhnikov, _The $\nu$MSM, Dark Matter and Neutrino Masses_, _Phys. Lett._ B631 (2005) 151–156, [hep-ph/0503065].
  * (10) T. Asaka and M. Shaposhnikov, _The $\nu$MSM, Dark Matter and Baryon Asymmetry of the Universe_, _Phys. Lett._ B620 (2005) 17–26, [hep-ph/0505013].
  * (11) J. Beacham et al., _Physics Beyond Colliders at CERN: Beyond the Standard Model Working Group Report_ , 1901.09966.
  * (12) J. P. Chou, D. Curtin and H. J. Lubatti, _New Detectors to Explore the Lifetime Frontier_ , _Phys. Lett._ B767 (2017) 29–36, [1606.06298].
  * (13) D. Curtin et al., _Long-Lived Particles at the Energy Frontier: The MATHUSLA Physics Case_ , _Rept. Prog. Phys._ 82 (2019) 116201, [1806.07396].
  * (14) J. L. Feng, I. Galon, F. Kling and S. Trojanowski, _ForwArd Search ExpeRiment at the LHC_ , _Phys. Rev._ D97 (2018) 035001, [1708.09389].
  * (15) V. V. Gligorov, S. Knapen, M. Papucci and D. J. Robinson, _Searching for Long-lived Particles: A Compact Detector for Exotics at LHCb_ , _Phys. Rev._ D97 (2018) 015023, [1708.09395].
  * (16) V. V. Gligorov, S. Knapen, B. Nachman, M. Papucci and D. J. Robinson, _Leveraging the ALICE/L3 cavern for long-lived particle searches_ , _Phys. Rev._ D99 (2019) 015023, [1810.03636].
  * (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, [1801.04207].
  * (18) SHiP collaboration, _A facility to Search for Hidden Particles (SHiP) at the CERN SPS_ , 1504.04956.
  * (19) SHiP collaboration, _A facility to Search for Hidden Particles at the CERN SPS: the SHiP physics case_ , _Rept. Prog. Phys._ 79 (2016) 124201, [1504.04855].
  * (20) SHiP collaboration, _Sensitivity of the SHiP experiment to Heavy Neutral Leptons_ , _JHEP_ 04 (2019) 077, [1811.00930].
  * (21) M. Shaposhnikov, _A Possible Symmetry of the $\nu$MSM_, _Nucl. Phys._ B763 (2007) 49–59, [hep-ph/0605047].
  * (22) J. Kersten and A. Yu. Smirnov, _Right-Handed Neutrinos at CERN LHC and the Mechanism of Neutrino Mass Generation_ , _Phys. Rev._ D76 (2007) 073005, [0705.3221].
  * (23) K. Moffat, S. Pascoli and C. Weiland, _Equivalence between massless neutrinos and lepton number conservation in fermionic singlet extensions of the Standard Model_ , 1712.07611.
  * (24) G. Anamiati, M. Hirsch and E. Nardi, _Quasi-Dirac neutrinos at the LHC_ , _JHEP_ 10 (2016) 010, [1607.05641].
  * (25) M. Drewes, J. Klarić and P. Klose, _On Lepton Number Violation in Heavy Neutrino Decays at Colliders_ , _JHEP_ 11 (2019) 032, [1907.13034].
  * (26) L. Canetti, M. Drewes, T. Frossard and M. Shaposhnikov, _Dark Matter, Baryogenesis and Neutrino Oscillations from Right Handed Neutrinos_ , _Phys. Rev._ D87 (2013) 093006, [1208.4607].
  * (27) C. Arbelaéz, C. Dib, I. Schmidt and J. C. Vasquez, _Probing the Dirac or Majorana nature of the Heavy Neutrinos in pure leptonic decays at the LHC_ , _Phys. Rev._ D97 (2018) 055011, [1712.08704].
  * (28) C. O. Dib, C. S. Kim and K. Wang, _Signatures of Dirac and Majorana sterile neutrinos in trilepton events at the LHC_ , _Phys. Rev._ D95 (2017) 115020, [1703.01934].
  * (29) P. Hernández, J. Jones-Pérez and O. Suarez-Navarro, _Majorana vs Pseudo-Dirac Neutrinos at the ILC_ , _Eur. Phys. J._ C79 (2019) 220, [1810.07210].
  * (30) G. Cvetic, C. Dib and C. S. Kim, _Probing Majorana neutrinos in rare $\pi^{+}\to e^{+}e^{+}\mu^{-}\nu$ decays_, _JHEP_ 06 (2012) 149, [1203.0573].
  * (31) A. B. Balantekin, A. de Gouvêa and B. Kayser, _Addressing the Majorana vs. Dirac Question with Neutrino Decays_ , _Phys. Lett._ B789 (2019) 488–495, [1808.10518].
  * (32) H. K. Dreiner, H. E. Haber and S. P. Martin, _Two-component spinor techniques and Feynman rules for quantum field theory and supersymmetry_ , _Phys. Rept._ 494 (2010) 1–196, [0812.1594].
  * (33) C. Giunti and C. W. Kim, _Fundamentals of Neutrino Physics and Astrophysics_.  2007\.
  * (34) T. Takagi _Japanese J. Math._ 1 (1927) 83.
  * (35) R. Sachs, _Interference phenomena of neutral K mesons_ , _Annals of Physics_ 22 (1963) 239 – 262.
  * (36) C. Giunti, C. W. Kim, J. A. Lee and U. W. Lee, _On the treatment of neutrino oscillations without resort to weak eigenstates_ , _Phys. Rev._ D48 (1993) 4310–4317, [hep-ph/9305276].
  * (37) M. Beuthe, _Oscillations of neutrinos and mesons in quantum field theory_ , _Phys. Rept._ 375 (2003) 105–218, [hep-ph/0109119].
  * (38) E. K. Akhmedov and J. Kopp, _Neutrino Oscillations: Quantum Mechanics vs. Quantum Field Theory_ , _JHEP_ 04 (2010) 008, [1001.4815].
  * (39) E. K. Akhmedov and A. Yu. Smirnov, _Neutrino oscillations: Entanglement, energy-momentum conservation and QFT_ , _Found. Phys._ 41 (2011) 1279–1306, [1008.2077].
  * (40) E. Akhmedov, _Quantum mechanics aspects and subtleties of neutrino oscillations_ , in _International Conference on History of the Neutrino: 1930-2018 Paris, France, September 5-7, 2018_ , 2019, 1901.05232.
  * (41) S. Eijima, M. Shaposhnikov and I. Timiryasov, _Parameter space of baryogenesis in the $\nu$MSM_, _JHEP_ 07 (2019) 077, [1808.10833].
  * (42) D. Boyanovsky, _Nearly degenerate heavy sterile neutrinos in cascade decay: mixing and oscillations_ , _Phys. Rev._ D90 (2014) 105024, [1409.4265].
  * (43) G. Cvetič, C. S. Kim, R. Kögerler and J. Zamora-Saá, _Oscillation of heavy sterile neutrino in decay of $B\to\mu e\pi$_, _Phys. Rev._ D92 (2015) 013015, [1505.04749].
  * (44) S. Antusch, E. Cazzato and O. Fischer, _Resolvable heavy neutrino–antineutrino oscillations at colliders_ , _Mod. Phys. Lett._ A34 (2019) 1950061, [1709.03797].
  * (45) G. Cvetič, A. Das and J. Zamora-Saá, _Probing heavy neutrino oscillations in rare $W$ boson decays_, _J. Phys._ G46 (2019) 075002, [1805.00070].
  * (46) P. S. Bhupal Dev, R. N. Mohapatra and Y. Zhang, _CP Violating Effects in Heavy Neutrino Oscillations: Implications for Colliders and Leptogenesis_ , _JHEP_ 11 (2019) 137, [1904.04787].
  * (47) G. Cvetič, A. Das, S. Tapia and J. Zamora-Saá, _Measuring the heavy neutrino oscillations in rare W boson decays at the Large Hadron Collider_ , _J. Phys._ G47 (2020) 015001, [1905.03097].
  * (48) S. Tapia and J. Zamora-Saá, _Exploring CP-Violating heavy neutrino oscillations in rare tau decays at Belle II_ , 1906.09470.
  * (49) A. Abada, C. Hati, X. Marcano and A. M. Teixeira, _Interference effects in LNV and LFV semileptonic decays: the Majorana hypothesis_ , _JHEP_ 09 (2019) 017, [1904.05367].
  * (50) LHCb collaboration, _Search for Majorana neutrinos in $B^{-}\to\pi^{+}\mu^{-}\mu^{-}$ decays_, _Phys. Rev. Lett._ 112 (2014) 131802, [1401.5361].
  * (51) ATLAS collaboration, _Search for heavy Majorana neutrinos with the ATLAS detector in pp collisions at $\sqrt{s}=8$ TeV_, _JHEP_ 07 (2015) 162, [1506.06020].
  * (52) CMS collaboration, _Search for heavy Majorana neutrinos in $\mu^{\pm}\mu^{\pm}+$ jets events in proton-proton collisions at $\sqrt{s}$ = 8 TeV_, _Phys. Lett._ B748 (2015) 144–166, [1501.05566].
  * (53) CMS collaboration, _Search for heavy Majorana neutrinos in e ±e±\+ jets and e± $\mu^{\pm}$\+ jets events in proton-proton collisions at $\sqrt{s}=8$ TeV_, _JHEP_ 04 (2016) 169, [1603.02248].
  * (54) LEBC-EHS collaboration, _$D$ Meson Production From 400 GeV/$cpp$ Interactions_, _Phys. Lett._ B189 (1987) 476.
  * (55) SHiP collaboration, _Measurement of associated charm production induced by 400 GeV/c protons_ , Tech. Rep. CERN-SPSC-2017-033. SPSC-EOI-017, CERN, Geneva, Oct, 2017.
  * (56) K. Bondarenko, A. Boyarsky, D. Gorbunov and O. Ruchayskiy, _Phenomenology of GeV-scale Heavy Neutral Leptons_ , _JHEP_ 11 (2018) 032, [1805.08567].
  * (57) P. Ballett, T. Boschi and S. Pascoli, _Heavy Neutral Leptons from low-scale seesaws at the DUNE Near Detector_ , 1905.00284.
  * (58) SHiP collaboration, B. Hosseini and W. M. Bonivento, _Particle Identification tools and performance in the SHiP Experiment_ , Jul, 2017.
  * (59) G. Ke, Q. Meng, T. Finley, T. Wang, W. Chen, W. Ma et al., _Lightgbm: A highly efficient gradient boosting decision tree_ , in _Advances in Neural Information Processing Systems 30_ (I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan et al., eds.), pp. 3146–3154.  Curran Associates, Inc., 2017.
  * (60) Microsoft Corporation, _LightGBM_ , _GitHub repository_ (2016–2019) .
  * (61) J. Neyman and E. S. Pearson, _On the Problem of the Most Efficient Tests of Statistical Hypotheses_ , _Phil. Trans. Roy. Soc. Lond._ A231 (1933) 289–337.
  * (62) S. S. Wilks, _The Large-Sample Distribution of the Likelihood Ratio for Testing Composite Hypotheses_ , _Annals Math. Statist._ 9 (1938) 60–62.
  * (63) M. Drewes, _On the Minimal Mixing of Heavy Neutrinos_ , 1904.11959.
  * (64) F. James, _Monte-Carlo phase space_ , 1968.  10.5170/CERN-1968-015.
  * (65) P. Ilten, _Electroweak and Higgs Measurements Using Tau Final States with the LHCb Detector_ , Ph.D. thesis, University Coll., Dublin, 2013-09-06.  1401.4902.
  * (66) S. Algeri, J. Aalbers, K. Dundas Morå and J. Conrad, _Searching for new physics with profile likelihoods: Wilks and beyond_ , 1911.10237.