Beyond lepton number violation at the HL-LHC: Resolving heavy neutrino-antineutrino oscillations

Collider testable low-scale seesaw models predict pseudo-Dirac heavy neutrinos, that can produce an oscillating pattern of lepton number conserving and lepton number violating events. We explore if such heavy neutrino-antineutrino oscillations can be resolved at the HL-LHC. To that end, we employ the first ever full Monte Carlo simulation of the oscillations, for several example benchmark points, and show under which conditions the CMS experiment is able to discover them. The workflow builds on a FeynRules model file for the phenomenological symmetry protected seesaw scenario (pSPSS) and a patched version of MadGraph , able to simulate heavy neutrino-antineutrino oscillations. We use the fast detector simulation Delphes and present a statistical analysis capable of inferring the significance of oscillations in the simulated data. Our results demonstrate that, for heavy neutrino mass splittings smaller than about 100 $\mu$eV, the discovery prospects for heavy neutrino-antineutrino oscillations at the HL-LHC are promising.


Introduction
The observed light neutrino flavour oscillations [1] can be explained by introducing at least two mass splittings between the light neutrinos.It follows that at least two of the neutrinos have to be massive.However, the Standard Model (SM) of particle physics lacks the corresponding mass terms.Therefore, it needs to be extended with a model able to generate these light neutrino masses, such as the type I seesaw mechanism [2][3][4][5][6][7].In order to generate the two light neutrino masses in this model at least two sterile neutrinos have to be added to the particle content of the SM.In general, the type I seesaw mechanism is governed by two sets of parameters, that can be adjusted to obtain light neutrino masses at the right order of magnitude.Besides their Yukawa couplings, the second set is given by their Majorana masses m M , which are only possible since the additional neutrinos are SM singlets.After electroweak symmetry breaking (EWSB) the Yukawa couplings result in Dirac mass terms m (i) D between each of the sterile neutrinos N i and the active neutrinos ν α . 1 For the case of two heavy neutrinos, the generated light neutrino mass matrix then follows from the seesaw equation M . (1.1) In order to obtain light neutrino masses, one could take the Yukawa couplings to be very small or choose the Majorana masses of the heavy neutrinos to be very large.Those two limits of the seesaw mechanism are referred to as the small coupling and the high scale limit, respectively.In both cases the parameters have to be taken to such extreme values that a direct detection of the sterile neutrinos at current collider experiments is not feasible.A third possibility, which leads to collider testable low-scale seesaw models, consists in heavy neutrino pairs forming a pseudo-Dirac structure.This leads to a cancellation between the two terms in (1.1) and can be justified by an approximate symmetry.The symmetry can be realised as an extension of the SM lepton number and we refer to it as a lepton number-like symmetry (LNLS) [8][9][10].It allows for the heavy neutrinos to be light enough, and at the same time have sufficiently large Yukawa couplings, to be directly observable at e.g. the large hadron collider (LHC) [11][12][13][14][15].When the LNLS is exact, the light neutrinos are massless and the two heavy neutrinos are precisely mass-degenerate, forming a single Dirac neutrino.However, if the LNLS is broken by small parameters, light neutrino masses are generated and typically a mass splitting between the heavy neutrinos is introduced.The corresponding seesaw scenario is referred to as symmetry protected seesaw scenario (SPSS) [8][9][10].The phenomenological effects of the slightly broken LNLS are captured by the phenomenological symmetry protected seesaw scenario (pSPSS) introduced in [10].Similar to the case of other neutral particles, this leads to particle-antiparticle oscillations [16][17][18][19], which are in this context called heavy neutrino-antineutrino oscillations (NN Os) [20].The quantum field theoretical (QFT) framework of external wave packets, which has been developed in [21][22][23][24] and applied to the case of NN Os in [25], allows to describe this phenomenon including decoherence effects [26].
Although the small breaking of the LNLS generates lepton number violating (LNV) processes, one can argue on principle grounds that for prompt heavy neutrino decays it is not possible to observe LNV at the LHC [27].In contrast, NN Os generate observable amounts of LNV when the heavy neutrinos are sufficiently long-lived compared to their oscillation period.When the heavy neutrinos have large enough lifetimes, such that their time of decay can be reconstructed, it can be possible to resolve the oscillation pattern of the NN Os.This allows to deduce the mass splitting of the heavy neutrinos, which might not be possible from measuring the amount of lepton number conserving (LNC) and LNV processes alone.The possibility of resolving the NN Os at colliders has been estimated to be feasible during the high luminosity phase of the LHC (HL-LHC) in [20].
The goal of the present paper is to explore under which conditions it is possible to resolve NN Os at the HL-LHC at the reconstructed level, employing the process shown in figure 2.
To this end, the FeynRules implementation of the pSPSS [28] is used to simulate NN Os in MadGraph.Subsequently, Pythia is employed to simulate QCD effects such as hadronisation, and the fast detector simulation Delphes is used to simulate the CMS phase II detector.A cut based analysis of the generated events is performed using custom C ++ code.Finally, the prospects to resolve oscillation patterns for several benchmark model (BM) scenarios of long-lived nearly mass degenerate heavy neutrinos at the HL-LHC are derived, using a detailed statistical analysis implemented in Mathematica.This paper is structured as follows: First, the pSPSS and NN Os are briefly reviewed in section 2. Subsequently, in section 3 the simulation of signal events containing these oscillations and relevant background processes are discussed.Afterwards, the statistical analysis is introduced in detail in section 4. Finally, the results are presented in section 5 and the paper is concluded in section 6.In appendix A we comment on the induction of residual oscillations through cuts on the transverse impact parameter.

Symmetry protected seesaw scenario
When considering two heavy neutrinos under the assumption of an intact LNLS, the only allowed additions to the SM Lagrangian are where N 1 and N 2 are sterile neutrinos written here as left-chiral four-component spinor fields.
The Higgs and lepton doublets of the SM are denoted by H and ℓ, respectively, and y 1 is the Yukawa coupling vector with components y 1α .The ellipses capture further contributions that can be generated by additional sterile neutrinos, but are assumed to be subdominant here.In this case the two terms in the seesaw formula (1.1) cancel precisely.In addition, the following symmetry breaking terms can be introduced When the Yukawa coupling vector y 2 times the vacuum expectation value (VEV) as well as the Majorana mass parameters µ M and µ ′ M are small compared to m M , the light neutrino masses are guaranteed to be small as well.Due to the approximate LNLS, it is possible for the heavy neutrinos to have a mass well below the W boson mass, while at the same time their coupling can be large enough to obtain a sizable number of events at e.g. the HL-LHC, without violating constraints from light neutrino experiments such as searches for neutrinoless double β (0νββ) decay.

Phenomenological symmetry protected seesaw scenario
In the SPSS with exact LNLS (2.1), the neutrino mass matrix for the interaction eigenstates , which results for all points in a decay width of Γ = 13.8 μeV.However, they vary in their mass splitting and consequently have different oscillation periods τ osc , which leads to different LNV to LNC ratios R ll , defined in (2.11).
where m D = y 1 v with the SM Higgs VEV v ≈ 174 GeV.The mass term of the Lagrangian after EWSB is given by This mass matrix does not generate neutrino masses.Furthermore, the two heavy neutrinos are mass degenerate with a relative phase of −i in their Yukawa couplings and can therefore be described as a single Dirac fermion.However, in the presence of the small symmetry breaking terms (2.2) the complete mass matrix reads where µ D = y 2 v.This mass matrix not only generates small masses for the light neutrinos, but additionally introduces a mass splitting between the two heavy neutrinos, that is also suppressed by the same small LNLS breaking parameters.Thus, the pair of two Majorana fermions can no longer be described as a pure Dirac particle, but as a pseudo-Dirac particle.The symmetry breaking terms can arise from specific low scale seesaw models, such as the linear or the inverse seesaw, which yield only the terms proportional to µ D and µ M , respectively.However, the complete mass matrix (2.5) can be generated from more complicated seesaw models.Since the symmetry breaking terms are very small, their phenomenological effects can often be neglected in collider studies with the notable exception of NN Os, which can be phenomenologically significant as they are an interference effect.
At leading order (LO), the NN Os are fully described by the mass splitting ∆m between the heavy neutrinos together with an additional damping parameter λ, capturing the potential decoherence effects discussed in section 2.3.Therefore, instead of adding the terms (2.2) to the Lagrangian (2.1), in the pSPSS, the mass splitting is directly introduced as a model parameter.Consequently, the masses of the heavy neutrinos are given by where θ = m D /m M is the active-sterile mixing parameter.A detailed description of the pSPSS can be found in [10].

Benchmark models
In order to obtain the significance with which NN Os could be observed at colliders, three BMs that differ only in their mass splitting are used and given in table 1.The Majorana mass parameter and the active-sterile mixing parameters for all BM points are chosen to be m M = 14 GeV and |θ µ | 2 = 10 −7 , leading to a decay width of Γ = 13.8 μeV.The active-sterile Figure 1: Contour lines for the number of expected displaced vertex events N as well as bands for the LNV over LNC event ratio R ll , cf. definition (2.11), for the three BMs defined in table 1.The contour lines for N apply to the CMS detector at the HL-LHC with L = 3 ab −1 and cuts as defined in sections 3.1 and 3.2.They have a nose-like shape and depend only on the heavy neutrino mass m and the active-sterile mixing parameter θ µ , with θ e = θ τ = 0. N is therefore identical for all three BMs, which differ only in their mass splittings ∆m.The position of the BMs in the parameter plane is indicated by a purple cross.The different ∆m of the BMs result in different R ll bands.The contour lines where R ll ∈ [0.1, 0.9] are shown in yellow and orange for BM1 and BM3, respectively.The relative position of the cross to the bands shows that the BMs have an R ll close to one.The grey area corresponds to the exclusion bounds from displaced vertex searches [14,15] and the shaded grey area to the bounds from searches for prompt LNV processes [11,12].Since the prompt searches rely on LNV signals, they apply only to models with an R ll close to one.mixing parameter corresponds to a Yukawa coupling of y µ = 2.55 × 10 −5 , while the Yukawa couplings to the electron and τ -lepton have been set to zero.The number of expected events with muons as final state leptons in conjunction with the cuts presented in sections 3.1 and 3.2, are shown in figure 1.It becomes clear that the chosen BMs lie comfortably beyond the current bounds [14,15].
In the minimal linear seesaw model, only one pseudo-Dirac pair of heavy neutrinos is added to the SM.This results in the lightest neutrino being massless and the mass splitting of the heavy neutrinos ∆m being identical to the mass splitting of the light neutrinos, cf.[20].The mass splitting of BM3 is chosen to represent this possibility, where the light neutrino mass splitting is taken from a recent global fit assuming inverse light neutrino mass hierarchy [29].However, in top-down realisations of low scale seesaw models, e.g. the low scale linear seesaw models in [30,31], it is common to have multiple pseudo-Dirac pairs of heavy neutrinos.Since the light neutrino masses get contributions from all these heavy neutrinos, it is possible for the pseudo-Dirac pairs to have smaller mass splittings than in BM3 without the need of cancellations in the mass matrix.For simplicity it might be assumed that the collider phenomenology is dominated by only one of the pseudo-Dirac pairs.We introduce two additional BM points reflecting this possibility.

Oscillations
Oscillations of neutral particles, and thus NN Os, can be described in the QFT framework of external wave packets.Compared to simple methods relying on plane waves, wave packets allow for a self-consistent description of oscillations.This is due to the finite uncertainty in energy-momentum and spacetime, that allows to produce a coherent superposition of nondegenerate mass eigenstates, while simultaneously makes it possible to introduce the notion of a travelled time and distance.The QFT framework also allows to calculate the possible suppression of oscillations due to the loss of coherence between the propagating mass eigenstates.For phenomenological studies these effects can be captured by a damping parameter λ, which is thus included in the pSPSS.
At LO, the probability to obtain a LNC or LNV event is given by Therefore, the oscillation period, defined in the proper time frame, is 2 Since the probability density of the heavy neutrino to decay is the probability for a heavy neutrino to decay between the proper times τ min and τ max and forming an LNC or LNV event is given by By integrating the oscillations from the origin to infinity, one derives the total ratio between LNV and LNC events [10,32,33] that can be directly deduced from cut and count based analyses.The expected number of events in a collider experiment is given by where the collider luminosity L , the heavy neutrino production cross section σ, and the branching ratio (BR) are used and D(ϑ, γ) is the probability density of the heavy neutrino to have a Lorentz factor γ and an angle ϑ with respect to the beam axis.The detector geometry is incorporated in the parameters τ min and τ max when boosting to the laboratory frame via The NN Os in the pSPSS have been discussed in detail in [10].We have checked that for the parameter points under consideration decoherence effects can be neglected and we thus take λ = 0 in the following [26].
Diagram depicting the production, oscillation, and decay of a heavy neutrino.The heavy neutrino interaction eigenstate N is produced in association with a prompt antimuon.Subsequently, the mass eigenstates oscillate such that finally a neutrino or antineutrino interaction eigenstate decays into a displaced muon or antimuon, respectively.We indicate that the process is initiated by proton collisions and that, for our parameter points, the two final quarks, originating from the hadronic W decay, immediately hadronise into a single jet.

Simulation
In order to simulate events exhibiting NN Os, we employ the FeynRules [34] implementation of the pSPSS [28].FeynRules is used to generate an UFO output, passed to Mad-Graph5_aMC@NLO 2.9.10 (LTS) [35] in order to generate events at parton level.We use the patch introduced in [10], that modifies the function calculating the particle's time of flight (TOF), to implement the oscillations.When evaluating the matrix element, MadGraph treats the process as prompt, adding the TOF information afterwards.Therefore, almost no LNV is obtained when interference between diagrams with different mass eigenstates are taken into account.In order to circumvent this, the process is initially specified in such a way as to prevent interference between the heavy neutrino mass eigenstates.This is achieved by writing the heavy neutrinos explicitly as intermediate states.As a consequence the total cross sections for the LNC and LNV process are identical and the amount of generated LNC and LNV events is, except for statistical variations, the same.The correct interference effects, including NN Os, are then included via the patch described in [10].The syntax to generate a process, as the one given in figure 2, is generate p p > ww, (ww > mu nn, (nn > mu j j)) where, in addition to the jet j and proton p multi particles, containing quarks and gluons, the multi particles define mu = mu+ mudefine ww = w+ wdefine nn = N4 N5 are used.Here, the initial W boson is additionally taken to be on-shell, since we focus on heavy neutrinos with masses far below the W boson mass.An additional hard initial jet is included by using add process p p > ww j, (ww > nn mu, (nn > mu j j)) where MLM matching was enabled with the standard parameter choices of MadGraph to prevent double counting.
MadGraph utilises Pythia 8.306 in order to hadronise and shower the parton level events [36].Finally Delphes 3.5.0 is used with the standard CMS phase II card CMS_PhaseII_0PU to simulate the detector effects [37].

Secondary vertex reconstruction and smearing
Since Delphes neither simulates displaced tracks properly, nor reconstructs secondary vertices, we implemented a smearing function affecting the displaced vertices, with the idea to capture experimental uncertainties.To our knowledge, no results for the overall precision of the vertex reconstruction have been published by CMS so far.Therefore, we vary it in our simulations between zero and 4 mm, which we assume to be a conservative parametrisation of our ignorance, relying on private communication with members on the experimental collaborations.We like to emphasise that it would be highly welcome if the experimental collaborations could provide such information and, ideally, if such uncertainties could be implemented in Delphes directly.In detail, it is assumed that each reconstructed vertex is distributed with a Gaussian, with a standard deviation of n mm, around its actual value.The true value of the displaced vertex is obtained from the displaced muon.The results are presented for different values of n, demonstrating the effects of the uncertainty in the vertex reconstruction.

Signal
In order to observe LNV via NN Os using the process given by the diagram presented in figure 2, the two leptons have to be measured.For the BMs given in table 1 and indicated in figure 1, the final state quarks are soft and immediately hadronise into a multi pion state, which then forms a single jet that can be captured with a cone radius of ∆R = 0.4.Furthermore, the BM points are chosen such that the heavy neutrinos are long-lived.Hence, the signal contains one prompt muon, one displaced muon, and one displaced jet.
A cut is used, which requires the muon to have an impact parameter of |d 0 | ≤ 100 μm, to ensure exactly one prompt µ.Furthermore, this muon is required to have a transverse momentum p T greater than 20 GeV, ensuring that the event can be triggered.Events that do not contain at least one displaced µ are excluded by a cut that requires at least one muon with an impact parameter of |d 0 | ≥ 2 mm.Furthermore, the displaced muon candidate is only valid if it emerges from the same vertex as the displaced jet.A jet is considered displaced if it contains at least two displaced tracks that originate from the same vertex within a radius of 100 μm.A track is considered displaced if it has an impact parameter of |d 0 | ≥ 2 mm.Signal events are required to have exactly one displaced j.
In addition to those basic cuts that define the signal, we demand µ isolation by requiring that the displaced muon and displaced jet have a ∆R larger than 0.4 in order to reject muons radiated from the jet.Furthermore, restricting the reconstructed heavy neutrino mass to the N mass window of ± 2 GeV around its theoretical value results in a cleaner signal.Generally, displaced particles are only considered if their origin falls inside a cylinder with dimensions given by half the tracker size in each direction.For the CMS phase II detector, this means that displaced tracks must have an origin with a transverse distance smaller than 60 cm, and a longitudinal distance smaller than 150 cm, from the primary vertex.This requirement ensures that the produced particles propagate through a volume of the tracker that should be large enough to facilitate detection.At the moment this is an optimistic assumption about the performance of the displaced track reconstruction in the CMS tracker, however, we think that pushing this capability is a worthwhile goal for the collaboration (see also [13] and references therein).Although the reconstruction of NN Os discussed here mostly relies on displaced vertices appearing in the inner tracker, it is also possible to reconstruct muon tracks and vertices using the muon chambers [38][39][40].2: Cut flow for the signal and the heavy hadron background events.The physical events are given at an integrated luminosity of L = 3 ab −1 with a signal cross section of σ = 951 ab.

Background
The cuts introduced above define the basic search strategy for the signal, which is based on reconstructing a displaced vertex.The main sources of background for such a process are given by heavy flavour SM processes generating long-lived heavy hadrons, interaction of particles with dense detector material, and cosmic ray muons [14,15].The heavy SM hadrons can travel macroscopic distances before they decay, potentially forming displaced vertices.However, this is only possible if those heavy flavour particles, and consequently their decay products, are highly boosted.Since the selection rules that define the signal already require the displaced muon and the displaced jet to have a |d 0 | > 2 mm, this background is strongly reduced [14,41].A W mass window cut around the reconstructed W boson mass of ± 20 GeV is employed to reduce the background even further.We simulated this background using the same programs as for the signal events.The MadGraph syntax we used is generate p p > mu all bb where the additional multi particle define bb = b bĩ s defined.An additional hard initial jet was included by using add process p p > mu all bb j where MLM matching was enabled with the standard parameter choices of MadGraph.The whole process was simulated using the model sm-no_b_mass, which employs the five flavour scheme in the definitions of the protons and jets.After the parton level events are passed to Pythia and Delphes, the above mentioned cuts and event selection rules were used, and it has been found that the background is eliminated entirely.The cut flow for the simulated signal and background is given in table 2.
Interactions of SM particles with detector material can also result in a displaced vertex signature and thus have to be considered as part of the background.A map of regions containing dense detector material is required in order to accept only displaced vertices that are outside that region.Simulating this background requires detailed knowledge about the detector structure and is beyond the scope of this work.However, requiring a displaced vertex to be reconstructed from enough good tracks reduces this background.In an experimental analysis, a sophisticated track reconstruction algorithm judges which tracks are good, in the sense that the track is likely to be produced from a charged particle rather than detector noise.With no specific insight about this algorithm, we define the tracks as being good that pass our cut and selection rules, which ensure that the charged particles traverse at least half the tracker.With the cuts introduced to define the signal, it is ensured that each event contains at least three displaced tracks from which the displaced vertex can be reconstructed.
In order to veto against cosmic ray muons, the vertex direction cut is implemented that requires the reconstructed heavy neutrino momentum to be in the same direction as the displaced vertex.For this the distance in (η, ϕ)-space between the reconstructed momentum p N and the displaced vertex direction d is limited to ∆R ≤ 1.5.We assume that this background is eliminated when using this cut in combination with additional timing information, therefore we have not attempted to simulate it.
With the cuts described above, the signal region can be assumed to be background free.Additionally we found that a prompt e veto, where prompt is again defined as |d 0 | < 100 μm, with a p T ≥ 20 GeV does not affect the signal, such that it could be used to eliminate further background if that becomes necessary.

Statistical Analysis
The number of events entering the cut based analysis is given by the luminosity times the cross section for the considered process.The number of expected events after the cut based analysis N exp is then given by where the efficiency factor captures how the cuts summarised in table 2 impact the number of signal events.The set of events surviving all cuts is labeled D all , and can be divided into LNC events D LNC and LNV events D LNV .These datasets are then binned in the proper time τ space and the number of events in the i-th bin are given by where | • | is the cardinality such that e.g.N exp = |D all |.The heavy neutrino TOF τ E is defined as the proper time at which it decays in the event E.

Hypotheses
The shape of the histograms describing the heavy neutrino TOF may be predicted by two hypotheses.In contrast to the null hypothesis, the alternative hypothesis features oscillations.

Null hypothesis
The hypothesis without oscillations M 0 is based on the assumption that, for each bin, the probability of the heavy neutrino superposition to decay in a LNC process is equal to the probability to decay in a LNV process.Therefore, the mean number of predicted events in the i-th bin is given by ) with 15 k events, demonstrating the proper time distribution after applying the cuts summarised in table 2. The inverse Gaussian distribution µ all i (M 0 ) described in (4.4) is depicted ( ) together with the prediction of the null hypothesis µ The inverse Gaussian is shown with the parameters given in step 1 on page 14, where the normalisation is approximated as N 0 = 15500.The distribution of the alternative hypothesis (4.5) is shown for LNC ( ) and LNV ( ) events of BM1, with example value α = 0.05 for the washout parameter.
The probability of a given bin count can then be computed assuming a Poisson distribution of bin counts around this mean.In principle, one would expect the predicted mean values of this hypothesis to follow an exponential due to the finite lifetime of the heavy neutrinos.However, the employed cuts change this distribution into a non trivial one, which can be approximated by the probability density function (PDF) of a generalized inverse Gaussian, described by four free parameters, yielding where N 0 denotes the overall normalisation, " is the index parameter, -is the mean of the distribution, andis a scale parameter.Additionally, τ i denotes the position of the middle point of the i-th bin, and K α (z) is the modified Bessel function of the second kind.The distribution is shown in figure 3, where the parameters correspond to the best fit point estimated in step 1 on page 14.

Alternative hypothesis
The second hypothesis M 1 includes oscillations, with an oscillation period given by (2.8), on top of the distribution described by the first hypothesis.Due to the imperfect reconstruction of the Lorentz factor, an additional washout effect obscures the oscillation pattern for larger τ .This effect is included by an exponential factor α, suppressing the oscillation.The prediction for the mean number of expected events in a bin is then given by where the two additional parameters ∆m and α are incorporated, resulting in a total of six free parameters.Again, a Poisson distribution of bin counts, around this mean, is assumed.One example of such oscillations is given in figure 3.

Likelihood ratio test
To test whether the hypothesis including oscillations is statistically preferred by the simulated data, we use a likelihood ratio test.The main idea is to decide if the null hypothesis, given by M 0 , can be discarded in favour of the alternative hypothesis, given by M 1 .The likelihood is denoted by P ({N i }, M) and describes the probability of finding the measured bin counts {N i } given the hypothesis M. It is given by the product of likelihoods for all bins in both, the LNC and LNV, cases.With the assumed Poisson distributed number of events around the mean value µ i in the i-th bin, the likelihood of a single bin is given by where, as introduced above, N i is the number of events in the i-th bin, which depends on the dataset and the binning.The best fit point for each hypothesis, given the bin counts {N i }, is evaluated by maximising the corresponding likelihood.A hypothesis at its best fit point is denoted with a hat, e.g.M 0 .Therefore, the likelihood ratio can be computed, given the bin counts, and yields Taking the washout parameter to infinity in the alternative hypothesis reproduces the null hypothesis.The two hypotheses are therefore nested and as a consequence the inequality3 holds and hence the likelihood ratio is restricted by A likelihood ratio close to zero means that the given binned data are much better fitted by the alternative hypothesis than by the null hypothesis.In contrast, a ratio close to one shows that there is no clear distinction between the two hypotheses for the given binned data.Since in practice the logarithm of the likelihood is better suited for numerical computations, we continue the discussion using the logarithm of the likelihood ratio (LLR) which is defined as The LLR ranges from zero to infinity, where now a value near zero states that both, the null hypothesis and the alternative hypothesis, produce an equally good fit of the binned data.Contrary, a high value states that the alternative hypothesis produces a better fit.
Probability Caused by statistical fluctuations of the bin counts around their predicted mean, the null hypothesis M 0 can also feature an oscillation pattern, as demonstrated in figure 4. Therefore, the value of a LLR alone does not contain enough information to decide whether oscillations in a given dataset are significant or not.To translate the LLR into a significance, it is crucial to know the probability that the null hypothesis produces the same LLR due to fluctuations.Per construction, the most likely LLRs produced by M 0 are small, and the larger the LLR the less likely it is to be produced.The distribution of LLRs, generated by statistical fluctuations, is the PDF of the LLR under the assumption that the null hypothesis is true.We label it P.The challenge is then to find a value k p , for which the probability of obtaining a Λ({N i }) ≥ k p through statistical fluctuations, is smaller than p.This means that p is the probability of rejecting the null hypothesis with respect to the alternative hypothesis, even though the null hypothesis is true.Given P, the value of k p can be obtained via While in the limit of large sample sizes the PDF is typically assumed to approach a χ 2 distribution, the form of the PDF for finite sample sizes is generally unknown.Therefore a simulation is performed, with the goal to sample the PDF, as follows 1.All simulated events that survive the cut based analysis are used to estimate the true distribution of the null hypothesis M 0 , which is depicted in figure 3.For those events, the best fit parameters of M 0 are given by -= 10.81, -= 17.28, and " = 0.71.The overall normalisation parameter N 0 is not relevant for the following steps as they only depend on the PDF, which is normalised to unity.
2. To obtain a set of events that follows the null hypothesis, taking statistical fluctuations into account, TOFs are taken according to the PDF from step 1.This is done until the set contains N exp signal events.
3. Using these TOF values, bin counts are computed using the binning parameters based on equation (4.2).
4. The LLR is computed based on the bin counts from step 3.
By repeating steps 2 to 4 the desired distribution can be obtained.The obtained distributions are well approximated by χ 2 distributions, where the degrees of freedom (DOFs) are treated as a free parameter.The values of the DOFs for the three BMs are given in table 3.
Significance With the so obtained PDF it is possible to compute the probability that a given LLR produced from statistical fluctuations of the null hypothesis M 0 .This probability p can BM χ 2 DOFs k( 5 be translated to a significance Z by comparing it to a standard normal distribution.A smaller probability translates to a larger significance and a higher threshold k p .A significance of 5 σ corresponds to a probability of p ≈ 2.87 × 10 −7 , which for BM1 corresponds to a LLR threshold of k p ≈ 30.94.See table 3 for the values of the other BMs.Therefore, in the case of BM1, a LLR greater than 30.94 can be interpreted as a discovery, where oscillations have been found with a significance larger than 5 σ.For very small LLRs, it is possible that the obtained probability is greater than 50 %.Since for the translation into a significance a standard Gaussian is used, cf.[42], the corresponding significance would be smaller than zero .However, in cases where the LLRs are so small, the result should just be interpreted in the way that no oscillations could be proven in the given data.The simulated events and the resulting statistical fit for two example oscillations are shown in figure 5.

Data pre-processing
Some pre-processing of the events is performed to stabilise the numerics.In principle, it is expected that for sufficiently large proper times alternative hypotheses with a washout parameter greater zero produce the same bin counts as the corresponding null hypotheses.Therefore, the LLRs are dominated by small proper times, for which the washout effect is also small.We have found that the fitting algorithm used by us gives more reliable results when we consider only these dominant oscillations.Therefore we do not use all events for the statistical analysis but consider only a window containing the first oscillations.The details of this restriction for each BM point are given in table 4. To help reduce noise in the data a Gaussian filter of radius r = 1, which corresponds to a standard deviation of σ = r/2 in bin space, is applied to the bin counts before the maximisation takes place.

Results
In this paper, the first ever analysis of NN Os at reconstructed level is performed.A detailed statistical analysis is employed to obtain a significance describing the feasibility to resolve the oscillations at the HL-LHC.To that end, three BM points with different oscillation periods, given in table 1, have been simulated.All BM points feature a mean mass of the heavy neutrinos of 14 GeV as well as an active-sterile mixing parameter of |θ µ | 2 = 10 −7 .This leads to a decay width of Γ = 13.8 μeV.The parameters are chosen such that the bounds of current collider searches are well evaded.While BM3 captures the mass splitting of the minimal linear seesaw model, producing the measured light neutrino data with an inverted hierarchy, BM1 and BM2 feature a smaller mass splitting as discussed in section 2.2.

Secondary vertex smearing Since
Delphes does not provide the experimental uncertainty to reconstruct secondary vertices, we have introduced a free parameter governing Gaussian smearing of the otherwise perfectly reconstructed secondary vertices.Figure 6a shows the The simulated events are based on BM3 with an oscillation period of 1.67 mm.The best fit parameters are given by an oscillation period of 1.63 +0.03 −0.04 mm and a washout parameter of α = 7.44 +17.76 −5.45 × 10 −2 .The LLR was found to be 5.29, which in this case yields a significance of 0.67 σ.
Figure 5: Examples for the best fit of the alternative hypothesis to the data.For each parameter point a luminosity of 3 ab −1 is used, resulting in a total of about 90 events after cuts.However, the number of events contributing to the fit is based on the range of τ values used as given in table 4 and therefore differs between the BMs.The fit has been performed based on the binning options in table 4. In these examples the secondary vertex smearing has been neglected.The bands around the oscillations depict the errors of one standard deviation, assuming a Poisson distribution for the event count in each bin.
expected confidence with which oscillations could be observed in the data for a given smearing of the secondary vertex.With our BM points, the number of signal events that survive the cuts, shown in table 2, is roughly 90.This assumes a luminosity of 3 ab −1 which is the expected total integrated luminosity of the HL-LHC [43].A factor of two in the number of events can easily be achieved by choosing BM points closer to the excluded region.Therefore we regard this analysis as conservative.For each data point in the plot, 100 LLRs have been computed to obtain a mean value and a standard deviation for the significance.The figure shows that for an oscillation period of 15 mm in proper time space, corresponding to BM1, a significance of 5.19 σ can be expected if no smearing is taken into account.Moreover, the significance is above 5 σ up to a smearing of 2 mm, dropping to 4.46 σ for a smearing of 4 mm.Parameter points with smaller oscillation periods are affected stronger by the smearing.This is expected since a smaller oscillation period in proper time space is related to a smaller oscillation length in lab space.A larger smearing therefore results in a stronger washout for smaller oscillation periods.
Lorentz factor reconstruction Since it is crucial to reconstruct the Lorentz factor on an event by event basis, we show in figure 6b how a better reconstruction of the Lorentz factor as well as a higher luminosity makes it possible to observe oscillations even for BM3.The reconstruction error directly effects the quality of the reconstructed oscillations in the proper time frame.Therefore, if the error is too large for the mass splitting one wants to resolve oscillations for, a higher event numbers only yields a limited improvement.This is shown by  4: Binning parameters for the different BM points, as well as the number of events used in the computation of the LLR.Note that for small samples of event surviving the cut based analysis, as it is the case here, the number of events participating in the computation of the LLR undergoes fluctuations, which contribute to the fluctuations in the obtained significances.
the lowest line in figure 6b.In contrast the higher lines, that represent a smaller reconstruction error of the Lorentz factor, benefit much more from a larger event number.
The quality of the reconstruction of the Lorentz factor is measured using the relative error of the Lorentz factor, which is defined as where γ gen is the true Lorentz factor of the heavy neutrino and γ reco is the reconstructed one.
The set of events used for the analysis yields an exponential distribution of the relative errors of the Lorentz factor.Most events have a small ∆γ, while only a small fraction of events have a large ∆γ.The overall quality in the reconstruction can be measured by the standard deviation of that exponential distribution.While a large standard deviation corresponds to many events with a large relative error, and therefore to a poor reconstruction, the opposite is true for a small standard deviation.Without any improvement, with respect to the Lorentz factor achieved in this analysis, the events of BM3 yield ∆γ's that result in an exponential distribution with a standard deviation of 0.16.By improving the reconstruction such that the standard deviation is reduced to 0.0096, the significance is improved from around zero to (3.37 ± 1.10) σ for 90 events.Additionally doubling the number of events yields a significance of (5.13 ± 2.28) σ.This scaling is justified since we assume that it is possible to improve the LHC analysis presented here, such that more signal events survive while the background is still eliminated.Additionally, it is possible to choose a BM point closer to the experimentally excluded region as mentioned earlier.Furthermore, future collider experiments, such as the ones at the FCC [44,45], might yield much higher luminosities as well as better reconstruction possibilities of the Lorentz factor of the heavy neutrino.

Mass splitting dependency
Oscillations can be used to resolve very small mass splittings.However, from the discussion above it is clear that if the mass splitting becomes too large, i.e. the oscillation period becomes too small, the reconstruction of oscillation patterns will be challenging.This is shown in figure 7 where the oscillation period has been varied, using a fast simulation described below.One can see that larger oscillation periods produce higher significances.It is expected that the significance drops again if the oscillation period reaches the mean lifetime of the heavy neutrinos, since then oscillations can not develop before the mass eigenstates decay.The fast simulation is based on the assumption that the kinematics of the events is independent of the oscillation period.With this assumption, the oscillation period has no impact on the cut based analysis.As a consequence, the sum of LNC and LNV events follows for all BM points the same distribution, given by the null hypothesis.Thus it is possible to obtain a large sample of valid signal events by combining the events passing all cuts of the three simulated BM points.It is then possible to give each event a new tag, describing if that events should be counted as  For this comparison no smearing has been taken into account.
LNC or as LNV.For this the TOF of the heavy neutrino is computed on a per event basis, using the relation where d is the position of the displaced vertex with respect to the primary vertex.The formula for the oscillation probability can then be used to tag the event based on the new oscillation period.At this point one has generated a sample of valid signal events with the new oscillation period.
After that, one can pick a random subset of this sample containing the physical number of events, that can be computed using (4.1).Subsequently, the analysis to obtain the significance can be applied to this subset of events.This strategy is orders of magnitude faster than performing the full Monte Carlo (MC) simulation and cut based analysis for each oscillation period separately.

Conclusion
In this paper, we have performed a first full analysis of NN Os at the reconstructed level.The simulations are based on the FeynRules implementation of the pSPSS introduced in [10].After the generation of events at parton level with a patched version of MadGraph, hadronisation and showers are simulated using Pythia.Subsequently, a fast detector simulation of the CMS detector has been preformed using Delphes.The uncertainty in the reconstruction of the displaced vertex has been implemented with a smearing function, that randomly selects a value of the displaced vertex around its true value, based on a Gaussian.The analysis of the events is performed using custom C ++ and Mathematica code.
In our analysis we have focused on three BMs within the pSPSS, defined in table 1, with heavy neutrino parameters conservatively chosen inside the region allowed by current experimental constraints.The BMs differ by the heavy neutrino mass splitting, which is largest for BM1 and smallest for BM3.Simulating events containing heavy SM hadrons, we have shown that with the event selection rules and cuts as defined in sections 3.1 and 3.2 the corresponding background is completely evaded.It has also been argued that with the given cuts other backgrounds that could not be simulated should also be evaded.Thus, the surviving signal events can be treated as background-free.The statistical method to obtain the significance with which oscillations can be found in the simulated data is described in section 4.
Our analysis shows that for small enough heavy neutrino mass splittings, corresponding to large enough oscillation periods, it is possible to discover NN Os with the CMS detector at the HL-LHC assuming 3 ab −1 integrated luminosity.The impact of the oscillation period, the displaced vertex smearing, the number of events, and the error in the reconstruction of the heavy neutrino Lorentz factor on the significance are shown in figures 6 and 7.For resolving the NN Os, it is important that the smearing is smaller than the oscillation length in lab space.Similarly, the variance of the TOF, due to the error in reconstructing the Lorentz factor, should be smaller than the oscillation period.This is the case for BM1, for which a significance of (5.01 ± 0.9) σ is obtained, assuming a smearing of 2 mm and around 90 total events relevant for the analysis, cf.table 4. For smaller oscillation periods, as in BM2 and BM3, the significance is below 3 σ even if smearing is not taken into account.
However, we like to stress that smaller mass splittings might also be resolved with higher significance if the reconstruction of the Lorentz factor is improved.This would not only increase the significance itself but also improve the effect for larger event numbers as shown in figure 6b.The event number could, e.g., be increased by choosing a parameter point closer to the experimentally excluded region, cf.figure 1.Additionally, it might be possible to increase the significance for smaller mass splittings by increasing the decay width of the heavy neutrinos, i.e. parameter points with increased Yukawa couplings or mass.Then the heavy neutrinos would decay faster and there would be more events in the first oscillation cycles such that resolving the pattern becomes more feasible.However, more events would be lost by the d 0 cut in this case.In order to study the interdependence between such considerations, a scan over a larger sample of benchmark parameters is necessary.In addition, the presented study might be improved by more sophisticated background reduction and by optimising the window of considered TOFs defined in table 4. In summary, we have shown that the HL-LHC offers the exciting possibility to not only discover the LNV induced by NN Os of long-lived heavy neutrinos, but also to resolve the NN O pattern.Reconstructing the oscillation period would allow to measure this mass splitting and therefore discover the pseudo-Dirac nature of the heavy neutrino pair.This would provide deep insight into the mechanism of neutrino mass generation and could help to shed light on whether leptogenesis is able to generate the baryon asymmetry of the universe (as discussed e.g. in [46]).discussed in this paper, the transverse impact parameter is given by [10] where d T is the transverse distance of the displaced vertex and the sine measures the angle in the transverse plane between the momenta of the heavy neutrino and the muon it decays into.This sine introduces an angular dependency, sensitive to spin correlations in the process under consideration.Since the LNC and LNV processes expose dissimilar spin correlations, the d 0 cut effects them differently.This leads to the observation of NN Os patterns in event samples that are a priory insensitive to the difference between LNC and LNV processes.As an example, figure 8 shows the residual oscillations appearing in a large sample of N LNC + N LNV events after introducing a d 0 cut.While the event sample with no cuts does not feature any oscillation pattern, it is shown that the d 0 cut results in residual oscillations, with peaks aligning with the ones of the LNV oscillation pattern.It can be concluded that the d 0 cut affect the LNC events more severely than the LNV ones.By contrast, a cut on d T is independent of spin correlations and thus no residual oscillations appear.We found that this effect is subdominant for smaller event samples, such as in this analysis, and therefore neglected it in the main part of the paper.

Figure 3 :
Figure 3: Example histogram () with 15 k events, demonstrating the proper time distribution after applying the cuts summarised in table 2. The inverse Gaussian distribution µ all i (M 0 ) described in (4.4) is depicted () together with the prediction of the null hypothesis µ

Λ = 5 .Figure 4 :
Figure 4: Three examples of statistical fluctuations in the null hypothesis, producing patterns that mimic NN O. From panel (a) to panel (c) the oscillatory pattern becomes more distinct while the probability that the depicted histogram is generated decreases.

cτ osc = 15
The simulated events are based on BM1 with an oscillation period of 15 mm.The best fit parameters are given by an oscillation period of 14.08 +0.85 −0.71 mm and a washout parameter of α = 3.66 +31.73 −3.66 × 10 −3 .The LLR was found to be 51.0, which in this case yields a significance of 6.66 σ.
Secondary vertex smearing for 90 signal events.Lorentz factor reconstruction error in BM3.

Figure 6 :
Figure 6: Panel (a): Significance of the three BM points at a luminosity of L = 3 ab −1 as a function of the secondary vertex smearing.Panel (b): Significance of BM3, as function of the number of events surviving the analysis, for three different relative errors of the reconstruction of the Lorentz factor (5.1).For this comparison no smearing has been taken into account.

Figure 7 :
Figure7: Significance as function of the oscillation period for 90 events.Three independently simulated BM points as well as eight points simulated using a fast simulation are depicted.Smearing of the secondary vertices is neglected.

Figure 8 :
Figure 8: Effects of the spin correlation sensitive cut d 0 , in comparison to the non-sensitive cut d T , on a set of N LNC + N LNV events.For each of the three datasets, represented by histograms, 5 × 10 5 generator level events, based on the values of BM3, have been simulated.All histograms are normalised with the same factor, that ensures that the area under the uncut histogram in (a) sums to unity.Panel (a): Comparison of the MC data before and after a d T cut of 4 mm.Panel (b): Impact of a d 0 cut of 4 mm overlayed by the pattern of LNV oscillations.Due to the angular dependence appearing in (A.1) the d 0 cut does not only generate a delayed onset but additionally imprints an oscillatory pattern originating in the LNV oscillations.

Table 1 :
All BM points have a mass of 14 GeV and an active-sterile mixing parameter satisfying 3)BM ∆m/ μeV cτ osc / mm

Table 3 :
Statistical fluctuations in the bin counts, following the null hypothesis, lead to a χ 2 distribution with given DOFs.The distributions are computed following the steps 1 to 4 on page 14.The LLR corresponding to a p value of 5 σ are found to be increasing from BM1 to BM3.