We study pseudo-Dirac pairs of two almost mass degenerate sterile Majorana neutrinos which generate light neutrino masses via a low-scale seesaw mechanism.
These pseudo-Dirac heavy neutral leptons can oscillate between interaction eigenstates that couple to leptons and antileptons and thus generate oscillations between lepton number conserving and lepton number violating processes.
With the we introduce a minimal framework capable of describing the dominant features of low-scale seesaws at colliders and present a FeynRules implementation usable in Monte Carlo generators.
Additionally, we extend MadGraph to be capable of simulating heavy neutrino-antineutrino oscillations and show results obtained from such simulations.

§ INTRODUCTION

The discovery of flavour oscillations between neutrinos <cit.> implies that they have nonzero masses.
One possible extension of the able to generate neutrino masses consists of adding sterile neutrinos, that are singlets under all symmetries, to its particle content <cit.>.
Testing whether these sterile neutrinos are indeed the missing piece that explains the light neutrino masses is one of the key questions towards a more complete theory of elementary particles.

With sterile neutrinos added to the , there are two main routes resulting in non-zero neutrino masses.
In the first case, only Dirac mass terms are present.
The Yukawa coupling term is formed by sterile neutrinos, the Higgs field and the lepton doublets, analogous to the mechanism generating charged fermion masses.
After this Yukawa term generates Dirac neutrino masses.
However, in order to yield the phenomenologically required small masses for the light neutrinos, the size of the Yukawa couplings has to be tiny $\vec y \lessapprox \order{\expno{-12}}$.

Since the sterile neutrinos are singlets, it is possible to add a Majorana mass term to the Lagrangian <cit.>.
Realising light Dirac neutrinos therefore requires a mechanism to eliminate this term and to enforce exact .
As soon as some is present, the light neutrinos are of Majorana-type <cit.>.
One way to probe these considerations are searches for decay <cit.>, since its observation would prove .
The two main routes for non-zero neutrino masses are depicted in the upper part of <ref> and in this work we focus on the case that the light neutrinos have Majorana masses.

As an alternative to decays, one can also search for induced by the sterile neutrinos at colliders.
While this sounds promising at first sight, one can argue on general grounds that , observable at the , would lead to too heavy light neutrino masses, and thus it should be impossible to observe <cit.>.
An effect that has not been taken into account in such considerations are <cit.>, see also <cit.>.
In particular, their interplay with the potential longevity of the heavy neutrinos are the main subject of this paper.
Despite the smallness of the terms in the Lagrangian, can introduce processes at the same order as processes, depending on the lifetime and oscillation length of the sterile neutrinos.
Since the oscillations are an interference phenomenon they are able to probe mass splittings of heavy neutrinos so small that they would otherwise be unobservable in collider studies.
Including the correctly in collider studies is thus a crucial aspect when simulating seesaw extensions of the with the goal to explain the light neutrino masses.

This paper is organised as follows: In <ref> we describe the possible domains of seesaw models and argue that collider testable seesaw models are protected by a resulting in pseudo-Dirac pairs of heavy neutrinos.
Afterwards, in <ref> we give the relevant results for the description of in the external wave packet formalism and derive the integrated effects of these oscillations.
Subsequently, we introduce the <cit.> in <ref>, first in the symmetric limit, then extended by small symmetry violating terms, and finally as the with the minimal set of parameters able to describe the dominant collider effects of low-scale type I seesaws.
In <ref> we introduce the FeynRules model file of the and describe the necessary steps to extend MadGraph to be able to simulate .
In <ref> we present selected results obtained from a study using this software implementation.
Finally, we conclude in <ref>.
Additionally, we comment in <ref> on the discussion about whether it is possible to distinguish Majorana and Dirac particles from each other using their decay width.
Last but not least, the code for the MadGraph patch is presented in <ref>.

§ SEESAW MODELS

[Comparison between the different possible choices for active and sterile neutrinos]
Comparison between the different possible choices for active and sterile neutrinos discussed in the text.
Note that for collider accessible heavy neutrinos generating Majorana masses for the light neutrinos via a low-scale seesaw with couplings far above the naive seesaw line only pseudo-Dirac pairs of two nearly mass degenerate Majorana $\DOFs$ are a viable option.

With $i = 1, \dots, n$ sterile neutrinos $N_i$ and both Dirac and Majorana mass terms extending the Lagrangian, one arrives at a general theoretical framework where the sterile neutrino Lagrangian
For ease of notation, the sterile neutrinos are introduced as left-chiral $\DOFs$.
Note that any left-chiral field can also be described by the charge conjugate of a right-chiral field and vice versa.

\begin{equation} \label{eq:generic seesaw}
\mathcal L_N =
- \sum_i y_\alpha^{(i)} \widetilde H^\dagger \widebar \ell_\alpha N_i^c
- \frac12 \sum_i m_M^{(i)} \widebar N_i^{} N_i^c
+ \HC \,,
\end{equation}

is added to the Lagrangian.
Here $H$ and $\ell$ are the Higgs and lepton $\SU(2)$ doublets, respectively, $\vec y^{(i)}$ is the neutrino Yukawa coupling vector $\row{y_1^{(i)},y_2^{(i)},y_3^{(i)}}^\trans$, and $m_M^{(i)}$ is the Majorana mass parameter.
We indicate unindexed vectorial quantities by using bold face font.
In order to illustrate how small neutrino masses can emerge from this framework, we consider the minimal case of two sterile neutrinos $n=2$.
It is minimal since neutrino flavour oscillations require at least two of the light neutrino masses to be nonzero, which in turn requires at least two sterile neutrino .
The following discussion can of course be generalised to larger $n$, which is necessary for models such as the inverse seesaw to be phenomenological viable.
After , the Lagrangian can be written in the diagonal basis for the Majorana masses $m_M^{}$,

\begin{equation} \label{eq:mass Lagrangian}
\mathcal L_N =
- m_{D\alpha}^{(1)} \widebar \nu_\alpha N_1^c
- m_{D\alpha}^{(2)} \widebar \nu_\alpha N_2^c
- \frac12 m_M^{(1)} \widebar N_1^{} N_1^c
- \frac12 m_M^{(2)} \widebar N_2^{} N_2^c
+ \HC \,,
\end{equation}

where $\vec m_D^{} = \vec y v$ with the Higgs $v \approx \unit[174]{GeV}$ describes the Dirac mass contribution.
When $m_M^{}$ is sufficiently larger than $\vec m_D^{}$, the light neutrino mass matrix is given by the seesaw relation, which takes the form

\begin{equation} \label{eq:seesaw}
M_\nu = \frac{\vec m_D^{(1)} \otimes \vec m_D^{(1)}}{m_M^{(1)}} + \frac{\vec m_D^{(2)} \otimes \vec m_D^{(2)}}{m_M^{(2)}} \,.
\end{equation}

With these ingredients, there are three limiting cases allowing to arrive at the observed small neutrino masses:

[Limiting cases of seesaw models]
Depiction of the three limiting seesaw regimes described in the text.
The phenomenologically accessible parameter space of low-scale seesaws with large enough coupling is connected to the symmetry protected corner and indicated by dashed lines.

* High-scale seesaw limit

Here neutrino masses are small because the Majorana masses $m_M^{}$ of the sterile neutrinos are large.
The neutrino Yukawa couplings $\vec y$ can be large in this case — up to $\order 1$ for $m_M^{}$ close to the scale of grand unification.
Due to the high mass scale involved this limit of the seesaw mechanism can not be probed directly at colliders.

* Small coupling seesaw limit

It is possible to simultaneously lower the sterile neutrino mass scale $m_M^{}$ and the size of the neutrino Yukawa couplings $\vec y$ without changing the generated light neutrino masses $M_\nu$ in <ref>.
Envisioned future collider experiments such as the -$ee$ can probe this limit for some choices of parameters, but for most cases the couplings are too small for direct tests.

* Symmetry protected seesaw limit

The third possibility emerges when the two terms in <ref> almost cancel.
Then $\vec y$ can be large enough and simultaneously $m_M^{}$ can be small enough, such that the heavy neutrinos are within the reach of collider experiments.
This cancellation can be protected by a that generalises the lepton number $L$ of the and ensures that $M_\nu$ is equal to zero in the symmetry conserving limit.

One can achieve this cancellation also by tuning of parameters, without relying on a symmetry.
However, this mechanism of cancellation is unstable under radiative corrections, <cit.>.
We will not consider this possibility further here, noting only that in this case it is expected that the effects would be generically unsuppressed.

The various manifestations of the seesaw mechanism are schematically depicted in <ref>.
The three corners represent the three limiting cases discussed above.
However, also the cases in between are viable options.
For the Majorana mass scales $m_M^{}$, there is a maximal value above which it is no longer possible to test this sterile neutrino directly at a given collider experiment.
This means the observable seesaw models are low-scale models opposite to the high-scale seesaw corner.
Among the low-scale seesaw models, the small coupling limit is also not testable at the or even the .
The potentially testable region of the (type I) seesaw mechanism is thus the area inside the dashed lines in <ref>, it is bound to have a certain degree of symmetry protection, if one wants to avoid tuning of the parameters.

Although the model of two sterile neutrinos with exact serves as a valuable starting point for the construction of viable seesaw models with small symmetry breaking, it is itself not able to generate neutrino masses due to the unbroken symmetry.
With the symmetry intact, the two Majorana are mass degenerate and combine precisely to a single Dirac particle.
However, a small amount of symmetry breaking not only generates small neutrino masses but at the same time causes a small mass splitting between the two Majorana .
Such pseudo-Dirac heavy neutrinos can exhibit potentially detectable at collider experiments.
For large symmetry breaking this feature vanishes as the mass splitting becomes too big and decoherence sets in, so that the two sterile neutrinos appear as two separate Majorana particles.
At the same time the generated neutrino masses become too large if the model has no additional mechanism to prevent this from happening.
We depict these considerations in the lower part of <ref>.

Although realistic low-scale seesaw models predict pseudo-Dirac heavy neutrinos the majority of searches for have been performed in either the pure Dirac or single Majorana scenario.
In an effort to distinguish between these two models it has sometimes been argued that one can discriminate Majorana from Dirac using their decay width.
In <ref> we describe in detail why this is not the case.
The main insight is that the factor of two appearing when comparing the decay widths counts the number of Majorana forming the observed .
However, two Majorana particles can only be described as a Dirac particle when their Yukawa couplings have a relative phase of $-i$, which can not be determined using the decay width, <ref>.

[Feynman diagram visualising the oscillations]
A neutrino interaction eigenstate is produced together with an antilepton.
After some time the superposition of mass eigenstates $n_i$ has oscillated into a heavy neutrino $N$ or antineutrino $\widebar N$ interaction eigenstate that decays into an antilepton $l^+$ or lepton $l^-$, such that the total process is or .

Short oscillation length.
Intermediate oscillation length.
Long oscillation length.

[Depiction of oscillations lengths]
Oscillations of the probability $P$ that the process shown in <ref> is or as a function of the proper time $\tau$.
The oscillation length depends on the mass splitting $\Delta m$, the overall decay depends on the decay width $\Gamma$, and the decoherence depends on the damping parameter $\lambda$, all of which are defined in <ref>.
Panel fig:short shows a short oscillation length with $\Delta m/\Gamma = 100$, panel fig:Intermediate shows an intermediate oscillation length with $\Delta m/\Gamma = 10$, and panel fig:long shows a long oscillation length with $\Delta m/\Gamma = 1$.
All oscillations are shown for $\Gamma = \frac{\unit{meV}}{25}$ and $\lambda = \nicefrac15$ and $R_{ll}$ is defined in <ref>.

Generically, neutral particles can oscillate into their antiparticles <cit.> as known from meson oscillations <cit.>.
The distinction between particle and antiparticle in the case of sterile neutrinos can be made by distinguishing a neutrino interaction eigenstate by the charge of the associated lepton it is produced with.
If it is produced together with a lepton $l^-$, it is an antineutrino and if it is produced together with an antilepton $l^+$ it is a neutrino.
The corresponding mass eigenstates $n_i$ interfere when propagating, which results in an oscillation between different interaction eigenstates as a function of travelling distance.
Heavy neutrinos $N$ and antineutrinos $\widebar N$ are defined as the projection of neutrino and antineutrino states onto the heavy mass eigenstates.
This mechanism leads to .
In the collider testable region of <ref> the mass splitting can be small enough for the to be observable and of macroscopic length.
If are present one expects processes such as the one shown in <ref> to occur, leading to patterns as depicted in <ref>.
Furthermore, in cases in which the parameters do not allow for macroscopic oscillations to be resolvable an integrated effect could still be measured.

§.§ External wave packets

In order to predict the oscillatory behaviour of new particles a holistic framework is necessary that, in contrast to the simplified framework, is able to not only capture potential oscillations but can also predict the potential damping of these oscillations due to the loss of coherence of the mass eigenstate superposition.
Compared to the plane wave description of particle oscillations, the external wave packet approach discussed in <cit.> and adapted to the case of in <cit.> allows the derivation of the phenomenon free of contradictions.
This is due to the fact that the inherent uncertainty of wave packets in momentum and position space allows the simultaneous production of several on-shell mass eigenstates, that can subsequently interfere to produce oscillations.
Furthermore, the approximate localisation of wave packets is necessary to introduce the notion of a travelled distance and time, which is not possible using infinitely extended plane waves.
Additionally, several effects potentially leading to the decoherence of the mass eigenstate superposition are included and discussed under the name of observability conditions <cit.>.
External refers to the fact that only the external particles are explicitly assumed to be wave packets, whereas in an intermediate wave packed approach the intermediate particles, such as in this case the heavy neutrinos, are directly described by wave packets.
However, even in the external wave packet approach it is possible to interpret the intermediate particles as wave packets where some clarifications about the causal nature of this interpretation can be found in <cit.>.

For the external wave packet models the assumptions about the shape and the widths of the involved wave packets are a major source of uncertainty when deducing results.
Also the formalism discussed in <cit.> was mainly developed having long-lived quasi mass degenerate particles in mind.
In particular the theorem of Jacob-Sachs <cit.> is only valid above a certain time threshold, below which additional correction have to be taken into account.
In <cit.> the regions as well as the effects of the above mentioned corrections are discussed in detail and it is shown that in phenomenological studies of nearly mass degenerate, quasi long-lived heavy neutrinos, the effects of the corrections can be neglected for a range of external wave packet widths spanning orders of magnitude.

§.§ Oscillations at

In the case of slightly broken governing the physics of pseudo-Dirac heavy neutrinos the complete contributions to the can be described using only two parameters in addition to the seesaw model in the limit <cit.>.
Hence instead of working with the whole family of realistic seesaw models and their parameters one can instead introduce, in addition to the heavy neutrino mass $m$ and the active-sterile mixing parameter $\vec\theta$, just the mass splitting of the heavy neutrinos $\Delta m$ and the damping parameter $\lambda$.
As an effective parameter $\Delta m$ captures the effects induced by the breaking parameters and $\lambda$ entails the decoherence effects appearing in certain parameter regions of the description of neutral particle oscillations.
To , the formulae for the oscillation probabilities as a function of the proper time $\tau$ are

\begin{equation}
P^{\nicefrac{\LNC}{\LNV}}_\text{osc}(\tau) = \frac{1 \pm \cos\left(\Delta m \tau \right) \exp(-\lambda)}{2}\,,
\end{equation}

hence the oscillation time is given by $\tau_\text{osc} = 2 \pi / \Delta m$.
The overall decay of the heavy neutrino is given by the usual probability density for the decay of unstable particles

\begin{equation}
P_\text{decay}(\tau) = - \dv \tau \exp\left(- \Gamma \tau\right) = \Gamma \exp\left(- \Gamma \tau\right) \,.
\end{equation}

as function of the decay width $\Gamma$.
Therefore, the total probability for an unstable and oscillating particle to decay in a proper time window is given by

\begin{equation} \label{eq:probability displacement}
P_{ll}^{\nicefrac{\LNC}{\LNV}}(\tau_{\min}, \tau_{\max}) = \int_{\tau_{\min}}^{\tau_{\max}} P^{\nicefrac{\LNC}{\LNV}}_\text{osc}(\tau) P_\text{decay}(\tau) \d \tau \,,
\end{equation}

and the number of expected events in a collider experiment is then

\begin{equation} \label{eq:event numbers}
N^{\nicefrac{\LNC}{\LNV}} = \mathscr L \sigma \BR \int D(\vartheta, \gamma) P^{\nicefrac{\LNC}{\LNV}}_{ll}(\tau_{\min}(\vartheta, \gamma), \tau_{\max}(\vartheta, \gamma)) \d\vartheta \d\gamma\,,
\end{equation}

where the factor $D(\vartheta, \gamma)$ accounts for the probability density that the heavy neutrino has Lorentz factor $\gamma$ and is produced with an angle $\vartheta$ with respect to the beam axis, $\mathscr L$ is the luminosity of the collider, and $\sigma$ and $\BR$ are the sterile neutrino production cross section and branching ratio of the process under consideration.
The parameters $\tau_{\min}$ and $\tau_{\max}$ are defined by the detector geometry when transitioning from proper time coordinates to the lab frame via $\tau(\vartheta,\gamma) = (\gamma^2-1)^{-\nicefrac12} L(\vartheta)$.

§.§ Integrated effect

In an experimental situation in which it is not possible to resolve the oscillation patterns, an integrated effect could still be measurable by comparing the number of events with opposite- and same-sign leptons originating from processes such as the one presented in <ref>.
Evaluating the integral (<ref>) in order to obtain the probability of and decays between the minimal decay length $\tau_{\min}$ and the maximal decay length $\tau_{\max}$ results in the difference

\begin{equation} \label{eq:finite detector probability}
P_{ll}^{\nicefrac{\LNC}{\LNV}}(\tau_{\min}, \tau_{\max})
= \Gamma \frac{P^{\nicefrac{\LNC}{\LNV}}(\tau_{\max}) - P^{\nicefrac{\LNC}{\LNV}}(\tau_{\min})}2 \,,
\end{equation}

where the indefinite integral is given by

\begin{equation}
P^{\nicefrac{\LNC}{\LNV}}(\tau) = P(\tau, \Gamma, 0) \pm \frac{P(\tau, \Gamma_-, \lambda) + P(\tau, \Gamma_+, \lambda)}2 \,,
\end{equation}

\begin{align}
P(\tau, \Gamma, \lambda) &= \int e^{- \lambda - \Gamma \tau } \d \tau = - \frac{e^{- \lambda - \Gamma \tau}}{\Gamma} \,, &
\Gamma_\pm &= \Gamma \pm i \Delta m \,.
\end{align}

$R_{ll}^\text{obs}(R_{ll},\tau_{\min}/\tau_\text{osc})$ while $\tau_{\max}\to\infty$.
$R_{ll}^\text{obs}(R_{ll},\tau_{\max}/\tau_\text{osc})$ while $\tau_{\min}\to0$.

[$R_{ll}$ in finite detectors]
Impact of the finite detector size on the measurement of $R_{ll}$.
The plots show the observable $R_{ll}^\text{obs}$ as a function of $\tau/\tau_\text{osc}$ and the theoretical $R_{ll}$, calculated by taking simultaneously the limits (<ref>).
In panel fig:Rll tau min $\tau_{\min}$ is variable and $\tau_{\max}\to\infty$ while in panel fig:Rll tau max $\tau_{\max}$ is variable and $\tau_{\min} \to 0$.
Note that only in the region where the contours are vertical $R_{ll}$ coincides with $R_{ll}^\text{obs}$ and that $R_{ll}^\text{obs}$ is not bound to be smaller or equal to one and even has in poles in panel fig:Rll tau min.

In the limit that the experiment can observe all decays from the origin to infinity and under the assumption that the parameter point under consideration does allow to neglect the decoherence effects, which is equivalent to

\begin{align} \label{eq:Rll limits}
\tau_{\min} &\to 0 \,, &
\tau_{\max} &\to \infty \,, &
\lambda &\to 0 \,,
\end{align}

this expression simplifies to

\begin{equation}
P_{ll}^{\nicefrac{\LNC}{\LNV}} = \frac12
\begin{cases}
\frac{\Gamma^2}{\Delta m^2 + \Gamma^2} + 1 & \text{for \LNC} \\
\frac{\Delta m^2}{\Delta m^2 + \Gamma^2} & \text{for \LNV}
\end{cases}
\,.
\end{equation}

Therefore, the ratio between the two decay modes is in this case given by

\begin{equation} \label{eq:ratio}
R_{ll} =
\frac{P_{ll}^{\LNV}}{P_{ll}^{\LNC}} =
\frac{\Delta m^2}{\Delta m^2 + 2 \Gamma^2} \,,
\end{equation}

which matches the result of <cit.>.
While a Dirac heavy neutrino would have $R_{ll} = 0$ and a Majorana heavy neutrino would have $R_{ll} = 1$, a realistic pseudo-Dirac heavy neutrino can have any value in between as is shown in <ref>, where the dependence of $R_{ll}$ on $\flatfrac{\Delta m}{\Gamma}$ is depicted and compared to a simulation.
Furthermore, the bands of $R_{ll} \in [0.1,0.9]$ for five introduced in <ref> are given in <ref>.

However, taking the finite size of the experimental setup into account by integrating only over the fiducial detector size in (<ref>) changes the picture drastically.
In this case the parameters $\tau_{\min}$ and $\tau_{\max}$ in

\begin{equation} \label{eq:observable ratio}
R_{ll}^\text{obs}(\tau_{\min}, \tau_{\max}) =
\frac{P_{ll}^{\LNV}(\tau_{\min}, \tau_{\max})}{P_{ll}^{\LNC}(\tau_{\min}, \tau_{\max})} \,,
\end{equation}

can not be taken to be zero and infinity, respectively.
The consequences of this effect in the limit of vanishing damping and in terms of the proper time are presented in <ref>.
In this representation vertical contour lines indicate negligible deviation of $R_{ll}^\text{obs}$ from $R_{ll}$.
For a $\tau_{\min}$ an oscillatory pattern becomes apparent that has its minima for cuts in the vicinity of multiples of the oscillation length.
The dependence on $\tau_{\max}$ is less sever and only occurs close to an $R_{ll}$ of one.
The two extreme limiting cases are

\begin{align}
R_{ll}^\text{obs} &=
\begin{dcases}
\tan^2{\frac{r^{}_{\min}}{2}} & \text{for } R_{ll} \to 0 \text{ and } \tau_{\max} \to \infty \\
\frac{r^{}_{\max} - \sin r^{}_{\max}}{r^{}_{\max} + \sin r^{}_{\max}} & \text{for } R_{ll} \to 1 \text{ and } \tau_{\min} \to 0
\end{dcases} \,, &
r^{}_m &= \Delta m \tau_m \,.
\end{align}

While the effect for $R_{ll} \to 1$ and $\tau_{\min} \to 0$ is damped for large $\tau_{\max}$ the effect for $R_{ll} \to 0$ and $\tau_{\max} \to \infty$ remains undamped when varying $\tau_{\min}$ as long as the Lorentz factor distribution and the decoherence can be neglected <cit.>.

§.§ Oscillations in the lab frame

[Experimentally reachable oscillations lengths]
Oscillation length $L_\text{osc}$ as a function of the heavy neutrino mass $m$ and the mass splitting $\Delta m$.
The Lorentz factor is estimated using $W$ bosons at rest (<ref>).
The horizontal lines correspond to the mass splitting appearing in the two minimal linear seesaw $\BMs$ given in <ref>.

When boosting the oscillations described in <ref> into the lab frame two main effects have to be considered.
Firstly, the decay length in the lab frame is on an event per event basis increased by the Lorentz factor.
This well known effect allows to measure oscillations with $c\tau < L_{\min}$ and helps to compensate the short oscillations length of some realistic .
For processes such as the one depicted in <ref> the Lorentz factor of the heavy neutrino $\gamma$ can be estimated to be

\begin{equation} \label{eq:boost estimate}
\gamma \approx \frac{m_W^2 + m^2}{2 m_W^{} m^{}}\,,
\end{equation}

which relies on the assumptions that the initial $W$ boson decays at rest and neglects the mass of the prompt muon.
The resulting resolvable oscillation length as a function of the heavy neutrino mass $m$ and the mass splitting $\Delta m$ is shown in <ref>.

Secondly, the oscillation pattern is washed out by the event-dependent Lorentz factor.
Therefore, any study hoping to resolve needs to measure enough observables to reconstruct the Lorentz factor.
This makes studies relying on processes with final state neutrinos extremely challenging and is the reason why we focus on semileptonic final states as shown in <ref>.

As detectable heavy neutrinos cannot be too heavy or too weakly coupled they must originate from a model with an approximate and form pseudo-Dirac pairs as discussed in <ref>.
In the following we assume that one pseudo-Dirac pair dominates the collider phenomenology.
We will introduce first the type I seesaw in the limit of exact symmetry conservation in <ref>.
This model builds on two Majorana that form one exact Dirac particle.
Subsequently, we add small symmetry breaking terms in <ref> recovering models such as the linear and inverse seesaw.
After comparing this model with the terms governing the summarised in <ref> we introduce a phenomenological model with a minimal number of parameters that is able to describe the relevant observables in <ref>.

§.§ Seesaw in the symmetric limit

In the symmetric limit of the <cit.> the two heavy Majorana neutrinos are protected by a .
One simple choice of charges for the protecting the lepton number $L$ is given by

\begin{equation} \label{eq:charges}
\begin{array}{*4c} \toprule
& \ell & N_1 & N_2 \\ \midrule
L & +1 & -1 & +1 \\
\bottomrule\end{array}
\,,
\end{equation}

with all other fields having a charge of zero.
In this case the general Lagrangian (<ref>) takes the form

\begin{equation} \label{eq:symmetric Lagrangian}
\mathcal L_{\SPSS}^L =
\widebar N_i \ispd N_i
- y_{\alpha1} \widetilde H^\dagger \widebar \ell_\alpha N_1^c
- \widebar N_1 m_M^{} N_2^c
+ \dots + \HC \,,
\end{equation}

where $N_1$ and $N_2$ are taken to be left-chiral sterile neutrinos.
The ellipses capture contributions from additional sterile neutrinos, which are assumed to be heavier or much weaker coupled than the explicitly denoted pair and are expected to contribute only sub-dominantly to the collider phenomenology of the model.
After the neutrino mass Lagrangian of the interaction eigenstates $n = \row{\nu_e, \nu_\mu, \nu_\tau, N_1, N_2}^\trans$ can be written as

\begin{equation}
\mathcal L_\text{mass} = - \frac12 \widebar n M_n n^c + \HC \,,
\end{equation}

where the mass matrix is given by

\begin{equation}
M_n =
\begin{pmatrix}
0 & M_{3\times2} \\
M_{3\times2}^\trans & M_N
\end{pmatrix}
\begin{pmatrix}
0 & \vec m_D^{} & 0 \\
\vec m_D^\trans & 0 & m_M^{} \\
0 & m_M^{} & 0
\end{pmatrix}
\,,
\end{equation}

with $\vec m_D^{} = \vec y_1 v$ being the Dirac mass term.
The mass matrix can be approximately diagonalised using a Takagi decomposition

\begin{align} \label{eq:neutrino diagonlisation}
D_n &= U_n^\trans M_n U_n \,, &
U^\dagger U = \mathbb 1 \,,
\end{align}

following the steps presented in <cit.>.
Up to second order in the active-sterile mixing parameter

\begin{equation} \label{eq:active-sterile mixing}
\vec \theta = \frac{\vec m_D^{}}{m_M^{}},
\end{equation}

where $\vec \theta = \row{\theta_e,\theta_\mu,\theta_\tau}^\trans$,
the light neutrinos are massless while the heavy neutrinos are mass degenerate with masses

\begin{equation}\label{eq:mass degeneracy}
m_4 = m_5 = m_M^{} \left(1 + \frac12 \abs{\vec \theta}^2\right) + \order*{\abs{\vec \theta}^4} \,,
\end{equation}

and the mixing matrix is given by

\begin{equation}
U_n =
\begin{pmatrix}
U_{\PMNS} & U_{\CL} \\
U_{2\times 3} & U_N
\end{pmatrix}
\begin{pmatrix}
\mathbb 1_{3\times 3} - \frac12 \vec \theta^* \otimes \vec \theta & -\frac{i}{\sqrt2} \vec \theta^* & \frac{1}{\sqrt2} \vec \theta^* \\
0 & \frac{i}{\sqrt2} & \frac{1}{\sqrt2} \\
- \vec \theta & -\frac{i}{\sqrt2}(1 - \frac12 \abs{\vec \theta}^2) & \frac{1}{\sqrt2}(1 - \frac12 \abs{\vec \theta}^2)
\end{pmatrix}
\,,
\end{equation}

where the upper left block is the matrix.
One can show that the upper right $3\times2$ block of the mixing matrix, which relates heavy neutrinos and active neutrinos, has the exact form

\begin{equation}
U_{\CL} = \frac{1}{\sqrt{1 + \abs{\vec \theta}^2}} \row*{-\frac{i}{\sqrt2} \vec \theta^*, \frac{1}{\sqrt2} \vec \theta^*} \,.
\end{equation}

It is therefore possible to absorb higher order corrections in $\vec \theta$ into a rescaling of the Yukawa coupling.
In particular the exact form of $U_{\CL}$ can be recovered from the approximation by rescaling of the mixing parameter

\begin{equation}\label{eq:rescaling}
\vec \theta^* \rightarrow \vec \theta'^* = \frac{\vec \theta^*}{\sqrt{1 + \abs{\vec \theta}^2}}\,.
\end{equation}

Consequently it suffice to expand $U_{\CL}$ to .

Due to the mass degeneracy and phase difference of the heavy neutrinos there are no in the symmetric limit and this theory of two Majorana particles is equivalent to a theory with a single Dirac particle.
For any process with external particles that contains a Feynman diagram with one of the heavy neutrino mass eigenstates, also the other heavy neutrino mass eigenstate contributes.
In an explicit calculation it can be shown that the amplitudes from the different mass eigenstates cancel each other, such that processes are not present, see also <ref>.
We recover the argument laid out in <ref> that a theory relying solely on one exact Dirac heavy neutrino is and incapable of generating neutrino masses.

§.§ Seesaw with small symmetry breaking

In order to introduce neutrino masses the symmetric limit of the can be perturbed by extending the Lagrangian (<ref>) with additional small terms

\begin{equation} \label{eq:broken Lagrangian}
\mathcal L_{\SPSS}^{\cancel L} =
- y_{a2} \widetilde H^\dagger \widebar \ell_a N_2^c
- \mu_M^\prime \widebar N_1^{} N_1^c
- \mu_M^{} \widebar N_2^{} N_2^c
+ \dots + \HC \,,
\end{equation}

where a Yukawa coupling $\vec y_2$ and two Majorana masses $\mu_M^{}$ and $\mu_M^\prime$ are introduced.
The resulting mass matrix is

\begin{equation} \label{eq:broken mass matrix}
M_n =
\begin{pmatrix}
0 & \vec m_D^{} & \vec \mu_D^{} \\
\vec m_D^\trans & \mu_M^\prime & m_M^{} \\
\vec \mu_D^\trans & m_M^{} & \mu_M^{}
\end{pmatrix}
\,,
\end{equation}

where $\vec \mu_D^{} = \vec y_2 v$ is a Dirac mass term.
These additional terms break the symmetry characterised by the charges given in <ref>.
In order to ensure a small breaking of the we require

\begin{equation}
\order*{\frac{\vec \mu_D^{}}{m_M^{}}} = \order*{\frac{\mu_M^{}}{m_M^{}}} = \order*{\frac{\mu_M^\prime}{m_M^{}}} = \epsilon \ll 1\,,
\end{equation}

Since the parameters $\vec \theta$ can be larger than $\order \epsilon$, an expansion in small parameters might contain powers in $\vec \theta$ higher than in $\epsilon$.
The Takagi decomposition (<ref>) of the mass matrix (<ref>) results now in two distinct heavy neutrino masses of

\begin{align} \label{eq:mass splitting}
m_{\nicefrac{4}{5}}^{} &= m_M^{} \left (1 + \frac12 \abs{\vec \theta}^2 \right) \mp \left(\cos(\phi) \abs*{\vec \mu_D^* \vec \theta} + \frac{\abs*{\mu_M^* + \mu_M^\prime}}2 \right) + \order*{\max(\epsilon^2, \epsilon \abs{\vec \theta}^2, \abs{\vec \theta}^4)} \,.
\end{align}

\begin{equation}
\phi =
\begin{cases}
0 & \text{for } \mu_M = \mu_M^\prime = 0 \\
\arg(\vec \mu_D^* \vec \theta (\mu_M^* + \mu_M^\prime)) & \text{otherwise}
\end{cases}
\,.
\end{equation}

The upper right $2\times3$ part of the mixing matrix is given by

\begin{equation}
U_{\CL} =
\row*{\frac{i}{\sqrt2} \left(\vec \theta^* - \vec \theta_\epsilon^*\right),
\frac{1}{\sqrt2} \left( \vec \theta^* + \vec \theta_\epsilon^*\right)
+ \order*{\max(\epsilon \vec \theta, \abs{\vec\theta}^2 \vec \theta)} \,,
\end{equation}

where $\vec \theta_\epsilon = \vec \mu_D^{} / m_M^{}$ is the equivalent to the active-sterile mixing parameter.
In cases where higher order terms in $\vec\theta$ are of the same order as terms of $\order{\epsilon\vec\theta}$, the rescaling (<ref>) can be used to absorb higher order terms into a redefinition of the coupling $\vec \theta$.
This shows that higher order terms in $\vec \theta$ yield no qualitative new effects and the expansion of mixing matrix is valid up to $\order{\epsilon\vec\theta}$.

2lSeesaw     Hierarchy     2c
2*Linear     2*$\Delta m = \Delta m_\nu$     Normal     $\Delta m_\nu$     $\unit[42.3]{meV}$

Inverted     $\Delta m_\nu$     $\unit[748]{\mu eV}$
3*Inverse     3*$\Delta m = m_\nu \abs{\vec \theta}^{-2}$         $m_\nu$     $\unit[1]{meV}$

$m_\nu$     $\unit[10]{meV}$

$m_\nu$     $\unit[100]{meV}$

[Benchmark models]
Five for the linear and inverse seesaw discussed in the text.
The linear seesaw represent the only two possible models in the minimal linear seesaw of a single pseudo-Dirac heavy neutrino.
With more pseudo-Dirac heavy neutrinos in the spectrum other models become feasible.
Since the inverse seesaw requires at least two pairs of pseudo-Dirac neutrinos it is not possible to uniquely connect it to the measured neutrino mass differences and we use a wide spectrum of neutrino masses as .

[Model dependent sterile neutrino mass splitting]
Heavy neutrino mass splitting as function of the active-sterile mixing parameter for the five of the minimal linear and inverse seesaw given in <ref>.

The mass matrix (<ref>) is the most generic tree level seesaw matrix one can generate from two sterile neutrinos.
It incorporates as limiting cases the linear and the inverse seesaw.
In particular the individual matrix elements have the following phenomenology:

* The Dirac contribution $\vec \mu_D^{}$ yields the linear seesaw for which the neutrino mass matrix is given by $M_\nu = \vec \mu_D^{} \otimes \vec \theta + \vec \theta \otimes \vec \mu_D^{}$ and has the two nonzero eigenvalues $m_\nu = \abs{\vec \mu_D^{}} \abs{\vec \theta} \mp \abs{\vec \mu_D^* \vec \theta}$.
The mass splitting of the heavy neutrinos is in this case $\Delta m = 2 \abs{\vec \mu_D^* \vec \theta}$.
When reproducing the measured light neutrino data in the minimal linear seesaw with a single pseudo-Dirac pair, the mass splitting of the heavy neutrinos is therefore equal to the mass splitting of the light neutrinos $\Delta m = \Delta m_\nu$, <cit.>.

* The Majorana contribution $\mu_M^{}$ in the third entry on the main diagonal yields the inverse seesaw for which the neutrino mass matrix is $M_\nu = \mu_M^{} \vec \theta \otimes \vec \theta$.
Since this expression results in a single finite eigenvalue of $m_\nu = \mu_M^{} \abs{\vec \theta}^2$ at least two pseudo-Dirac pairs are necessary to describe the observed neutrino oscillation data.
The mass splitting for the pseudo-Dirac pair considered here is $\Delta m = \abs{\mu_M^{}}$ and can therefore be expressed in terms of the neutrino mass eigenvalue $\Delta m = m_\nu \abs{\vec \theta}^{-2}$, <cit.>.

* The Majorana contribution $\mu_M^\prime$ in the second entry on the main diagonal introduces a mass splitting between the heavy neutrinos of $\Delta m = \abs{\mu_M^\prime}$ without a tree level impact on the light neutrino masses.
It does, however, generate a contribution to $M_\nu$ at one loop <cit.>.
Furthermore, it can be exploited to reduce the mass splitting necessary to generate realistic neutrino masses when added to the linear or inverse seesaw.
If additionally the linear seesaw term is present the subdominant contribution $\Delta M_\nu = \mu_M^\prime \vec \theta_\epsilon \otimes \vec \theta_\epsilon$ to the neutrino mass matrix is generated <cit.>.

* The $3\times 3$ block in the upper left, which contains the light neutrino masses is zero at tree level.

From the observed neutrino oscillation data it is possible to extract two mass differences which can be combined in two different mass hierarchies, the normal and inverted light neutrino mass ordering.
Whether the third neutrino is also massive can not be determined with current data, hence we use as minimal for the linear seesaw the two possible mass differences between the two massive neutrinos $\Delta m_\nu^{\text{normal}} = \unit[42.27]{meV}$ and $\Delta m_\nu^{\text{inverted}} = \unit[747.9]{\mu eV}$ <cit.>.
At the same time an upper bound on the sum of the neutrino masses $\sum m_\nu < \unit[120]{meV}$ has been established using cosmological observations <cit.>.
Since models relying on the inverse seesaw require at least two pseudo-Dirac pairs, the mass splitting of the heavy neutrinos is not uniquely determined by the observed light neutrino parameters.
However, when defining masses for the light neutrinos as , the heavy neutrino masses can be deduced.
We use as corresponding light neutrino masses of $m_\nu = \unit[1,\,10,\,100]{meV}$.
These are summarised in <ref>.

The relation between the mass splitting and the mixing parameter for the five of the linear and inverse seesaw are visualised in <ref>.
For with a single entry in the mass matrix (<ref>), the minimal value for the mass splitting of $\Delta m = \unit[0.748]{meV}$ corresponding to a maximal oscillation length of $c\tau = \unit[1.04]{cm}$ is reached in the linear seesaw with inverted ordering.
In order to generate pseudo-Dirac pairs with a smaller mass splitting and therefore larger oscillation length one can on the one hand rely on cancellations between the linear and inverse seesaw terms or cancellations involving $\mu_M^\prime$ or on the other hand go beyond the minimal model consisting of a single pseudo-Dirac pair.

Given the small amount of symmetry breaking permissible for consistent low-scale seesaw models, most of their effects can be neglected in phenomenological studies <cit.>.
One important exception are , as they are an interference phenomenon and as such strongly enhanced beyond the typical expectation for an $\order \epsilon$ effect.
As shown in <ref> the at can be described with just two parameters in addition to the symmetric limit of the , which are the mass splitting of the heavy neutrinos $\Delta m$ governing oscillation length in the proper time frame and the damping parameter $\lambda$ governing the effects of decoherence.

Therefore, we introduce the  optimised for phenomenological studies.
Instead of adding the symmetry breaking terms from Lagrangian (<ref>) to the Lagrangian (<ref>) one can directly add the mass splitting without specifying if it originates from a linear seesaw, an inverse seesaw, or a more complicated model.
The masses of the two heavy neutrinos in the are derived from the mass  in the symmetric limit and the single additional parameter $\Delta m$

\begin{equation} \label{eq:Delta m mass splitting}
m_{\nicefrac45}^{} = m_M^{} \left(1 + \frac12 \abs{\vec \theta}^2\right) \mp \frac12 \Delta m \,,
\end{equation}

capturing the entirety of sources contributing to the mass splitting in (<ref>).
Since oscillations are an interference phenomenon they are observable, even though the mass splitting is orders of magnitude below the energies that can be hoped to be resolved in experimental setups.

Additionally the value of the damping parameter $\lambda$ has to be determined for each parameter point.
As shown in <cit.> the effects of decoherence described as observability conditions, and therefore the damping parameter $\lambda$, can be neglected for the part of the parameter space where it is likely that can be reconstructed.
The relevant parameters of the are thus the three active-sterile mixing parameters $\vec \theta$, the heavy neutrino Majorana mass $m_M^{}$, its mass splitting parameters $\Delta m$, and the damping parameter $\lambda$.

§ SOFTWARE IMPLEMENTATION

In order to simulate at colliders we have implemented the in FeynRules.
We provide this implementation online <cit.> and explain its details in <ref>.
Furthermore, we have patched MadGraph as described in <ref> in order to be able to simulate .

§.§ FeynRules model file

3lBLOCK PSPSS #
1    1.000000e+02    # mmaj

2    1.000000e-12    # deltam

3    0.000000e+00    # theta1

4    1.000000e-03    # theta2

5    0.000000e+00    # theta3

6    0.000000e+00    # damping

[FeynRules model file parameter]
The parameter of the implemented in the FeynRules model file as they appear in the MadGraph param.card.

In the following we introduce the quantities used in the FeynRules <cit.> model file <cit.>.
For each quantity we give, besides its mathematical symbol, also its variable name used in the model file denoted in brackets.

According to the discussion in the previous section the Lagrangian has to be extended by two sterile Majorana neutrinos denoted by $N_1$ (N1L) and $N_2$ (N2L).
In the model file an additional parameter for the Higgs denoted by vevNP has been introduced.
The only difference to the usual Higgs (vev) is the interaction order, which has been set to NP, -1 instead of QED, -1 in order to provide MadGraph with the correct power counting rules, but for all physical considerations they can be thought of as equal.
Along with vevNP the field PhiNP is introduced, which corresponds to the usual Higgs doublet Phi where the vev is replaced by vevNP.
The masses of the heavy neutrinos are given by (<ref>) where the Majorana mass parameter is denoted by $m_M^{}$ (Mmaj), the active-sterile mixing parameter given in (<ref>) by $\vec \theta$ (theta1, theta2, theta3), and the mass splitting of the heavy neutrinos is implemented via the parameter $\Delta m$ (deltaM).

The physical fields, the mass eigenstates, are extended by two self-conjugate neutrinos labeled $n_4$ (n4) and $n_5$ (n5).
The mixing matrix relating the mass eigenstate neutrinos to the interaction eigenstates is denoted by $U_n$ (Un) and satisfies

\begin{equation}
U_n^\trans M_n U_n = \diag(0, 0, 0, m_4^{}, m_5^{}) \,,
\end{equation}

up to second order in $\vec \theta$.
The matrix $U_{\CL}^\prime$ (UnCL) is the upper right $3\times5$ part of $U_n$ and is necessary for the automatic index contraction in FeynRules.
A field that contains all mass eigenstate neutrinos is introduced as $n_i$ (nL).
With this the active neutrino interaction eigenstates of the can be rotated to the mass eigenstate neutrinos using the relation

\begin{equation} \label{eq:neutrino rotation}
n_\alpha = (U_{\CL}^\prime)_{\alpha i} n_i \,.
\end{equation}

The neutrinos interact with the only via the lepton doublet $\ell$ (LL).
The mixing is implemented using <ref> by replacing the neutrino $\nu$ (vl) with the transformed neutrino (UnCL nL).

The kinetic terms of the sterile neutrinos

have been added to the LFermions Lagrangian of the .
The mass term

where sp1 is an internal spin index and the Yukawa term

where yvn is the Yukawa vector $\vec y_{\alpha1}$, ff1 is an internal family index and ii and jj are internale $\SU(2)$ doublet indices, have been introduced in the Lagrangian (LNP).
Finally the complete Lagrangian describing interactions of the is denoted by LpSPSS.
Since the implemented mixing matrix $U_n$ is only valid up to second order in $\vec \theta$, the full Lagrangian must be expanded in $\vec \theta$ and terms smaller that $\order{\abs{\vec \theta}^2}$ have to be neglected.
This is achieved using the function RemoveHigherOrder defined in the model file that expands its arguments up to second order in $\vec \theta$ and is automatically applied to the relevant Lagrangian parts.

The only additional free parameters in this model collected in the PSPSS block are the Majorana mass $m_M^{}$ (Mmaj), the mass splitting $\Delta m$ (deltaM), and the three active-sterile mixing parameters $\vec \theta$ (theta1, theta2, theta3).
The damping parameter $\lambda$ (damping) has also been implemented in the FeynRules model file and can be adjusted in the MadGraph param.card.
It is set to zero by default, when it is non-zero, it has the effect of damping the oscillations as discussed <ref>.
All parameters as they appear in the param.card are collected in <ref>.

§.§ Oscillations in MadGraph

[Flowchart of the MadGraph lifetime calculation]
Usual routine of MadGraph assigning the proper lifetime (vtim) to a particle (particle) as long as it larger than the threshold value (threshold) given in the run.card.
The process flows from the node vtim = 0 to the node write event.

[Flowchart of the patched MadGraph lifetime calculation incorporating oscillations]
Patched MadGraph routine assigning a proper lifetime (vtim) in the presence of oscillations to a particle (particle).
The process flows from the node vtim = 0 to the nodes write event and discard event.

The folder exported from FeynRules can be imported it into MadGraph <cit.> using the import model command.
In order to generate the process depicted in <ref> one defines the new multi-particles

As explained below it is necessary to force the heavy neutrino to be on-shell.
This can be achieved by generating a process that contains as well as events via

If the heavy neutrino mass is significantly below the $W$ mass it is also possible to additionally ignore off-shell effects from the $W$ boson by using

A directory with the code for this process can then be obtained using the output command.
In the following the path of this directory is denoted by [pSPSS].

There are two reasons why it is desirable to ensure that the heavy neutrino is produced on-shell.
First, the derivation of the oscillation formulae is based on the Jacob-Sachs theorem, see <cit.>, which takes the heavy neutrinos on-shell for large enough times.
In practice the constraint to be large enough does not give any restriction stronger than the ones already obtained in the observability conditions <cit.>, such that in cases where there are oscillations, the heavy neutrino can be taken to be on-shell.
For this paper it has been assumed that decoherence can be neglected, and therefore effects from times smaller than the threshold are also neglected.
The second reason to force the heavy neutrino to be on-shell is due to the implementation details of the MadGraph patch given here.
Since parameter points that feature oscillations are close to the symmetry limit one expect almost no events for a prompt decay.
Since MadGraph simulates the process as if it where prompt and adds displacement afterwards, one would not get events if interference between different mass eigenstates is taken into account, see also the diagrammatic explanation in <ref>.
Taking the heavy neutrinos on-shell destroys the interference between different mass eigenstates such that and events are created with equal probability.
The correct interference patterns, featuring oscillations and the correct ratio of and events, are then added with the implementation of the patch given here.

The patch has been developed and tested with [2.9.10 (LTS)]MadGraph5_aMC@NLO and is given in <ref>.
The idea behind it necessitates to change the behaviour of the function do_add_Time_of_flight located in the file [pSPSS]/bin/internal/madevent_interface.py,
If one wants to patch MadGraph globally, instead of just patching the individual process, one can apply the patch to the file located in [MadGraph]/madgraph/interface/madevent_interface.py instead.
which derives the decay vertex of a long-lived particle.
The main routine that is executed for each particle in each event and which assigns the lifetime to the particles is depicted in <ref>.
This routine is changed as described in <ref>.
The main idea is to compare the oscillation probability at a given proper time $\tau$, to a pseudo random variable in the range $[0, 1)$.
Based on this it is decided if the heavy neutrino should decay in a or process.
Comparing this with the actual lepton number of the process (leptonnumber), which is zero in the case and different from zero in the case, it is then decided if the event should be kept or discarded.
This algorithm ensures that any spin correlations that are present in the generated events are kept, but the drawback is that half of all events containing heavy neutrinos are discarded.
A possible alternative solution that does not keep the spin correlation but on the other hand does not discard any event, could be to change the charge of the second lepton depending on whether the event should be or .

In order to activate the computation of time of flight one also has to set the parameter time_of_flight in the run.card to an appropriate threshold value (threshold) value, to zero.
Events can then be generated in the usual way using MadEvent.
As described in <ref> we have carefully checked that the generated data faithfully represents the physical processes.

§ EXAMPLE RESULTS

We have calculated various observables within the using the model file and code modification described in <ref>.
For all of the results gathered here we have assumed that the effects of decoherence on the can be neglected, $\lambda = 0$.
In order to demonstrate that we reproduce prior results we show exemplary the heavy neutrino decay width and the expected number of events in a -like detector as a function of the heavy neutrino mass and its coupling strength in <ref>.
Furthermore, we derive the most likely Lorentz factors as well as the fraction of events with an higher Lorentz factor in <ref>.
The main goal of this work is presented in <ref> where we show an example of .
Afterwards, in <ref> the integrated effect is discussed.
In <ref> the dependence of the transverse impact parameter on the angular dependent spin correlation is presented.
We conclude the results in <ref> by commenting on how displaced and prompt searches can be interpreted as bounds on the pseudo-Dirac heavy neutrino of low-scale seesaw models.

§.§ Observable events at an experiment

Decay width in eV.
Number of expected events $N$.

[Heavy neutrino decay width and number of expected events at the $\HLLHC$]
Panel fig:decay width: Scan over the decay width of the heavy neutrino as a function of its mass and the active-sterile mixing parameter $\abs{\vec \theta}^2$.
Panel fig:event number: Scan over the number of expected displaced vertex events in the detector for the $\HLLHC$ with $\mathscr L = \unit[3]{\inv{ab}}$ and cuts as defined in <ref>.

We present a scan over the parameter space resulting in the decay width and the number of expected events for a given experimental setup and luminosity.
For each parameter point the decay width $\Gamma$ and the cross section $\sigma$ is computed and the event obtained from MadGraph is used for a toy analysis.
We present the decay width as a function of the heavy neutrino mass $m$ and the active-sterile mixing parameter $\abs{\vec \theta}^2$ in <ref>.
In order to derive the expected number of events we employ the most important cuts used for displaced vertex analyses at the .
For this example we use cuts inspired by the Phase II detector:

* Minimal transverse momentum $p_T^{\min}(\mu_\text{prompt}) = \unit[20]{GeV}$ of the prompt muon

* Minimal transverse momentum $p_T^{\min}(f_\text{disp}) = \unit[1]{GeV}$ of the displaced muon and quarks

* Maximal pseudorapidity $\eta^{\max} = 4$ for leptons and quarks

* Minimal impact parameter $d_0^{\min}(\mu_\text{disp}) = \unit[2]{mm}$ of the displaced muon

* Maximal distance of half the tracker size for displaced muon and quarks

Where the $p_T^{\min}(\mu_\text{prompt})$ cut is used to ensure that the event is triggered, the $p_T^{\min}(f_\text{disp})$ and $\eta^{\max}$ cuts ensures that the particles are captured by the detector, and the $d_0^{\min}(\mu_\text{disp})$ ensures that the secondary muon is indeed reconstructed as a displaced muon.
Finally, the cut on the maximal displaced distance is used to ensure that the tracks of the displaced muons and of the daughters of the quarks can be measured such that the displaced vertex can be reconstructed.
Using the simulated events and the above mentioned cuts an efficiency factor and the number of expected events is computed using

\begin{align}
N_\text{exp} &= \sigma \mathscr L f_\text{eff} \,, &
f_\text{eff} &= \frac{N_\text{after cuts}}{N_\text{all events}}\,,
\end{align}

The resulting expected number of events as function of the heavy neutrino mass and active-sterile mixing parameter are shown in <ref>.

§.§ Maximal Lorentz factors

Example distribution of Lorentz factors.
Resolvable mass splittings for $L_\text{osc} = \unit[2]{mm}$.

[Lorentz factor dependence of the resolvable mass splitting]
Panel fig:Lorentz factor distribution: Histogram of the simulated Lorentz factor for heavy neutrinos with a mass of [7]GeV.
Panel fig:resolvable mass splittings: Resolvable mass splittings for a fixed oscillation length of $L_\text{osc} = \unit[2]{mm}$ for event samples with different Lorentz factor thresholds $\gamma^f$.
The horizontal lines corresponds to the mass splitting appearing in the two minimal linear seesaw $\BMs$ given in <ref>.

Most likely Lorentz factors as function of mass.
$f$     $\gamma_W^{}$     $\alpha$
$1$     $1$     $0$

$\nicefrac12$     $3.66$     $\expnumber{3.839}{-7}$

$\nicefrac13$     $7.92$     $\expnumber{1.545}{-5}$

$\nicefrac1{10}$     $17.36$     $\expnumber{6.525}{-1}$

Fit values

[Fit of the most likely Lorentz factor thresholds]
Panel fig:most likely Lorentz factors plot: Most likely Lorentz factors $\widehat \gamma$ of heavy neutrinos as function of their masses.
The points correspond to $\MC$ data and the lines to a fit.
Panel fig:tab: fit parameter of the function (<ref>) for the most likely Lorentz factors $\widehat \gamma$ after applying a threshold $f$.

Based on the simulated events, the estimation of the Lorentz factor of the heavy neutrino in <ref> can be tested.
Since this relies on the probability distributions for the Lorentz factor of the heavy neutrinos, one example for such a Lorentz factor distribution is shown in <ref>.
In the following we define the most likely Lorentz factor $\widehat \gamma$ as the one where the distribution has its maximum.

It is also possible to consider only a fraction of events with the highest boost in order to probe the smaller mass splittings that become accessible since the most likely Lorentz factor will be increased in such a data sample.
We define the Lorentz factor threshold $\gamma^f$ such that only the fraction $f$ of events with the highest Lorentz factors is kept.
The mass splittings still resolvable for different thresholds and a fixed oscillation length of $L_\text{osc} = \unit[2]{mm}$ are depicted in <ref>.

In addition to the information in <ref> the boost of the initial $W$ boson can also be taken into account.
To that end the momentum of the heavy neutrino in the rest frame of the $W$ boson is multiplied by a fitting parameter $\alpha\in[0,1]$ estimating the component parallel to the $W$ boson momentum.
Finally, the energy and parallel momentum are boosted with $- \beta_W$ from the $W$ rest frame to the lab frame, yielding

\begin{equation} \label{eq:boost estimate 2}
\gamma \approx
\gamma_W^{} \left(
\frac{m_W^2 + m^2}{2 m_W^{} m}
+ \alpha \frac{m_W^2 - m^2}{2 m_W^{} m} \sqrt{1 - \frac1{\gamma_W^2}}
\right)
\,.
\end{equation}

The parameters $\alpha$ and $\gamma_W^{}$ are obtained by fitting this equation to simulated data and the results are presented in <ref>.
Comparing this most likely Lorentz factor to the estimate of the Lorentz factor <ref> shows that $\gamma_W^{} \approx 1$ as used in <ref> is a good approximation when considering the full event sample as can be seen in <ref>.
The values for the constrained event samples have been used in the derivation of <ref>.

§.§ Oscillations

Lab frame.
Proper time frame.

[Heavy neutrino-antineutrino oscillations]
of the in the lab frame are shown in panel fig:lab frame and the ones in the proper time frame are shown in panel fig:proper time frame.
The parameter point used for these plots features heavy neutrinos with a Majorana mass of $m_M = \unit[20]{GeV}$, an active-sterile mixing of $\abs{\vec \theta}^2 = 10^{-8}$, and a mass splitting of $\Delta m = \unit[748]{\mu eV}$.
The latter corresponds to the linear seesaw point with inverted hierarchy.
The oscillations are simulated by MadGraph using the model file <cit.> presented in <ref> after applying the patch explained in <ref> and explicitly given in <ref>.
While the oscillations in the lab frame are washed out they are clearly visible in the proper time frame.

The main purpose of this work is to study at colliders.
To that goal we use the FeynRules model file of the together with the MadGraph patch presented here to explore if such oscillations are potentially detectable at a like detector.
In <ref> the oscillation pattern as simulated by MadGraph is depicted in the lab frame as well as in the proper time frame.
As discussed in <ref> the oscillation pattern in the lab frame is washed out, since heavy neutrinos with different Lorentz factors overlap.
However, the proper time frame of the heavy neutrino is independent of the Lorentz factor and the oscillation pattern is visible.
The detectability of at the detector during the is studied in detail in <cit.>.

§.§ Integrated effect

[Comparison between simulation and calculation of $R_{ll}$]
Comparison between the analytic calculation (<ref>) and simulated events for the $\LNV$ over $\LNC$ ratio $R_{ll}$ as function of the ratio $\Delta m / \Gamma$.
The three indicated special cases for $R_{ll}(\Delta m/\Gamma)$ correspond to the oscillations shown in <ref>.

$R_{ll}^\text{obs}(R_{ll},\tau_{\min}/\tau_\text{osc})$ while $\tau_{\max}\to\infty$.
$R_{ll}^\text{obs}(R_{ll},\tau_{\max}/\tau_\text{osc})$ while $\tau_{\min}\to0$.

[Simulation of $R_{ll}$ in finite detectors as function of $\tau$]
simulation of the reliability of the observation of $R_{ll}$ when considering finite detectors.
The variables are the same as in <ref>.
The plot in panel Rll obs tau min agrees only in the lower right half with the theoretical results <ref>, since for large cut values and small $R_{ll}$ no events are generated.
The plot in panel Rll obs tau max reproduces the analytic behaviour from <ref>.

$d_0$ cut with $m=\unit[10]{GeV}$
$d_0$ cut with $m=\unit[100]{GeV}$
$d_0$ cut with $m=\unit[1]{TeV}$
$d_\text{min}$ cut with $m=\unit[10]{GeV}$
$d_\text{min}$ cut with $m=\unit[100]{GeV}$
$d_\text{min}$ cut with $m=\unit[1]{TeV}$
$d_\text{max}$ cut with $m=\unit[10]{GeV}$
$d_\text{max}$ cut with $m=\unit[100]{GeV}$
$d_\text{max}$ cut with $m=\unit[1]{TeV}$

[Simulation of $R_{ll}$ in finite detectors as function of $d_0$, and $d$]
Simulation of $R_{ll}$ under consideration of a finite $d_0$ in Rll obs d0 10 GeV, Rll obs d0 100 GeV, and Rll obs d0 1 TeV, of a finite $d_{\min}$ in Rll obs dmin 10 GeV, Rll obs dmin 100 GeV, and Rll obs dmin 1 TeV as well as of a finite $d^{\min}$ in Rll obs dmax 10 GeV, Rll obs dmax 100 GeV, and Rll obs dmax 1 TeV.
Each cut is shown for three different masses $m= \unit[10,\,100,\,1000]{GeV}$.

In the cases where oscillations can not be directly resolved there might still be non trivial relations between and events.
The integrated effect presented in <ref> depends on the interplay between the decay width and the oscillation length as shown in <ref>.
Counting the number of opposite- and same-sign events while scanning the parameter space over different values for the Majorana mass $m_M^{}$ and the mass splitting $\Delta m$ reproduces the analytic dependence of (<ref>) as shown in <ref> <cit.>.

Furthermore, the simulation reproduces the analytic calculation leading to <ref> when taking the finite detector size into account as shown in <ref>.
The simulation agrees with the analytic calculation for $\tau_{\min} \to 0$.
However, for $\tau_{\max} \to 0$ the upper left part of the plot are inaccessible no event is generated for these parameter points.
In particular this effects regularises the divergent parts in the corresponding analytic calculation.
Although the cut in $\tau$ corresponds to the cleanest theoretical description of $R_{ll}$ with finite detector effects, it does not correspond to realistic physical cuts.
Therefore, we present additional result using a $d_0$ cut, a minimal distance lab frame cut $d_{\min}$, and a maximal distance lab frame cut $d_{\max}$ in <ref>.
Since, the boost to the lab frame introduces a dependence on the Lorentz factor we report the results for the three mass points $m = \unit[10,\,100,\,1000]{GeV}$.

For the masses above the $W$ threshold the $d_0$ cut reproduces approximately the first maximum of the analytic calculation.
Furthermore, the area in which no events are generated becomes larger and for smaller masses the first maximum is pushed to larger values of the oscillations length.
When applying a minimal distance cut in the lab frame $d_{\min}$, parts of the behaviour described for a $\tau_{\min}$ and $d_0$ cut are washed out and the vertical lines that indicate trivial relation between $R_{ll}$ and $R_{ll}^\text{obs}$ are recovered.
However they are shifted and compressed towards lower values of $R_{ll}$.
Nonetheless, the maximum remains although it is shifted to larger values of the oscillation length.
A maximal distance cut in the lab frame $d_{\max}$ does not show the oscillatory patter and leads only to small distortions for cuts larger than the oscillations length.
However, for cuts below one to two oscillations lengths large deviations can be observed.

We conclude, that although the Lorentz factor distribution smears the finite detector effect observed <ref> it remains important and can not be neglected, especially, when applying a $d_0$ cut of order of the oscillations length.

§.§ Transverse impact parameter and spin correlation

[Angular dependence of heavy neutrino decays]
Comparison between the angular dependence of the process compared to the process in the observable (<ref>) using a simulation.

The transverse impact parameter $d_0$ of particles emerging in a secondary vertex is defined as

\begin{equation} \label{eq:impact parameter}
d_0 =
\frac{\vec d_T^\prime \wedge \vec p_T^\prime}{p_T^\prime} =
\frac{\epsilon_{ij} x_i^\prime p_j^\prime}{p_T^\prime} =
\frac{x_{}^\prime p_y^\prime - y_{}^\prime p_x^\prime}{p_T^\prime}
\,,
\end{equation}

where $\vec d_T^\prime = \row{x^\prime, y^\prime}$ is the position of the point with the smallest distance to the $z$-axis and $\vec p_T^\prime = \row{p_x^\prime, p_y^\prime}$ is the transverse momentum at this point.
[The two dimensional wedge product results in a scalar quantity.]
In cases where the transverse momentum of the particle is large or the magnetic field is small, such that the radius of the trajectory of the particle is much larger than the relevant length scales of the detector, the trajectory is well approximated by a straight line.
In this case the value of $d_0$ is constant for every point on the trajectory.
It is therefore possible to replace $\vec d_T^\prime$ with the coordinates of the production vertex $\vec d_T = (x,y)$ which coincides with the transverse distance of the production vertex.
In the case of the displaced muon, whose production vertex is the decay vertex of the heavy neutrino, one can furthermore use the relation

\begin{equation}
\vec d_T^\prime \rightarrow \vec d_T = d_T \frac{\vec p_T^N}{p_T^N}\,,
\end{equation}

where $\vec p_T^N$ is the transverse momentum of the heavy neutrino.
This yields a simplified definition of the impact parameter for displaced muons, under the assumption that their trajectories can be approximated as straight lines

\begin{align} \label{eq:impact parameter simplified}
d_0 &\simeq d_T \frac{\vec p_T^N\! \wedge \vec p_T^\mu}{p_T^N p_T^\mu} = d_T \sin \phi(N,\mu)\,, &
\sin \phi(N,\mu) &= \frac{\vec p_T^N\! \wedge \vec p_T^\mu}{p_T^N p_T^\mu} \,.
\end{align}

Therefore, the impact parameter consists, besides the transverse decay length of the invisible particle, of an angular depend part.
The angular dependent part encodes a spin dependence which can be exploited to distinguish the and processes as shown in <ref>.
One can see that just by measuring this or a related quantity such as the angle between the two charged leptons in <ref> it is possible to compute probabilities for the event being or .
On the other hand, imposing a minimal $d_0$ cut is the standard strategy to reduce the background of displaced searches as seen in <ref>.
The angular dependence of $d_0$ leads after cuts to a residual oscillation in the observable $N_{\LNC}+N_{\LNV}$ that naively would be expected to be free of oscillations.

§.§ Experimental bounds on low-scale seesaw models

[Model dependent comparison between decay length and $R_{ll}$]
Contour lines for the heavy neutrino decay length $c\tau$ in m overlayed by contour bands of $R_{ll} \in [0.1,0.9]$ for the summarised in <ref>.
While the two bands generated by the linear seesaw are well separated and reach far below the $W$ mass threshold the bands generated by the inverse seesaw are very close to each other and only reach below the $W$ threshold for comparably small values of the active-sterile mixing parameter.
For comparison the regions excluded by the and experiments via displaced searches are shown as grey areas <cit.>.
The areas shaded in grey at the center top are excluded by prompt searches <cit.>, as long as the prompt same-sign dilepton signature is valid.
This is the case only for seesaw models in which $R_{ll}$ is close to one in this region, hence to the bottom left of the $R_{ll}$ bands of the under consideration.
The most optimistic reach of the $\HLLHC$ <cit.>, -$ee$ <cit.>, and -$he$ <cit.> are shown as black lines.
When producing the figure we have neglected potential effects from a non-zero damping factor $\lambda$ as well as from times shorter than the Jacob-Sachs threshold, for details see <ref>.

Although neither pure Dirac nor single Majorana particles are sufficient to describe the phenomenology of low-scale seesaw heavy neutrinos, almost all prior searches are performed using these two models.
Therefore it is necessary to correctly interpret these searches in order to extract the proper bounds on the pseudo-Dirac heavy neutrinos of low-scale seesaw models.
Customarily, the searches relying on displaced vertex signatures are reported for both, pure Dirac and single Majorana .
Since these searches usually do not rely on , the major difference affecting these searches is the factor of two that appears in the decay width and is described in detail in <ref>.
Since the Dirac particle has the same decay width as the pseudo-Dirac particle it is prudent to interpret the exclusion extracted for a displaced Dirac as the correct exclusions for displaced pseudo-Dirac heavy neutrinos.
The situation is more complicated in the case of prompt searches.
The majority of prompt searches rely on the very clean same-sign dilepton signature of decays.
Since these are not generated for Dirac particles the searches are only reported for single Majorana particles.
Therefore, these searches assume an $R_{ll}$ of one, which depending on the mass splitting might not be a good approximation of reality.
Hence the exclusion bound depends on the details of the model under consideration.
In practice one has to check the intersection between the $R_{ll}$ band of one specific model and the reported exclusion bounds for the same-sign dilepton search in a representation such as the one presented in <ref>.
For $R_{ll}$ close to one (to the bottom left of the bands) the reported exclusion bounds are valid up to the factor of two between the simulated single Majorana particle and the realistic pseudo-Dirac particle discussed in <ref>.
For $R_{ll}$ close to zero (to the top right of the bands) there will be no , rendering the searches inconsequential.
Therefore, prompt searches can only be interpreted model-dependently, as a function of $\Delta m$.
The model file presented in <ref> together with the example implementation of in <ref> provide a unified framework for the correct simulation of prompt and displaced processes of pseudo-Dirac heavy neutrinos in low-scale seesaw models.

§ CONCLUSION

In this paper we have reviewed the arguments that collider testable seesaw models, that do not rely on tuning, must be protected by an approximate .
The breaking of the symmetry by small parameters simultaneously ensures small , tiny neutrino masses, and generically generates a small mass splitting between the heavy neutrinos.
The latter leads to that can have macroscopic oscillation lengths such that they are potentially resolvable at high energy colliders.

In order to simulate the effects of these oscillations we have reviewed the that builds on a systematic expansion in the small parameter governing the breaking of the .
We have then reduced it to the that relies on the minimal number of parameters necessary to describe pseudo-Dirac heavy neutrinos and in phenomenological studies.
In addition to the parameters, the introduces only the Majorana mass parameter $m_M$, the active-sterile mixing parameters $\vec \theta = \row{\theta_e,\theta_\mu,\theta_\tau}^\trans$, the mass splitting $\Delta m$ governing the oscillation length, and a damping parameter $\lambda$ that captures possible decoherence effects.
The last of which can be neglected in the parts of the parameter space where prospects are good to resolve the oscillations, as discussed in <cit.> and further elaborated in <cit.>.

In order to facilitate studies using this model, we have published a FeynRules implementation of the together with a patch extending MadGraph in such a way that it is capable of simulating .
Using these tools we have calculated some example results:
We have demonstrated that our implementation recovers typical prior results such as the heavy neutrino decay width.
Additionally, we have provided the maximal mass splitting $\Delta m$ for which it may be feasible to resolve the oscillation pattern, when only a fraction of events containing the largest Lorentz factors is considered.
Furthermore, we have shown that the integrated effect of the oscillations is reproduced by our model file.
We have emphasized in this context that care has to be taken when measuring $R_{ll}$, as the finite detector geometry can have a major impact on the derived values.
Finally, we have demonstrated the dependence of the impact parameter on the angular dependent spin correlation, and have shown that this observable can in principle be used to distinguish between and decays.

§.§ Acknowledgements

The work of J.H. was partially supported by Portuguese Fundação para a Ciência e a Tecnologia (FCT) through the projects CFTP-FCT unit UIDB/00777/2020 and UIDP/00777/2020.
S.A. and J.R. acknowledge partial support from the Swiss National Science Foundation grant 200020/175502.

§ DECAY WIDTHS OF MAJORANA, DIRAC, AND PSEUDO-DIRAC PARTICLES

In the main part of this article we have discussed that in order to explain light neutrino masses the corresponding heavy neutrino can not be a Dirac particle.
Additionally, in order to generate the observed light neutrino oscillation pattern at least two neutrinos have to be added to the .
Nevertheless, it is standard practice to study the phenomenology of a single heavy neutrino that is either of Majorana or Dirac type.
In the following we compare their phenomenology with the one of pseudo-Dirac heavy neutrinos with different mass splittings and comment on the decay widths of each of these particles.

Single Majorana

The relevant part of the Lagrangian describing a single Majorana particle $N_1$ reads

\begin{equation}
\mathcal L_\text{Majorana} = y_\alpha \widebar L_\alpha N_1^c  H + \frac{1}{2} m_M^{} \widebar N_1 N_1^c + \HC\,.
\end{equation}

After the coupling to the $W$ bosons is obtained from the kinetic terms of the active neutrinos via mixing of light $\nu_\alpha$ and heavy $n_4$ neutrino mass eigenstates and reads

\begin{align}\label{eq:coupling_majorana_W}
\mathcal L_\text{Majorana}^W &= \theta_\alpha' \widebar l_\alpha n_4 W  + \HC \,, &
& \coupling{n_4}{l_\alpha}{\theta_\alpha^\prime}\,.
\end{align}

where the mixing parameter is given by $\vec \theta' = \vec y v m_M^{-1}$.
The expected number of and events at a collider experiment can then be schematically computed using

\begin{equation}
\abs{\mathcal A}^2 = \abs*{\scatterLNC{\theta_\alpha'}{m}}^2 + \abs*{\scatterLNV{\theta_\alpha'}{m}}^2 \propto 2 \abs{\theta_\alpha^\prime}^2 \,,
\end{equation}

where the diagram encodes the production and decay of a heavy neutrino from and to charged leptons and $W$ bosons.
In this case one expects an equal number of and events and therefore $R_{ll}=1$.
As long as $\Gamma \ll m$, the magnitude of the (greyed out) coupling in the decay vertex is not relevant since it is cancelled by the decay width dependence of the propagator of the heavy neutrino, which can be verified using the narrow width approximation.

In general the value of the active-sterile mixing parameter can be obtained from measuring the production cross section and the decay width of the heavy neutrino.
In the case of a single Majorana particle the production cross section is proportional to $\abs{\vec \theta^\prime}^2$.
Generically, the decay width can be obtained from scanning the $p^2$ distribution of any process with the heavy neutrino mass eigenstate on the $s$-channel and compare it to the Breit-Wigner shape.
Furthermore, for long-lived particles it can be extracted from the decay exponential reconstructed by counting the number of displaced vertices per distance interval.
In the case of a single Majorana particle the decay width is proportional to $\abs{\vec \theta^\prime}^2$.

A heavy Dirac neutrino can be described as two mass degenerate Majorana neutrinos $N_{\nicefrac12}$ with an exact phase difference of $-i$ in their couplings.
The relevant part of the Lagrangian is given by

\begin{equation}
\mathcal L_\text{Dirac} = y_\alpha \widebar L_\alpha N_1^c H + m_M^{} \widebar N_1 N_2^c + \HC\,.
\end{equation}

After the mass matrix of the neutrinos can be diagonalised and the active neutrino interaction eigenstates can be expressed as a linear combination of the heavy neutrino mass eigenstates.
These linear combinations are the heavy neutrino $N$ and heavy antineutrino $\widebar N$ introduced in <ref>

\begin{align}
N &= \frac{\theta}{\sqrt2} \left(-i n_4 + n_5 \right) \,, &
\widebar N &= \frac{\theta}{\sqrt2} \left(i n_4 + n_5 \right)\,,
\end{align}

which leads to the coupling to the gauge bosons via

\begin{align} \label{eq:Dirac mass eigenstate}
\mathcal L_\text{Dirac}^W &= \theta_\alpha \widebar l_\alpha \frac{1}{\sqrt2} \left(- i n_4 + n_5 \right) W  + \HC \,, &
& \coupling{n_4}{l_\alpha}{\pm i \frac{\theta_\alpha}{\sqrt2}} \,, &
\coupling{n_5}{l_\alpha}{\frac{\theta_\alpha}{\sqrt2}}\,,
\end{align}

where $\vec \theta = \vec y v m_M^{-1}$.
Compared to the case of a single Majorana particle the coupling of the gauge bosons to each heavy neutrino mass eigenstate obtains an additional factor of $\flatfrac{1}{\sqrt2}$, which originates from the diagonalisation of the heavy neutrino $2\times2$ block in the mass matrix.
The number of and events can be derived using the Feynman diagrams

\begin{equation}\label{eq:feynman_dirac}
\begin{aligned}
\abs{\mathcal A}^2 &=
\abs*{\scatterLNC[\frac{(\pm i\theta_\alpha)^*}{\sqrt 2}]{\frac{\pm i\theta_\alpha}{\sqrt 2}}{m_4} + \scatterLNC[\frac{\theta_\alpha^*}{\sqrt 2}]{\frac{\theta_\alpha}{\sqrt 2}}{m_5}}^2 +
\abs*{\scatterLNV{\frac{\pm i\theta_\alpha}{\sqrt 2}}{m_4} + \scatterLNV{\frac{\theta_\alpha}{\sqrt 2}}{m_5}}^2 \\
& = \abs*{\scatterLNC[\frac{(\pm i\theta_\alpha)^*}{\sqrt 2}]{\pm i\frac{\theta_\alpha}{\sqrt 2}}{m_4} + \scatterLNC[\frac{\theta_\alpha^*}{\sqrt 2}]{\frac{\theta_\alpha}{\sqrt 2}}{m_5}}^2 = 4 \abs*{\scatterLNC[\frac{\theta_\alpha^*}{\sqrt 2}]{\frac{\theta_\alpha}{\sqrt 2}}{m_{\nicefrac45}}}^2 \propto 4 \abs*{\frac{\theta_\alpha}{\sqrt2}}^2 = 2 \abs{\theta_\alpha}^2 \,.
\end{aligned}
\end{equation}

where the Dirac condition ensures mass degeneracy $m_4 = m_5$.
In contrast to the Majorana heavy neutrino the Dirac heavy neutrino generates only decays and therefore $R_{ll} = 0$.
In the picture where the heavy Dirac neutrino is represented by two Majorana this is due to the phase difference of $-i$ between the couplings of the mass eigenstates, ensuring that the diagrams cancel exactly.
Both, the Dirac and the Majorana heavy neutrino, yield the same total number of events, such that it is not possible to distinguish them just by the observation of the overall number of events.

The production cross section of a single mass eigenstate $n_{\nicefrac45}$ in (<ref>) is proportional to $\abs{\vec \theta}^2/2$.
However, since it is intrinsically impossible to distinguish the two identical mass eigenstates in a single Dirac particle the measured production cross section is proportional to $\abs{\vec \theta}^2$.
Nonetheless, the decay width is only proportional to $\abs{\vec \theta}^2/2$ as one does not sum over the two in order to calculate the total decay width.
Therefore, in comparison to the single Majorana heavy neutrino the decay width is reduced by a factor of two.

Naively one might expect to distinguish Dirac from Majorana heavy neutrinos in this way.
However, this is impossible as we will demonstrate using different limits of a model consisting of a pair of two heavy Majorana neutrinos with finite mass splitting.

Majorana pair

The Majorana and Dirac heavy neutrinos can both be described as limiting cases of a model with two Majorana particles encompassing also the pseudo-Dirac heavy neutrino.
Starting from the Dirac heavy neutrino a small perturbation of the , such that the two masses of the Majorana as well as the modulus of their couplings are no longer exactly degenerate, generates a pseudo-Dirac heavy neutrino, as explicitly shown in <ref>.
Conversely, the Dirac heavy neutrino can be seen as the limit of the pseudo-Dirac heavy neutrino in which the perturbation of the goes to zero, such that the symmetry is restored.
In the limit of large mass splitting the phenomenology of two well separated Majorana heavy neutrinos is recovered.
In the following, when we refer to the size of the mass splitting we imply that the Yukawa couplings are also perturbed from the degeneracy given in the Dirac case.

For small mass splittings as described in <ref> the mass eigenstates $n_{\nicefrac45}$ in (<ref>) are still produced coherently describing the phenomenology of a pseudo-Dirac pair.
However, the diagrams do not cancel exactly anymore, leading to processes.
While the total number of produced events stays the same, are possible such that a part of generated events are now leading to $0< R_{ll} < 1$.
However, the arguments for the production and decay of a pure Dirac particle apply also to this case.

Intermediate mass splitting

While intermediate mass splitting can still be too small to experimentally distinguish different mass eigenstates it is still possible that decoherence already sets in.
This effect is governed by the observability conditions <cit.> and by the damping factor included in the patch.
In this case the interference between diagrams with different mass eigenstates in (<ref>) is destroyed.
Therefore the and processes become both equally likely while the total number of events is proportional to

\begin{multline}
\abs{\mathcal A}^2 =
\abs*{\scatterLNC[\frac{(\pm i\theta_\alpha)^*}{\sqrt 2}]{\frac{\pm i\theta_\alpha}{\sqrt 2}}{m_4}}^2 + \abs*{\scatterLNC[\frac{\theta_\alpha}{\sqrt 2}]{\frac{\theta_\alpha}{\sqrt 2}}{m_5}}^2 +
\abs*{\scatterLNV{\frac{\pm i\theta_\alpha}{\sqrt 2}}{m_4}}^2 + \abs*{\scatterLNV{\frac{\theta_\alpha}{\sqrt 2}}{m_5}}^2 \\
\propto 4 \abs*{\frac{\theta_\alpha}{\sqrt2}}^2 = 2 \abs{\theta_\alpha}^2
\end{multline}

and thus remains the same.
The production cross section corresponds to the sum of the production cross sections of both mass eigenstates, it is proportional to $\abs{\vec\theta}^2$, while the decay width is the one of each individual mass eigenstate and therefore proportional to $\abs{\vec \theta}^2/2$.
Although the ratio between and processes is identical to the case of a single heavy Majorana neutrino, one still obtains a decay width that is a factor two smaller compared to the case of a single Majorana particle.
Therefore, from the comparison of the decay width and the production cross section of a newly discovered heavy neutrino one can gain information about the number of Majorana taking part in the production cross section.
However, it is not possible to distinguish a Dirac particle from two Majorana particles with a small or intermediate mass splitting
just by the comparison of the production cross sections, decay widths or total event numbers.

Large mass splitting

For large mass splitting the individual mass eigenstates can be experimentally distinguished.
As before one finds for each mass eigenstate the production cross section as well as the decay width to be approximately proportional to $\abs{\vec \theta}^2/2$.
Considering that it is likely that the lighter of the two mass eigenstates is discovered first it is insightful to express the coupling in the language of a single Majorana, which would lead to $\abs{\vec \theta^\prime}^2 =: \abs{\vec \theta}^2/2$.
Therefore, one obtains exactly the same phenomenology as if the two mass eigenstates are two distinct Majorana particles, avoiding the unnecessary interpretation as a quasi pseudo-Dirac particle with large mass splitting.

Pure Dirac     3cMajorana Pair     Single Majorana
Mass splitting         Small     Intermediate     Large
Production     $\abs{\vec \theta}^2\hphantom{/2}$     $\abs{\vec \theta}^2\hphantom{/2}$     $\abs{\vec \theta}^2\hphantom{/2}$     $\abs{\vec \theta^\prime}^2$     $\abs{\vec \theta^\prime}^2$

Decay     $\abs{\vec \theta}^2/2$     $\abs{\vec \theta}^2/2$     $\abs{\vec \theta}^2/2$     $\abs{\vec \theta^\prime}^2$     $\abs{\vec \theta^\prime}^2$

$R_{ll}$     $0$     $0<R_{ll}<1$     $1$     $1$     $1$

[Production ratios and decay widths for different types of heavy neutrinos]
Production ratios and decay widths for different types of heavy neutrinos.
The primed mixing parameters occur in models without mass mixing between two heavy states, while the unprimed mixing parameter occur in models with mass mixing between the two heavy neutrinos.

We summarise the discussion above in <ref> and conclude be noting that the comparison of the decay widths yields insight to whether there are two Majorana , as in the case of Dirac particles, or only one, as in the case of a single Majorana particle, but it can not prove the Dirac property of the heavy neutrino.
This is already clear from the fact that the factor of two in the decay widths approximately remains when a small term is added which turns the pure Dirac heavy neutrino into a pseudo-Dirac heavy neutrino.
To test the Dirac property of the heavy neutrino it is unavoidable to test for .
From a measurement of the mass(es) and lifetime(s) one can at best conclude that there are two Majorana particles with almost identical properties.
But whether these two can be combined to form a Dirac heavy neutrino depends on the phase between the Yukawa couplings of the two mass degenerate Majorana particles.
Only when the phase is $-i$, then on the one hand cancels out and on the other hand the two Majorana particle can be combined into a Dirac heavy neutrino.
As long as this phase (related to ) is not tested, one can only speak of two Majorana .

float,captionpos=b,caption=[Original MadGraph code]
Original MadGraph code calculating particle lifetime.
for event in lhe:
for particle in event:
id = particle.pid
width = param_card['decay'].get((abs(id),)).value
if width:
vtim = c * random.expovariate(width/cst)
if vtim > threshold:
particle.vtim = vtim
#write this modify event

caption=[Patched MadGraph code]
Patched MadGraph code calculating particle lifetime and .
Note that the sterile neutrino particle IDs and the position of the  parameter in the param.card are model file dependent and might need to be adjusted.

, label=lst:patch
mass_splitting = param_card.get_value('PSPSS', 2)||
damping = param_card.get_value('PSPSS', 6)||
for event in lhe:
leptonnumber = 0
write_event = True
for particle in event:||
if particle.status == 1:
if particle.pid in [11, 13, 15]:
leptonnumber += 1
elif particle.pid in [-11, -13, -15]:
leptonnumber -= 1
for particle in event:
id = particle.pid
width = param_card['decay'].get((abs(id),)).value
if width:
if id in [8000011, 8000012]:||
tau0 = random.expovariate(width / cst)||
if 0.5 * (1 + math.exp(-damping)*math.cos(mass_splitting * tau0 / cst)) >= random.random():||
write_event = (leptonnumber == 0)
write_event = (leptonnumber != 0)
vtim = tau0 * c||
vtim = c * random.expovariate(width / cst)
if vtim > threshold:
particle.vtim = vtim
# write this modify event
if write_event:||

§ MADGRAPH PATCH

The patch presented here is used to include in MadGraph.
It is applied to the [pSPSS]/bin/internal/madevent_interface.py
In other versions the relevant section might be found in the [pSPSS]/bin/internal/common_run_interface.py file.
generated by [2.9.10 (LTS)]MadGraph5_aMC@NLO.
To apply the patch the main event loop, responsible for adding time of flight information to the .lhe file given in <ref> has to be replaced by the code presented in <ref>.
The relevant lines can be found by searching for the keyword vtim.
The loop over the particles in the event, starting in <ref> is used to determine the leptonnumber of the event.
Eletrons, muons and $\tau$-leptons are counted as a lepton, whereas their antiparticles are counted as antileptons.
The masssplitting is extracted from the PSPSS block in the param.card in <ref> and the damping factor damping similar in <ref>.
In <ref> the program executes the new code of the patch if the particle is a heavy neutrino with particle ID 8000011 or 8000012 as defined in the model file presented in <ref>.
First the proper time at which the heavy neutrino decays is obtained in <ref>.
Next the formula for the oscillation probability in the proper time frame is used together with a random number to decide if the event should be or .
Subsequently if the leptonnumber agrees with the result from <ref> the switch write_event is set to True.
Otherwise it is set to False.
In <ref> the proper time is converted to mm as in the unpatched code from MadGraph.
Finally using the write_event switch from above it is decided if the event should be kept or discarded in <ref>.

Note that when using the automated width computation of MadGraph for some parameter points there are issues with the computation of the phase space volume in the decay $n_5 \rightarrow n_4 \nu \nu$, where $n_5$ is the heaviest neutrino and the $\nu$ are massless neutrinos, due to the unusually small mass splitting.
Since we found this decay channel to be negligible, we fixed the issue by replacing the argument of the square root of the return value of the function calculate_apx_psarea in the file [MadGraph]/mg5decay/decay_objects.py by its absolute value.
However, this can be a drastic change  of the MadGraph behaviour, therefore we have checked that the results are physically meaningful <cit.>.
In our case that was found to be true, since the problematic decay channel is discarded by MadGraph at a later step independently of the change, so there is no impact on any physical results.
Additionally, we also compared the automated computed width for a parameter point with a small mass splitting of the heavy neutrinos to a parameter point with a large mass splitting, changing the other parameters in such way that the mass of $n_4$ as well as its couplings constant is stays constant.
The two computations yield the same results, showing that there is no impact from the extremely small mass splitting in the computation performed by MadGraph.

Missing 'biblatex' package
The bibliography requires the 'biblatex' package.

journaltitlePhys. Rev. Lett.
titleDirect evidence for neutrino flavor transformation from neutral current interactions in the Sudbury Neutrino Observatory

journaltitlePhys. Lett. B
title$\mu \to e\gamma$ at a Rate of One Out of $10^9$ Muon Decays?
reportnumberPrint-77-0182 (BERN)

journaltitleConf. Proc. C
titleComplex Spinors and Unified Theories

journaltitlePhys. Rev. D
titleNeutrino Masses in SU(2)${}\times{}$U(1) Theories

journaltitlePhys. Rev. Lett.
titleNeutrino Mass and Spontaneous Parity Nonconservation

journaltitleProg. Theor. Phys.
titleHorizontal Symmetry and Masses of Neutrinos

journaltitlePhys. Rev. D
titleNeutrino Decay and Spontaneous Violation of Lepton Number

journaltitleNuovo Cim.
titleTeoria simmetrica dell’elettrone e del positrone

journaltitlePhys. Rev. D
titleNeutrinoless Double beta Decay in SU(2)${}\times{}$U(1) Theories

journaltitlePhys. Rev.
titleOn transition probabilities in double beta-disintegration

journaltitlePhys. Rev. D
titleRight-Handed Neutrinos at CERN LHC and the Mechanism of Neutrino Mass Generation

titleHeavy Neutrino-Antineutrino Oscillations in Quantum Field Theory

journaltitlePhys. Rev. D
titleOscillation of heavy sterile neutrino in decay of $B \to \mu e \pi$

titleQuasi-Dirac neutrinos at the LHC

journaltitleMod. Phys. Lett. A
titleResolvable heavy neutrino–antineutrino oscillations at colliders

titleTesting sterile neutrino extensions of the Standard Model at future lepton colliders

journaltitleInt. J. Mod. Phys. A
titleSterile neutrino searches at future $e^-e^+$, $pp$, and $e^-p$ colliders

journaltitlePhys. Rev.
titleBehavior of neutral particles under charge conjugation

journaltitleSov. Phys. JETP
titleMesonium and anti-mesonium

journaltitleSov. Phys. JETP
titleNeutrino Experiments and the Problem of Conservation of Leptonic Charge

journaltitlePhys. Rev.
titleObservation of Long-Lived Neutral V Particles

journaltitlePhys. Lett. B
titleObservation of $B^0-\bar B^0$ Mixing

journaltitlePhys. Rev. Lett.
titleObservation of $B^0_s - \bar{B}^0_s$ Oscillations

journaltitlePhys. Rev. Lett.
titleObservation of $D^0 - \bar D^0$ oscillations

journaltitlePhys. Rept.
titleOscillations of neutrinos and mesons in quantum field theory

notein preparation

titleNeutrino wave packets in quantum field theory

journaltitlePhys. Rev.
titleMass and Lifetime of Unstable Particles

titleNon-unitary Leptonic Mixing and Leptogenesis

journaltitleZ. Phys. C
titleRadiatively induced neutrino masses and large Higgs neutrino couplings in the standard model with Majorana fields

journaltitlePhys. Rev. D
titleElectroweak-scale resonant leptogenesis

titleRadiative Corrections to Light Neutrino Masses in Low Scale Type I Seesaw Scenarios and Neutrinoless Double Beta Decay

titleThe fate of hints: updated global analysis of three-flavor neutrino oscillations

journaltitleAstron. Astrophys.
titlePlanck 2018 results. VI. Cosmological parameters

subtitlePhenomenological symmetry protected seesaw scenario
titleaddonFeynRules model file

journaltitleComput. Phys. Commun.
subtitleA complete toolbox for tree-level phenomenology
titleFeynRules 2.0

subtitleGoing Beyond
titleMadGraph 5

notein preparation
titleHeavy neutrino-antineutrino oscillations at the HL-LHC

titleOn lepton number violation in heavy neutrino decays at colliders

titleHNL mass degeneracy: implications for low-scale seesaws, LNV at colliders and leptogenesis

journaltitlePhys. Rev. D
titleSearch for long-lived, massive particles in events with a displaced vertex and a muon with large impact parameter in $pp$ collisions at $\sqrt{s} = 13$ TeV with the ATLAS detector

titleSearch for long-lived heavy neutral leptons with displaced vertices in proton-proton collisions at $ \sqrt{\mathrm{s}} = 13$ TeV

titleSearch for heavy Majorana neutrinos with the ATLAS detector in pp collisions at $ \sqrt{s} = 8$ TeV

journaltitlePhys. Rev. Lett.
titleSearch for heavy neutral leptons in events with three charged leptons in proton-proton collisions at $\sqrt{s} = 13$ TeV

titleHeavy Neutrinos in displaced vertex searches at the LHC and HL-LHC

journaltitleEur. Phys. J. C
titleHECATE: A long-lived particle detector concept for the FCC-ee or CEPC

titleSearches for Long-Lived Particles at the Future FCC-ee

titleLepton-Trijet and Displaced Vertex Searches for Heavy Neutrinos at Future Electron-Proton Colliders

titleAutomatic Width Computation with Small Mass Gap

journaltitleZh. Eksp. Teor. Fiz.

journaltitleZh. Eksp. Teor. Fiz.
journaltitleAstron. Astrophys.