Decoherence effects on lepton number violation from heavy neutrino-antineutrino oscillations

We study decoherence effects and phase corrections in heavy neutrino-antineutrino oscillation (NN¯O\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ N\overline{N}O $$\end{document}s), based on quantum field theory with external wave packets. Decoherence damps the oscillation pattern, making it harder to resolve experimentally. Additionally, it enhances lepton number violation (LNV) for processes in symmetry-protected low-scale seesaw models by reducing the destructive interference between mass eigenstates. We discuss a novel time-independent shift in the phase and derive formulae for calculating decoherence effects and the phase shift in the relevant regimes, which are the no dispersion regime and transverse dispersion regime. We find that the phase shift can be neglected in the parameter region under consideration since it is small apart from parameter regions with large damping. In the oscillation formulae, decoherence can be included by an effective damping parameter. We discuss this parameter and present averaged results, which apply to simulations of NN¯O\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ N\overline{N}O $$\end{document}s in the dilepton-dijet channel at the HL-LHC. We show that including decoherence effects can dramatically change the theoretical prediction for the ratio of LNV over LNC events.


Introduction
The origin of the observed neutrino masses is one of the great open questions in current particle physics.When the new particles involved in the neutrino mass generation have masses close to the electroweak (EW) scale, it is possible to investigate this question at the Large Hadron Collider (LHC) and future accelerators.One possible extension of the Standard Model (SM) of elementary particles that explains the observed light neutrino masses is based on the introduction of sterile neutrinos, i.e. fermions which are uncharged under the gauge symmetry of the SM [ ], see also [ -].
When they form Yukawa interaction terms with the lepton and Higgs doublets and, in addition, have Majorana mass terms [ ], light neutrino masses can be generated, which are then of Majorana-type.However, when the sterile neutrino masses are around the EW scale, care has to be taken not to exceed the bounds on the light neutrino masses [ ].When the Yukawa couplings are not tiny, the smallness of the light neutrino masses is realised by an approximate lepton number-like symmetry.The sterile neutrinos then form pseudo-Dirac pairs of nearly mass-degenerate heavy neutrinos.Although lepton number violation (LNV) is significantly suppressed for prompt heavy neutrino decays, cf.[ ], it can lead to observable effects via the phenomenon of heavy neutrino-antineutrino oscillations (NN Os) [ , ], see also [ -].Since the light neutrino masses become zero in the limit of exact lepton number conservation (LNC), observing LNV processes is crucial for probing the origin of neutrino masses.
Due to the NN Os, the number of LNC and LNV events in a given process depends on the time difference between the production and decay of the heavy neutrinos.Recently it has been shown for a selected benchmark model (BM) point consistent with present constraints, featuring a long-lived pseudo-Dirac heavy neutrino pair, that NN Os could be resolved during the high luminosity phase of the LHC (HL-LHC) [ ], see also [ ].However, even when the oscillations are not resolvable, they can induce LNV.The total ratio of LNV over LNC events, R ll , can be used to quantify the effect.
Decoherence and phase correction effects on NN Os are so far unexplored at the quantitative level.Previous studies have used estimates to verify that decoherence effects can be neglected for the considered BM parameters, e.g.[ , ], or have assumed this to be the case.While decoherence can, in principle, depend on various parameters, it has to be a function of the mass splitting of the heavy neutrinos.This can be argued from the fact that for experimentally resolvable mass splittings, the pseudo-Dirac pair must reproduce the phenomenology of two separate Majorana neutrinos.In such cases, NN Os are expected to vanish.Thus, in regions where decoherence effects are relevant, the simple NN O formulae have to be modified.Phase corrections for NN Os have not yet been discussed.
One can calculate the possible decoherence and phase correction effects in NN Os using quantum field theory (QFT) with external wave packets.The formalism is discussed in [ ] and has been adapted to the case of NN Os in [ ].In [ ], the effective damping parameter λ is introduced, which contains the collective effects of decoherence onto the leading order (LO) oscillation formulae.In the present work, we explore how decoherence and a time-independent phase shift affect NN Os as well as the quantitative prospects for observing LNV.
The remainder of this publication is organised as follows: In section , we introduce the external wave packet formalism.Afterwards, in section , we describe the derivation of the damped oscillation probability for the general case and subsequently apply the results to the case of NN Os in the symmetry protected seesaw scenario (SPSS).We show that the effects of decoherence can be summarised by a damping parameter λ, leading to a simple extension of the oscillation formulae.Results for λ, including its impact on R ll and searches for LNV, are discussed in section .Finally, we conclude in section .Details of the analytical derivations of the oscillation probabilities are presented in the appendices.The detailed steps necessary to integrate the transition amplitude over the intermediate particles' momentum and travelled distance are presented in appendices A and B, respectively.The constant phase shift is discussed in appendix C. The algorithm to compute the damping parameter λ numerically is discussed in appendix D, where the kinematics of the considered process is simulated using the phenomenological symmetry protected seesaw scenario (pSPSS) introduced in [ , ].

External wave packet formalism
In this section, we derive the transition amplitude between two external states that are prepared as wave packets.This essential quantity is the main ingredient to derive an oscillation probability following the arguments made in [ , ].It is defined as a function of a distance (t, x) in spacetime where H(t , x ) is the interaction Hamiltonian and T is the time ordering operator.In comparison to the usual QFT approach, in which plane wave states | Φ(p) with momentum p are used, the initial | Φ(t , x ) and final Φ(t , x ) | states are wave packets centred at the indicated points in spacetime and can be written as a function of a plane wave state using where E(p) is the energy of the particle, and ψ(t, x, p, p 0 ) is the wave packet envelope which describes the shape of the wave packet and is centred around the momentum p 0 .Assuming that the external wave packets are Gaussian and approximating the matrix element at the peak of those Gaussian functions, it is possible to evaluate the momentum integration over the external wave packets, yielding the transition amplitude for a mass eigenstate i [ ] Here G i (s) is the denominator of the renormalised propagator with s = E 2 − |p| 2 , M i (E, p) denotes the interaction amplitude, defined as the matrix element without the denominator of the propagator, and N is a normalisation constant.The imaginary part of the exponent contains the phase Here x is the distance, and t is the time difference between the production and detection point.
The real part of the exponent contains the energy-momentum envelope (EME).This name is derived from the fact that it describes the shape of the intermediate particle's wave packet as a function of its energy and momentum.Its shape is defined by the shapes of the external particles' wave packets at the production P and detection D vertices V and is given by [ ] Quantities with a suppressed vectorial index are indicated by boldface.
The precise form of the normalisation constant N changes throughout the paper.However, the normalisation constant can always be evaluated using an appropriate normalisation condition, as discussed in section . .where Here E 0 and p 0 are the energy and momentum of the intermediate particle obtained from the peaks of the external particles' wave packets using energy-momentum conservation either at the detection or production vertex.They are thus called reconstructed energy and momentum.The (P → D) is a shorthand notation where quantities at production are replaced by similar quantities at detection.If the EME is approximated as a Gaussian, its width can be interpreted as the effective width σ eff of the intermediate particles' wave packet.The total energy and momentum widths are given by the reciprocal sum of the respective widths at the production and detection vertices 1 Each of these widths can be expressed in terms of the widths of the external particles in position space.The widths of the external particles in position space are parameters of the theory and are determined by the experimental situation under consideration.In the following, we only explicitly write definitions for quantities at production, while analogous definitions hold for quantities at detection.The energy and momentum widths at this vertex are given by where σ xPn is the width of the external particle n in position space and The velocity of the production region is defined by The particle with the smallest width dominates these terms unless its velocity is much smaller than the velocities of the other particles.Since it holds that [ ] one can calculate that the energy and momentum widths obey the inequality From the EME ( .), it follows that energies E and momenta p far from the reconstructed energy E 0 and momentum p 0 are exponentially suppressed, where far is defined according to the energy or momentum width, respectively.Additionally, damping from the propagator is expected when the reconstructed energy and momentum are such that the reconstructed mass is far from the intermediate particles' masses m i .This damping defines the shape of the resonance, which in plane wave QFT would be given by the Breit-Wigner distribution.
f (E, p) φ(t, x, E, p) Flowchart depicting the integration steps from the spacetime-dependent transition amplitude A i (t, x) of the mass eigenstate i to the time and distance-dependent oscillation probabilities P(t) and P(x), respectively.The terms appearing in the exponential part of the transition amplitude and the oscillation probability are the energy-momentum envelope (EME) f , the phase φ, the decay term γ, the dispersion term F , and the localisation term Λ.The last two terms give the main contribution to decoherence.Thus they dominate the damping parameter λ.

Derivation of the damped oscillation probability
The transition amplitude ( . ) needs to be integrated over the energy and momentum of the intermediate particle.The strategy to perform these integration steps is depicted in figure .While the energy integral is performed using the Jacob-Sachs theorem in the following section, the subsequent momentum integral is evaluated in appendix A. The final step is a distance average performed in appendix B, contrasting the time average used in [ ].In this publication, we express the oscillation probability as a function of elapsed time instead of distance since the relevant observables naturally depend on the proper time of the oscillating particles rather than their travelled distance.
For example, since the heavy neutrinos, once they are detected at a collider experiment, will exhibit a range of Lorentz boosts, the oscillation pattern has to be translated into the proper time frame of the heavy neutrino in order to be reconstructable, see [ , , ].Therefore, an oscillation probability as a function of time, averaged over the distance, is more suited for this purpose.The same applies to measurements of the R ll ratio, which is sensitive to the interplay between NN Os and the decay of heavy neutrinos and, hence, naturally defined in the proper time frame of the neutrinos.Furthermore, a distance average is more straightforward from a technical point of view since there are fewer distance-dependent terms than time-dependent terms in the relevant exponential, as can be seen in figure . .

Energy integration via Jacob-Sachs theorem
To further evaluate the transition amplitude ( .), the energy integral is evaluated using the Jacob-Sachs (JS) theorem [ ].This theorem states that for times larger than a threshold time t JS , which is estimated in section ., and for functions Ψ (E, p) that are non-zero only for a finite range of s = E 2 − |p| 2 , the energy integral can be approximately evaluated using where the complex pole energy is defined in terms of the complex pole of the propagator as while m i and Γ i are the mass and decay width of the mass eigenstate i, respectively.After this approximate integration, the transition amplitude reads The pole energy can be rewritten as where the decay width expansion parameter, i (p), and the mass eigenstate energy, E i (p), are defined as Under the assumption that the decay width is small compared to the mass eigenstate energy, the phase, ( .), can be expanded in the decay width expansion parameter, which yields The real part results in the phase of the mass eigenstate i, while the imaginary part generates an exponential decay term After the energy integration, the EME ( . ) of the mass eigenstate i can be approximated to be such that the LO term reads where The derivation from which follows that higher orders in the decay width expansion parameter can generically be neglected is presented in appendix A. . .However, for very short times, the O( i (p)) terms can lead to a time-independent phase shift, discussed in appendix C. In the numerical calculation presented in appendix D, these corrections are explicitly taken into account by identifying the imaginary part as a correction to the phase.Finally, the transition amplitude ( . ) after the energy integration ( . ) takes the form where next-to leading order (NLO) terms in the decay width expansion of the interaction amplitude M i (p) = M i (E i , p) are neglected.The remaining integrals are the three-momentum integral and an integral that averages over distance or time. .

Applicability of the formalism
The Jacob-Sachs (JS) theorem used in the previous section is only valid for times larger than the Jacob-Sachs threshold time t JS .It is defined via the diameter of the support of the intermediate particle's wave packet Since the interaction amplitude, and therefore the wave packet envelope, must vanish for values of √ s − m i larger than the uncertainties, this support is estimated by the experimental uncertainty in reconstructing the mass m i in [ , ].In the case of NN Os, where the mass of the heavy neutrino has to be reconstructed from semi-leptonic decay products, assuming this uncertainty to be of the order of one per cent of the heavy neutrino mass yields a time threshold of In order to have a fraction f of particles decaying beyond that time requires decay widths of where γ denotes the Lorentz boost factor.Demanding % of all particles decaying later than that time results in a JS width of when assuming a Lorentz boost factor of γ ≈ 10, which is a reasonable estimate for the parameter region considered in this work.
In contrast, the width of the wave packets of the external particles is estimated in section using the size of the silicon atom radius and the proton-proton distance in a beam bunch for final and initial states, respectively.This line of argument suggests that the neutrino wave packet should be zero outside a range defined by the width of the wave packet in the squared four-momentum Using the approximations derived in section ., the numerical values given in table and further approximating √ s = m, this estimate leads to a time threshold of where n is the number of standard deviations which is taken to define the support of the Gaussian distribution.Requiring a fraction f = 99 % of particles decaying later than that time leads to a decay width of Taking the σ range leads to a JS decay width of about .eV.In order to be conservative, we use this more restrictive value in the following. .

Dispersion regimes of the momentum integration
The integration over the three-momentum is carried out differently in three separate regimes depending on how fast the phase varies over the effective width.For slowly varying phases, the integral is evaluated using Laplace's method.This regime is called the no dispersion regime (NDR) since time-dependent dispersion effects can be neglected.In this regime, the argument of the exponential in the transition amplitude ( . ) is expanded up to second order in the momentum p around the position of the minimum of the EME at p i .Therefore, the momentum p i maximises the exponential of the EME The Hessian of the EME ( . ) with respect to the momentum is given at LO, i.e. neglecting the mass splitting, by (A. ) and defines the inverse of the effective width of the intermediate particle. 2ff The matrix structure of ˚0 is defined by the two vectors u P and u D .Therefore, there exists a vector which is orthogonal to both, and the corresponding eigenvalue is Due to the inequality ( .), this is the smallest eigenvalue leading to the largest effective width.The other two eigenvalues, which are dominated by the energy width, are therefore larger and approximately given by This approximation is justified when the velocity vectors are almost aligned u :≈ u P ≈ u D and the inequality ( . ) is large σ E σ p .The NDR applies to times shorter than the short-time threshold (A. ) Since, in its derivation, the phase is required to vary slowly over the effective width of the EME in all directions, the short-time threshold depends on the largest effective width and, therefore, the smallest eigenvalue of the Hessian.The detailed computation is described in appendix A. .When wave packets travel longer, the phase oscillates more rapidly as a function of p, such that Laplace's method, used in the short time regime, becomes unsuitable.Since the wave packets are broader in directions transversal to the reconstructed momentum p 0 , an intermediate regime exists in which Laplace's method can only be used for the longitudinal direction.In contrast, transversal directions are integrated using the method of stationary phase.This intermediate regime is called the transverse dispersion regime (TDR).The method of stationary phase Quantities with suppressed matrix indices are indicated by sans-serif font.yields p x = p y = 0, assuming that the longitudinal component, indicated by hatted variables, is p 0 = p 0 e z .Laplace's method in the longitudinal direction results in The Hessian at LO is given by (A. ) where the last approximation holds for u :≈ u D ≈ u P and when the inequality ( . ) is large, i.e.
σ E σ p .Similar to the definition in the NDR ( .), the effective width in the TDR is defined as 2 σ 2 eff = Σ −1 0 .
( . ) and the long-time threshold, which forms the upper bound for this regime (A. ), is defined by The computation leading to this result is presented in detail in appendix A. .Longer times t t long are not relevant for the discussion of heavy neutrinos in the parameter space of interest for this paper.However, for the respective regime, called the longitudinal dispersion regime (LDR), the distance-dependent formulae derived in [ , ] can be used.
The short-and long-time thresholds ( . ) and ( . ) are given in the lab frame.Using τ = t m 0 /E 0 they can be reexpressed in the proper time frame as For heavy neutrinos appearing in the process presented in figure , these regimes are depicted in figure after averaging over 100 events.The partition based on distances is presented in figure a.It shows that for experimental length scales smaller than about km, only the NDR and TDR are relevant, and the short-time threshold is of O(dm).The regimes in the proper time frame are shown in figure b.For decay widths leading to lifetimes comparable to the short-time threshold, it becomes relevant to quantify the fraction of events that fall into the NDR and the TDR, respectively.Decay widths Γ ≈ τ −1 short result in a fraction of 1 − e −1 events decaying before the threshold, and therefore inside the NDR.For decay widths Γ 10τ −1 short practically all events decay before the threshold.Contrary, for decay widths Γ 10 −1 τ −1 short , practically all events decay beyond the threshold, and therefore in the TDR.For decay widths Γ 10 peV, the LDR, becomes relevant.

. Time dependent oscillation probability
The probability for a superposition of mass eigenstates i and j to yield the transition between the given initial and final states, defined in the amplitude ( .), is given by The normalisation constant N can be evaluated using the condition outgoing where the sum is understood to include all possible processes, i.e. decay channels of the intermediate particle.Since mass eigenstates acquire a complex phase while propagating, the superposition of distinct eigenstates depends on a phase difference that varies with time and distance, leading to a periodic fluctuation of the probability.Therefore, we refer to the probability ( . ) as an oscillation probability.The position space integral in this probability is performed in appendix B. for the NDR and in appendix B. for the TDR.After this integration, the oscillation probability reads, according to results (B. ) and (B. ), where the product of the interaction amplitudes with their momenta evaluated at the peak of the intermediate particle's wave packet is From the definition of the oscillation probability ( .), it can be seen that it depends on the sum over the two mass eigenstates of the absolute value squared transition amplitudes.Since the EME ( . ) and the decay term ( . ) are real-valued, the probability depends on their sum where their values at the minimum of the EME are given by Feynman diagram depicting the heavy neutrino production at a hadron collider with subsequent oscillation and semi-leptonic decay.In the external wave packet formalism, each external particle's width is a free parameter and needs to be adjusted according to the experimental setting.We present our estimates for the external widths in table .Contrary, the exponential term describing the phase ( . ) is imaginary, such that the probability depends on the phase difference calculated in (B. ) and (B. ), where E 0 /m 0 is the Lorentz factor of the intermediate particle and τ (t) denotes the proper time the intermediate particle travels between production and decay.It vanishes when the mass splitting becomes zero, e.g. if i = j.This expression is modified by subdominant NLO terms [ , ] and augmented by a time-independent shift, see appendix C.
From the derivation of the localisation term Λ ij in (B. ) and (B. ) as well as the dispersion term F ij (t) in (B.a), it can be seen that they inherit this dependence on the mass splitting and that they are given by where the inverse of the effective width Σ 0 is defined in ( . ) and ( .).Both the localisation and the dispersion term are decoherence terms and thus lead to a damping of the oscillations.While the time-dependent dispersion term is absent in the NDR, it becomes relevant in the TDR. .

Heavy neutrino-antineutrino oscillation probability
From here on, we restrict to the LO phenomenology of symmetry-protected low-scale seesaw models and heavy neutrino-antineutrino oscillations (NN Os) appearing in processes such as the one presented in figure .The discussion is based on the symmetry protected seesaw scenario (SPSS), recently introduced with its minimal phenomenological version, the pSPSS, in [ ].The generation of light neutrino masses in seesaw models is directly related to the presence of LNV.
The process in figure is LNC if the two charged leptons have opposite charges and LNV if they have equal charge.In this scenario the oscillation probability ( . ) for these two possible processes takes the form In comparison to reference [ , equation .], the exponential has been replaced by the one of the oscillation probability ( . ) containing, apart from the phase, additional terms due to the wave packet nature of the involved particles.An additional term in [ ], which summarises the effects of the mass splitting in the interaction amplitudes, is neglected here since we treat the oscillations at LO.The factors of the leptonic mixing matrix at production α and decay β are collected in the terms For the SPSS this results at LO in [ ] for LNC and for LNV with i = j , for LNV with i = j .
( . ) where the active-sterile mixing angle is defined by with the SM Higgs vacuum expectation value v ≈ 174 GeV and the Yukawa coupling of one sterile neutrino labelled y, see [ ].The normalisation condition ( . ) for this scenario is evaluated for each flavour at production and yields 1 = ( . ) In the last step, it has been used that the sum of leptonic mixing matrix factors over LNC and LNV processes vanishes for i = j.Since the dispersion term, the localisation term, and the phase difference vanish for i = j, they are absent in the last line.
The oscillation probability ( . ) between the two mass eigenstates N 4 and N 5 is then given by where the damping parameter takes the form and can be expressed as Here sech(x) denotes the hyperbolic secant function, which is equal to one at the origin, and decays exponentially for values |x| 1.For the parameter region and time scales considered in this work, it is justified to assume that the two decay parameters are approximately equal such that From the EME ( .), it can be seen that its minimum goes to zero if m i = m 0 .However, heavy neutrinos with distinct masses cannot have m 4 = m 0 and m 5 = m 0 simultaneously.Therefore, we consider two more limiting cases: On the one hand, in cases where the reconstructed mass is near the mean of the heavy neutrino masses, with respect to the energy and momentum widths, the values of the EMEs are approximately equal f 4 ≈ f 5 .The normalisation then cancels these contributions, such that the damping factor becomes On the other hand, configurations in which either the EMEs or the decay terms are significantly different between the two mass eigenstates can lead to a damping of the oscillations.For example, for mass splittings much larger than the energy and momentum widths, the minima f 4 and f 5 are very different.The result is that one of the mass eigenstates is favoured by the available energy and momentum of the process, such that each event is dominated by one of the two Majorana particles, and the phenomenology is that of a pair of Majorana neutrinos without 1 but f 4 f 5 the damping parameter is given by This leads to significant damping if f 5 1.A similar argument holds for γ 4 (t) γ 5 (t).The interpretation, in this case, is that if one of the mass eigenstates decays much faster than the other, oscillations are significantly suppressed, and damping is large.
The reconstructed mass m 0 has to be near to one of the pole masses m 4 or m 5 since otherwise, the whole process is suppressed.This effect is similar to the resonant scattering, described by an s-channel process with an intermediate particle of mass m.For cases in which |s − m| Γ , the process is suppressed compared to |s − m| Γ .In the present case, s is labelled m 0 , and the width of the resonance is dominated by the energy and momentum widths.
While the definition of the damping parameter is derived in the context of NN Os in the SPSS, the presented strategies for its evaluation also apply to more general processes.

Damped heavy neutrino-antineutrino oscillations
For the simulation of the damped oscillations discussed in the previous section, the parameters that can impact the damping are • The masses of the heavy neutrinos.
• The decay widths of the heavy neutrinos, correlated with the time the heavy neutrinos propagate between production and decay.
• The momentum configuration of the external particles.
• The wave packet widths of the external particles.
The masses of the heavy neutrinos can be described in terms of their mean mass and their mass splitting Table : Position space widths assumed for the external wave packets of the incoming and outgoing particles appearing in the process presented in figure .For simplicity, we assume that the widths fall into three distinct classes.The widths of the incoming particles σ p are estimated by the LHC beam bunch proton-proton distance.The widths of the outgoing leptons σ l are estimated by the atom radius of the silicon in the detector, and the width of the outgoing quarks and W bosons σ j are estimated to be ten times as large.Panel (a) shows the different classes the widths fall into, and panel (b) shows our baseline estimates assumed for the simulations performed in this work.
The considered process and the heavy neutrinos' mean mass restrict the external particles' momentum configuration.Since the exact momentum configuration changes on an event-perevent basis, a general result is obtained by averaging the computed damping parameter over several events.Realistic momentum configurations are generated using the general purpose Monte Carlo (MC) generator MadGraph _aMC@NLO [ ] together with the FeynRules [ ] implementation of the pSPSS defined in [ , ].The numerical computation, obtained using the algorithm presented in appendix D, takes the decay time of the heavy neutrino into account and is accordingly performed either in the NDR or in the TDR.
In order to simulate the process shown in figure , the widths of the external particles' wave packets in position space need to be estimated.When the heavy neutrino is lighter than the W boson, the first W boson can be on-shell such that its width can be directly estimated.In contrast, the second W boson is off-shell, and the external widths of its decay products must be estimated.The situation is reversed if the heavy neutrino is heavier than the W boson.The wave packet widths in configuration space of the incoming particles σ p are assumed to be defined by the average distance between two protons in a beam bunch [ ].The wave packet width of the outgoing leptons σ l is assumed to be defined by the atom radius of silicon present in the detector material.Final quarks and the final W boson are expected to have a larger uncertainty than final leptons, such that their width σ j is given by 10σ l .See table for more details. .

Decay width dependence of the damping parameter
The mean decay width of the heavy neutrinos determines the time range the neutrino can propagate before it decays.It is thus possible to examine the time dependence of the damping parameter λ = λ 45 (t) by studying its dependence on the decay width.A numerical computation of the damping parameters for fixed mean masses of heavy neutrinos is presented in figure a.
The shape of the contours depicting constant damping consists of three regions: • To the right is a plateau stretching over several orders of magnitude.The plateau demonstrates that the effects due to varying decay widths, and with it, the time dependence of λ, are not significant in this region.The plateau can be understood from the result ( .), noting that neither the momentum differences p ij nor the matrix Λ ij contains any terms in Γ or t at LO.
• For smaller decay widths Γ 0.1 μeV the damping increases with decreasing decay width.This effect is due to non-identical group velocities of wave packets of different mass  Very small decay widths of Γ 10 peV are governed by the longitudinal dispersion regime (LDR), which is not calculated in this work.Therefore, the predictions for these events are simulated using the same techniques as for the events falling into the TDR and hence are not reliable.Note that in the m = 500 GeV case, the analytical TDR line and the numerical line coincide.
eigenstates, which causes the wave packets to separate over time and, in turn, causes decoherence.The effect becomes larger for heavy neutrinos that live longer.However, if only the first 100 oscillation cycles are considered, the effect vanishes, and the plateau in the central section continues for small decay widths.Alternatively, the effect also vanishes if decays inside a sphere of radius cm are considered.Since these two restrictions cover most phenomenologically interesting cases, the effects of decoherence due to the separation of wave packets can be neglected in the parameter region under consideration.In the following discussions, we assume these restrictions.They imply that the damping parameter depends, in addition to the width of the external wave packet in position space, only on the mean mass and the mass splitting of the heavy neutrinos, i.e. λ = λ(m, ∆m).
In figure the plateau extends beyond Γ = 0.1 eV, until where our numerical calculations are applicable, as we discussed in section . .Since the physics leading to the damping as a function of ∆m is independent of Γ at LO according to our analytical derivations in appendix C, we conjecture that we can extrapolate the plateau also to larger Γ .We make use of this conjecture when we analyse the consequence of damping on R ll in sections .and . .The analytical damping formula ( .), together with the approximated expressions for decoherence ( .), reproduces the plateau found in the numerical evaluation for both regimes, Note that the cm represents a somewhat randomly chosen value for which we have checked that the effects can be neglected.It does not represent a boundary at which those effects become relevant.
as shown in figure b.Since the time-dependent dispersion is disregarded in the NDR, the respective formulae do not feature the increased damping for small decay widths.In the region where they are applicable, the analytical formulae are in good agreement with the numerical results for m = 10 GeV and in almost perfect agreement for m = 500 GeV. .

Mass dependence of the damping parameter
With the restrictions established in the last section, the damping parameter is time-independent and can be studied as a function of the mean mass m and the mass splitting ∆m of the heavy neutrino.This dependency is presented in figure .The numerical results are shown in figure a, the results for the time-dependent analytic formulae ( . ) are presented in figure b, and the results for a distance-dependent oscillation probability, as obtained in [ , ], are shown in figure c.While the numerical results in the NDR and the TDR are approximately equal, the results for the approximated analytic formulae derived in appendices A and B differ between those regimes.Therefore, the results obtained from the analytical computation are presented for each regime individually, while the presented numerical results are valid in both regimes.
Since the time-dependent formulae for the damping factor are similar in the NDR and TDR, the LO effects of the momentum dependence can be explained by studying the time-independent part of the damping factor ( . ) where ( . ) and ( . ) are used, and the last approximation is obtained by observing that, in both regimes, the energy width at production dominates the reciprocal sum in the localisation term for the baseline estimate of external widths, defined in table b.Although the exact dependence of the damping factor on the mass is complicated since the process-dependent orientation of momenta and velocities change with varying mass, the sudden decrease of damping around the W boson mass can be explained by a change in the energy width.The energy width at production is given by a sum of all external particles at the production vertex.For heavy neutrinos lighter than the W mass, this includes the initial W boson and the initial charged lepton.For heavy neutrinos above the W mass, the initial W boson is off-shell, and thus, the relevant particles are given by the two incoming quarks and the initial charged lepton.An increase in the number of external particles participating in the production process results in a sudden increase in the energy width, which results in a sudden decrease in damping.
The shape of the contours describing constant damping is very similar for all computation methods in all regimes, except for the distance-dependent formulae in the NDR.This difference can be traced back to the localisation term in the NDR, which takes the form where ).The time-dependent formula is given in ( .), and the distance-dependent one can be found in [ , section . ., equation ].While the upper formula is proportional to the eigenvalues of ˚0, the lower formula is proportional to the inverse of the eigenvalues of ˚−1 0 .Therefore, as long as the reconstructed neutrino velocity v 0 is not parallel to one of the eigenvectors of the effective width, this lead to significantly different damping behaviour.
The absence of the dispersion term in the NDR can yield different damping for each regime.However, for the restrictions discussed in section ., the effects of this term are expected to be  negligible.Together with the fact that all other additional effects considered in the numerical derivation compared to the analytic one are small, the numerical and time-dependent analytical results in both regimes are expected to be approximately identical.This is confirmed by the results shown in figure .Therefore, the time-dependent results feature a smooth transition of the damping between the NDR and the TDR.In contrast, for the distance-dependent results, the  damping in the NDR differs significantly from the results in the TDR.This smooth transition can be seen as a further advantage of the time-dependent formulae, as derived in this work, over the distance-dependent ones. .

Wave packet widths dependencies
In order to understand the impact of the estimates for the widths of the external wave packets, we rederive the previous results after rescaling individual widths by a factor of one hundred and present their dependence on this scaling in figure .The effects on the relevant time thresholds are presented in figure a.They can be understood by considering their dependence on the external width given in ( . ) and ( .).While the momentum width σ p depends on the smallest width in configuration space at production and detection, which is given by the widths of the charged leptons σ l , the energy width σ E depends, according to the following argument, mainly on the proton width σ p .
The energy width depends on the two smallest widths at production and detection.While the smallest width is cancelled by the global factor of σ pP or σ pD , respectively, the second smallest width dominates the energy width ( .).For the baseline of external widths, see table , the energy width at production is much smaller than the corresponding energy width at detection.The reciprocal sum of those energy widths, precisely σ E , is thus dominated by the proton width.
In figure b, we show that the main impact on the damping parameter is due to the proton width by individually varying the external widths.This relation can be traced back to the fact that the damping is dominated by the energy width, see approximation ( .).
The effect of the non-monotonous behaviour of the contour around the W boson mass in figure b is reduced when varying the jet width.The non-monotonous behaviour can be explained by a change in the number of particles participating in production and detection, as described in section . .The change in the number of particles is the opposite for production and detection, such that it weakens if the energy widths at production and detection are similar.Multiplying the jet width by a factor of 100 has precisely the effect of equalising those energy widths.
From the above results, it becomes clear that the value of the external widths plays an essential role in the prediction of decoherence and merits further dedicated studies.

. Decoherence effects on R ll
The ratio between the number of LNV and LNC decays is called R ll .Since it is calculated by integrating over the NN Os, it is affected by decoherence.Therefore, it is necessary to know the amount of damping to predict its value as a function of, e.g., the mean mass and the active-sterile mixing parameter.
The probability of obtaining an LNC or LNV event, between proper times τ min and τ max , is given by the integral [ ] where τ = m 0/E 0 t is the proper time.Here, the decay probability density is given by and the oscillation probability is given by ( . ) where the time dependence of the damping parameter ( . ) has been made explicit by defining the parameter µ using the dispersion term in the TDR ( .).In the limit τ min → 0 and τ max → ∞ the integral and the ratio of LNV over LNC events is given by where the function that appears in both quantities is Here erfcx(x) is the scaled complementary error function, which decays exponentially for negative x approaches one for small x and is inversely proportional to x for large x.For a subleading time dependence in the damping parameter, the function can be approximated using The LO term corresponds to the limit µ → 0, which captures the NDR.In this limit the equations simplify to where the damping independent term f (0, 0) corresponds to the term that appears when taking furthermore the limit that also the constant localisation term can be neglected, λ → 0, which recovers for the ratio of LNV over LNC events the known result [ , ] ( . ) In cases where the damping is large λ(τ ) 1, the coherence between the propagating mass eigenstates is lost, and the phenomenology is that of two independent Majorana neutrinos.Increasing λ, therefore, increases the observed R ll compared to the naive case that does not take damping into account.This behaviour is depicted for time-independent damping in figure a.Therefore, when considering R ll (λ) as a function of ∆m and Γ , parameter regions with large damping must have a large R ll and lines representing a smaller R ll cannot penetrate those regions.As depicted in figure b, the contours representing the naive R ll are given by constant ratios between Γ and ∆m.However, damped R ll contours are bound from above by regions of large damping.Therefore, once damping becomes relevant, the R ll bands follow the contour lines of the damping parameter shown in figure . .

Decoherence effects on prompt searches for lepton number violation
For the part of the parameter space where heavy neutrinos decay promptly, direct discovery of oscillations by resolving them as proposed in [ ] is not possible.However, the integrated R ll ratio introduced in the last section can still be measured.Hence, it is relevant to predict this ratio for the realistic benchmark models (BMs) of the linear and inverse seesaw, introduced and motivated in [ ] and summarised in table .While the heavy neutrino mass splitting in the inverse seesaw scenarios depends on the active-sterile mixing parameter, defined in ( .), the mass splitting in the linear seesaw scenarios is fixed for each BM point.Therefore, in the inverse seesaw, it is always possible to restore coherence and the naive value of R ll by considering larger values of the active-sterile mixing parameter and, accordingly, smaller mass splittings.However, in the linear seesaw, the damping solely depends on the mass of the heavy neutrinos.When it becomes relevant, coherence is lost independently of the value of the active-sterile mixing parameter.The effects of decoherence onto the R ll bands of the linear and inverse seesaw BM scenarios are depicted in figure .For our baseline estimates of external particle widths, the damping becomes relevant for mass splittings ∆m 1 eV as shown in figure a.The presented BM points for the inverse seesaw feature such mass splittings for values of the active-sterile coupling in the range of 10 −4 < |θ| 2 < 10 −1 .Therefore, the R ll bands deviate from the naive ones in that region.As discussed above, the bands then follow the contour lines of the damping parameter, such that large damping results in an R ll of one.The mass splittings for the BM points for the linear seesaw are such that decoherence effects can be neglected for our baseline estimate of external particle widths, defined in table .Consequently, the R ll bands do not deviate from the naive ones in this case.
However, the situation changes if different widths of the external wave packets are assumed.As shown in figure b, the most significant effect on the damping is given by varying the proton width σ p .Hence the effect of this variation on the R ll bands are depicted in figure .If the proton width is divided by 100, the damping becomes relevant only above mass splittings ∆m 120 eV.This results in the deviation from the naive R ll only at smaller values of the active-sterile mixing parameter in the inverse seesaw, while the linear seesaw is still not affected as depicted in figure a.On the contrary, if the proton width is multiplied by 100, the damping already becomes relevant for mass splittings ∆m 10 meV.The R ll bands of the inverse seesaw BM points now deviate from the naive results already for larger values of the active-sterile mixing parameter, and the effects on the linear seesaw BM points are also relevant, as shown in figure b.The mass splitting of the linear seesaw with inverted ordering of the light neutrino masses is still too small for decoherence to become relevant.However, for the normal ordered The genuine results, taking damping into account as discussed in section ., are depicted as coloured bands, while the naive results that neglect damping are shown as grey bands.The mass splittings of the two BMs of the linear seesaw are too small to be affected by decoherence, as seen at the lowest two bands.However, the three BMs for the inverse seesaw, forming the uppermost three bands, deviate from the grey bands beneath them.Hence, for a fixed prediction of R ll , the value of the squared active-sterile mixing parameters between the genuine and the naive results can vary by up to four orders of magnitude.The grey area indicates displaced searches by CMS, and ATLAS [ , ] linear seesaw, damping must be considered.From figure b, it can be seen that the mass splitting of the normal ordered linear seesaw results in a small damping for masses m 10 GeV and 120 GeV m.However, for masses in the range 10 GeV m 120 GeV, corresponding to the local minimum observable in figure a, decoherence becomes relevant.The largest damping is around m ≈ 50 GeV, and thus the observed R ll deviates from the naive one in this region.Since, as discussed above, coherence cannot be restored by varying the active-sterile mixing parameter as in the case of the inverse seesaw, the observed R ll is close to one in this mass range for all values of the active-sterile mixing angle.Therefore, the R ll bands form a pole around m ≈ 50 GeV as shown in figure b.
The effects of the damping parameter are crucial when reinterpreting prompt searches for LNV signals.Regions that result in an R ll close to zero for a considered model do not yield any LNV events and, thus, are not restricted by these searches.However, if damping is significant, the true value of R ll might significantly deviate from the naive one, as discussed above.Therefore, decoherence effects can result in prompt searches for LNV becoming applicable in regions of the parameter space that seem unconstrained when neglecting decoherence.With enhanced damping, decoherence becomes relevant also for the linear seesaw with normal ordered light neutrino masses.Since the mass splitting of the heavy neutrinos in the linear seesaw is fixed, the damping parameter only varies as a function of the mass.Therefore, damping is relevant for masses in the range 10 GeV m 120 GeV, where the contours in figure a show a local minimum.

Conclusion
Low-scale symmetry-protected seesaw models generically predict the appearance of pseudo-Dirac pairs of heavy neutrinos.It is usually expected that for such models, LNV is severely suppressed [ ].However, these considerations omit the crucial possibility of heavy neutrinoantineutrino oscillations (NN Os) which can lead to a sizable number of LNV events if the oscillation period is shorter or of the same order as the lifetime of the heavy neutrinos [ ].In particular, for long-lived heavy neutrinos, the oscillation pattern between LNC and LNV events may be reconstructible in collider experiments [ ]. Apart from NN Os, decoherence can yield observable amounts of LNV by reducing the destructive interference between propagating mass eigenstates and, with it, the suppression of LNV.While this obstructs the reconstruction of the oscillation pattern, it also enhances the chances of observing LNV when the heavy neutrino's lifetime is smaller than its oscillation period.In this work, we have quantitatively studied this effect for the first time.
To this end, we have derived oscillation probabilities for NN Os in the framework of QFT with external wave packets, extending previous results, see [ ], with a reformulation that depends on the time difference between production and decay of the heavy neutrinos.The derivations presented here are not only technically simpler than results that depend on the distance between production and decay of the heavy neutrinos but also more readily applicable to cases where an interplay between the oscillation period and lifetime is relevant, such as the LNV over LNC event ratio R ll and the reconstruction of the NN O pattern, which requires a translation of the oscillations into the proper time frame of the heavy neutrinos [ , , ].The analytical calculation relies on expansions in small parameters and differentiates between the NDR and the TDR that apply for short and longer-lived particles.The numerical comparison shows that the time-dependent results for the NDR and TDR are almost identical, such that a smooth transition connects the two regimes.On the contrary, in the distance-dependent results, the differences between the two regimes are significant such that no smooth transition is possible.We have compared the time-dependent analytical results with a numeric calculation using MC data and confirmed a broad range of applicability.
In [ ], we proposed using a single damping parameter λ that encodes all decoherence effects in NN Os.Here, we provide the formulae necessary to calculate this damping parameter from first principles as a function of the heavy neutrinos' mean mass and mass splitting.To that end, we provide the conditions under which the dependence on the decay width Γ and, therefore, the dependence on the propagation time between the production and detection of the heavy neutrinos can be neglected.Based on our analytical results, we conjecture that the dependence of λ on m and ∆m also extends to Γ larger than the ones we calculated numerically.
Furthermore, we discussed a novel time-independent contribution to the phase and derived analytical formulae for this phase shift in the no dispersion regime and transverse dispersion regime.However, in the parameter region under the shift can be since it is either small or damping is large, which results in a suppression of the phase shift effects.
The employed framework depends on the widths of the external particles' wave packets in position space.These input parameters need to be adjusted to the described experimental situation.We have picked well-motivated values in order to present our results.Additionally, we have discussed the dependence of relevant quantities on changes in these parameters.We have also discussed possible limitations of applicability of the formalism from the JS time threshold.In our computations, we followed a conservative approach restricted to Γ < Γ JS = 0.1 eV.
We illustrate the impact of decoherence by presenting bands of R ll in the (m, |θ| 2 )-parameter plane and show significant deviations from the predictions when decoherence is neglected.Hence, we quantify for the first time how the transition of a coherently oscillating pseudo-Dirac pair to two independently acting Majorana particles affects the phenomenology and, therefore, the discovery prospects of symmetry-protected low-scale seesaw models.
From the results of this work, it is clear that the possibility of decoherence has to be considered when studying LNV signatures.Large decoherence not only suppresses the oscillation pattern but also disables the mechanism that suppresses LNV for a pseudo-Dirac pair.Therefore, the phenomenology of a pseudo-Dirac pair changes significantly in regions of parameter space that exhibit sizable decoherence.the phase varies slowly such that the integral is approximated around the maximum of the intermediate wave packet up to the second order in the momentum.This method is referred to as Laplace's method.If the phase varies rapidly over the intermediate wave packet, the method of stationary phase is used, where the largest contribution to the integral is obtained near the point for which the phase has an extremum.In the transverse dispersion regime (TDR), the phase varies slowly in directions transversal to the reconstructed momentum p 0 , while in the longitudinal direction, Laplace's method can still be used.The third regime, called the longitudinal dispersion regime (LDR), in which the method of stationary phase has to be used for all directions of the momentum p, is not relevant to the phenomenology considered in this work.In the following, we treat the NDR and the TDR separately and derive an oscillation probability in each of them.

A. No dispersion regime
For short times the phase ( . ) varies slowly, as a function of the momentum, over the size of the intermediate particles' wave packet given by the EME ( .).Thus the momentum integration can be performed using Laplace's method, where the integral is approximated around the momentum p i for which the EME together with the decay term ( . ) is minimal.Hence a necessary condition that needs to be fulfilled in order to apply Laplace's method is While the EME yields a Gauss-like shape for the amplitude with a maximum near the reconstructed momentum, the decay term favours large momenta.Therefore we are interested in the minima of the EME.For the analytic computation, it is assumed that the impact of the decay term on the position of the maximum is negligible, which can be justified if the decay term varies slowly over the width of the intermediate particles' wave packet.The condition for this assumption to be valid is derived in (A. ).However, in the numerical computation in appendix D, the effects of the decay term are taken into account.With the decay term neglected, the position of the maximum is expanded around the maximum of the reconstructed momentum . ) where the mass splitting expansion parameter and the NLO term are The matrix ˚−1 0 appearing in the NLO term of the expansion is the inverse of the Hessian of the EME at LO in the momentum expansion, which is given in (A. ).Its eigenvalues can be interpreted as effective widths of the wave packet of the intermediate particle mass eigenstates in different directions, see section . .The corresponding expansion of the energy of the mass eigenstate at the minimum of the EME around the reconstructed energy yields Finally, the expansion of the velocity of the mass eigenstate at the minimum reads A. .

Expansion
In order to evaluate the momentum integral in the amplitude ( .), each exponential term is expanded up to second order in the momentum p around the minimum of the EME at p i resulting in a function that can be evaluated using Gaussian integration as demonstrated in appendix A. . .

Energy-momentum envelope
The expansion of the EME ( . ) around its minimum at p i results in where the linear term vanishes since the evaluation is performed at the minimum.The constant term and the Hessian of the EME at the minimum read . ) where the constant term to LO in the mass splitting expansion defines the mass width where . ) and the LO term in the mass splitting expansion of the Hessian is From the inequality ( . ) follows that the Hessian can be estimated to take values within . ) where |˚i| denotes that the considered inequality has to hold for all eigenvalues of ˚i.
Jacob-Sachs integration Solving the energy integral via the JS theorem as demonstrated in section .introduces the complex pole energy ( .).The dependence on the mass eigenstate energy ( . ) can then be used to estimate the effects of the decay width expansion parameter ( . ) on the EME ) where . ) and the ellipses denote terms that do not depend on the decay width expansion parameter, and for simplicity, we only consider the term at production, keeping in mind that the same arguments hold for the equivalent term at detection.Furthermore, assuming that the energy and mass splitting expansion parameters are of the same order, the constant term and the Hessian of the momentum expanded EME (A. ) are to LO The effects of the decay width expansion parameter on the Hessian are subleading and can therefore be neglected.While the effects on the EME are relevant, the term itself does only contribute to the damping of the amplitude ( . ) if it differs greatly between different mass eigenstates.This is due to the normalisation discussed in section . .The effects of the decay width expansion parameter are therefore neglected in the analytical derivation by using ( . ) while they are taken into account in the numerical calculations in appendix D.

Phase
The expansion of the phase ( . ) around the minimum of the EME results in where the constant term, the linear coefficient, the Hessian, and the velocity of the i-th mass eigenstate are Since the inverse of ˚i can be interpreted as an effective width of the intermediate wave packet ( . ) and ( .), it can be used to quantify the condition that the phase varies slowly over the width of the intermediate wave packet.Replacing |p − p i | 2 in the expansion series of the phase (A. ) with 2|˚i| −1 and requiring that the linear and quadratic terms are small results in These conditions can be approximated to read Since the averaging over the distance (or time) in a later step ensures that v i t ≈ x, the linear condition is automatically satisfied.The quadratic condition can be approximated as . ) where the eigenvalue of ˚0 containing σ 2 p is used since it imposes the most restrictive condition, see ( .).For the same reason, the velocity-dependent parts in the contribution from the phase can be neglected.This condition defines the short-time threshold up to which this integration method is valid.It is given by . ) and it is usually sufficient to work with the LO approximation.

Decay term
The decay term ( . ) is also expanded up to second order in its momenta around the minimum of the EME yielding . ) where the constant term, the linear coefficient, the Hessian and the parameter appearing in ( . ) evaluated at the minimum of the EME are given by Similar to the conditions appearing in the evaluation of the phase (A. ), two conditions can be derived, ensuring that the decay term varies slowly over the width of the wave packet . ) For times earlier than the short-time threshold (A. ), these conditions become . ) The linear condition holds as long as the decay width is of the same order or smaller as the momentum width, while the quadratic condition is satisfied for the assumptions used in the decay width expansion during the JS integration in section . .

A. . Integration
The momentum integral can then be evaluated using a standard technique for multidimensional Gaussian integrals over the coordinates x with a symmetric positive definite matrix A and a linear term b The vector b contains the linear order coefficients appearing in the momentum expansion of the phase (A. ) and the decay term (A. ) and reads For the analytical derivation, we neglect the correction from the decay term and proceed solely with the linear coefficient from the expansion of the phase around the minimum of the energymomentum integral.In the numerical calculation presented in appendix D the minimum p i is computed for the sum of the EME and the decay term, such that the linear contribution i (t) is absent all together.The matrix A collects the Hessian matrices resulting from the momentum expansion of the EME (A. ), the phase (A. ), and the decay term (A. ) around the minimum of the EME and reads . ) For times earlier than the short-time threshold (A. ), it is justified to approximate this sum with just the contribution from the EME, see the conditions (A. ) and (A. ).Hence the Hessian from the EME expansion (A. ) together with the linear term of the phase expansion (A. ) integrated over the momentum using the Gaussian integral (A. ) yield the spacetime envelope (STE) For a given time and velocity, this term leads to an exponential damping of the transition amplitude for values of x far from v i t.Here, far is defined by the eigenvalues of ˚i.We have explicitly discussed the effects of the decay width expansion parameter after (A. ) and found the STE to agree with [ ].
The oscillation amplitude in the NDR after momentum integration is therefore given by . ) where the terms in the exponential are the constant EME (A. ), the time dependent decay term (A. ), the spacetime dependent STE (A. ), and the spacetime dependent phase (A. ).
For times later than the short-time threshold (A. ), derived in the previous section, the oscillations are fast compared to the width of the wave packet such that Laplace's method is no longer the preferred method for evaluating the momentum integral.However, this argument depends on the direction of the heavy neutrino momentum.The second order term of the phase (A. ) yields a contribution |p| 2 t/E for momenta orthogonal to v i , while momenta in the direction of v i obtain an additional Lorentz contraction factor that leads to (1 − v 2 i )|p| 2 t/E.Since at LO v i = v 0 , this factor slows down the oscillations in the direction of the reconstructed momentum.Laplace's method is therefore still preferred in the direction along p 0 , while in directions orthogonal to it, the method of stationary phase is used.

A. . Stationary phase
The linear term of the expanded phase (A. ) averages oscillations to zero for momenta transversal to x. Therefore it can be assumed that p is parallel to x.Additionally, the EME requires that p is parallel to p 0 , which yields that x has to be parallel to p 0 .Assuming that p 0 is in the z direction the method of stationary phase can be used for x| x and x| y which yields . ) resulting in p x = p y = 0. Using those, the argument of the exponential of the amplitude ( . ) does not depend on p x and p y anymore such that only the integration over p = p| z is left.This integration is done using Laplace's method, as in appendix A. .The argument of the EME ( . ) can then be expressed as . ) where . ) Similar to the NDR, the phase, the EME, and the decay term are all expanded around the momentum p i for which the EME is minimal.The effects of the decay term onto the position of the maximum are neglected in the analytical derivation, and the exact conditions for this approximation to hold are given in (A. ).The position of the maximum up to linear order in the mass splitting expansion parameter (A. ) is . ) the mass splitting expansion of the mass eigenstate energy reads . ) and the mass splitting expansion of the mass eigenstate velocity is The effects of the decay width expansion parameter are neglected here for the same reasons as in appendix A. . .

A. . Expansion
Just as in the NDR, the terms in the exponent of the amplitude ( . ) need to be expanded up to second order in the momentum in order to perform a Gaussian integration in appendix A. . .

Energy-momentum envelope
The expansion of the EME (A. ) around the momentum of the minimum is given by where the linear term vanishes since the expansion is evaluated at the minimum, while the constant term and the Hessian are given by . ) and their LO contributions in the mass splitting expansion are given by Phase The expansion of the phase ( . ) around the momentum of the minimum of the EME yields where the constant term, the linear coefficient, the Hessian, and the velocity of the mass eigenstate are given by . ) Similar to the case of the NDR, a time threshold can be obtained by requiring that the phase varies slowly over the width of the wave packet.The wave packet width is approximated by Σ − 1 /2 i and used to reexpress the momentum deviations.The conditions resulting from the requirement that the linear and quadratic terms are small are The linear condition is ensured by the distance average performed in the next section, while the quadratic condition defines the long-time threshold.A time threshold independent of the mass eigenstates can be defined by approximating it at LO in the mass splitting expansion (A. ), leading to the validity range of the TDR

Decay term
The expansion of the decay term ( . ) is the same as in the NDR (A. ) where all vector quantities are replaced by their corresponding longitudinal component.
The constant term, the linear coefficient, the Hessian and the terms appearing in ( . ) are given by . ) The conditions for the decay term to vary slowly over the width of the neutrino wave packet read . ) assuming t t long i , these conditions reduce to . ) Similar to the situation in the NDR (A. ), the quadratic condition is automatically satisfied.The linear condition requires the decay width to be of the same order or smaller than the effective momentum width of the wave packet, which is typically of the order of the energy width σ E and thus much smaller compared to the width in the NDR.Additionally, the factor E i /m i can result in a violation of the condition for ultra-relativistic particles.This reflects the fact that for such particles, the TDR is valid up to arbitrarily large times, such that eventually, the decay term becomes important.In the numerical estimation of the damping, the effect of the decay term is taken into account.

A. . Integration
The integral in the transition amplitude ( . ) can be solved using the general result (A. ) where the vector b, now scalar, contains the first order terms from the expansion of the phase (A. ) and the decay term (A. ) and reads . ) In the following, the decay width expansion parameter correction is neglected for the analytical derivation but taken into account for the numerical calculation.The Hessians resulting from the momentum expansion of the EME (A. ), the phase (A. ), and the decay term (A. ) are collected in A which is now a scalar and reads . ) For times earlier than the long-time threshold (A. ), the LO contribution is given by the Hessian of the EME, see conditions (A. ) and (A. ).For the numerical results, the contribution from the Hessian of the phase and the decay term are taken into account.The integration (A. ) over the momentum yields the STE . ) and the amplitude ( . ) after integration is given by . ) where the terms in the exponent are the EME (A. ), the decay term (A. ), the STE (A. ), and the phase (A. ).

B Distance integration
In order to obtain an oscillation probability, the amplitudes for different mass eigenstates are coherently summed over An average over the oscillation distance is performed to obtain a formula dependent on the time.For the NDR as well as for the TDR the only distance-dependent components in the amplitudes (A. ) and (A. ) are the phases (A. ) and (A. ) and the STEs (A. ) and (A. ).Since the STE restricts the values of t and x to a range, similar to how the EME ( . ) restricts the values of p and E, the probability is expanded up to second order around x ij , the position of the minimum of the STE, before the resulting Gauss-like integral is evaluated.Typically, neither the oscillation distance x nor the oscillation time t are measured with perfect precision.Therefore, the oscillation probability (B. ) is either integrated over a span of time, as in [ ], or over a region of space in order to address the experimental uncertainty.In the following, the second possibility is employed, which yields an oscillation probability as a function of time.The distance-integrated oscillation probability reads where ∆x labels the experimental uncertainty in determining the distance travelled by the intermediate particle between its production and decay.This distance has to be larger than the width of the intermediate particle in spacetime, which describes the minimal uncertainty due to the wave packet nature of the intermediate particle.Due to the exponential damping stemming from the STEs (A. ) and (A. ) for values |x − x 0 | |˚0|, the distance integration can be taken to infinity ∆x → ∞.

B. No dispersion regime
Since the oscillation probability (B. ) depends on the absolute value square of the transition amplitude, the EME (A. ) and (A. ) need to be summed The same holds true for the STE (A. ) Since dispersion effects are neglected in the NDR, the STE can be approximated by its LO term in the mass splitting expansion ) Since it is quadratic in the spacial coordinates, there are no constant or linear terms.The position of the minimum of the STE is thus at LO given by Expansion The expansion of the STE in x around x ij (t) reads at LO in the mass splitting expansion Since dispersion is neglected, there is no constant term.The phase (A. ) is rewritten in terms of the deviation from the maximum in position space x − x ij (t) yielding Since the oscillation probability (B. ) is proportional to the absolute value square of the transition amplitude, the phase difference needs to be considered.The constant term encodes the phase difference in time, while the linear coefficient consists of a momentum difference Using the mass splitting expansion of the energy (A. ) and momentum (A. ) yields and the appearing difference in the mass splitting expansion parameter (A. ) can be approximated to be It is used that the reconstructed mass m 0 cannot be far from the mean of the mass eigenstate masses m, since otherwise, it is not possible to have p i ≈ p 0 and E i ≈ E 0 for both mass eigenstates at the same time, and thus one of the two f i terms in the amplitude (A. ) would lead to large damping of oscillations.Therefore, the phase difference is at LO in terms of the proper time Integration The integral (B. ) can be evaluated as a Gaussian integral, using relation (A. ) with the linear term from the phase expansion (B. ) and the Hessian from the STE expansion (B. ) resulting in the time-independent localisation term The transition amplitude after this distance integration is, therefore, where the STE and the decay term are given by (B. ), the phase is given by (B. ), and the localisation term is given by (B. ).

Transverse dispersion regime
The calculations for the TDR are very similar to the ones described in appendix B. .However, compared to the NDR dispersion is not neglected, and therefore, the approximation v i ≈ v j does not apply, which leads to a separation of wave packets of different mass eigenstates over time.
The sum of the EME (A. ) and the decay terms (A. ) are The distance average is performed by extending the terms in the exponent of the transition amplitude (A. ) around the minimum of the sum of the STEs (A. ) of the two mass eigenstates appearing in the oscillation probability The minimum is given at the position Expansion The expansion of the STE yields where the constant term and the Hessian are using the expansion (A. ) the velocity difference can be written as Since the expansion is performed around the minimum, the linear term is still absent.However, since dispersion is not neglected, the constant term does not vanish.The phase can be rewritten similarly to the NDR in terms of a constant and a linear term Since the phase contains the imaginary contributions to the amplitude, the difference between the two mass eigenstates appears in the oscillation probability the constant term and the linear coefficient are Using the expansions (A. ) and (A. ) yielding together with the approximation (B. ), leads for the phase difference to . ) which is identical to the phase difference in the NDR (B. ).

Integration of the oscillation probability
The integration P(t) ∝ Finally, the oscillation probability reads where the EME and the decay term are given by (B. ), the phase difference is given by (B. ), the dispersion term is given by (B.a), and the localisation term is given by (B. ).

C Phase shift
For short times, corrections to the LO expression of the phase ( .), which scale linearly with the decay width expansion parameter i become relevant.The expansion in the mass splitting expansion parameter δ i is given in (B. ) and (B. ) in detail.Therefore, we derive an analytical correction to the phase for parameter points in which i is non-negligible and test how the numerical phase as obtained via the algorithm described in appendix D deviates from the LO expression for the phase ( .).
As presented in section ., the NLO contribution in i to the phase results in the usual decay term γ i (t, p) ( .).However, there is also a contribution from the EME, which yields While the direct contribution of the correction term to the phase φ i (t, x, p) is negligible, it has a significant effect on the STE (B. ) since the corrections appear in the linear term This correction to the STE results in a shift of the position of its minimum (B.), which is now at The correction becomes relevant for large decay widths since then the particle's lifetime becomes small, while simultaneously, the correction term becomes larger.The resulting phase is then given by φ ij (t) = φ ij + m ij τ (t) , (C. )  where the constant phase shift is Since ˚0 is symmetric, this can be simplified to read Numerical results for the phase shift are shown in figure .As can be seen from figure a, in the considered parameter region, the total value of the modulus of the phase shift is small for most of the parameter space except for part of the region with large ∆m, where it can get large.However, we remark that in this parameter region, many oscillations take place before the heavy neutrino decays, and thus the phase shift per oscillation is still small.This is illustrated in figure b, which shows the modulus of the relative phase shift Furthermore, comparing with figure , one can see that in the region with a large total phase shift, the damping λ is very large, such that the oscillation term is strongly suppressed, making the phase shift practically unobservable.In summary, the numerical results show that, for current collider simulations, one can safely neglect the phase shift in the considered parameter region.However, outside the applicability region the phase shift can become significant, which can be seen from (C. ).The numerical results match this behaviour.Since these results are based on values of the decay width larger than Γ JS , and since we do not see a physical reasons for the phase shift to become arbitrarily large, we neglect the phase shift for the main part of this work.for i in {4, 5} do External wave packet formalism Derivation of the damped oscillation probability .Energy integration via Jacob-Sachs theorem . . . . . . . . . . . . . . . . . . . . . .Applicability of the formalism . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Dispersion regimes of the momentum integration . . . . . . . . . . . . . . . . . . .Time dependent oscillation probability . . . . . . . . . . . . . . . . . . . . . . . . .Heavy neutrino-antineutrino oscillation probability . . . . . . . . . . . . . . . . .Damped heavy neutrino-antineutrino oscillations .Decay width dependence of the damping parameter . . . . . . . . . . . . . . . . . .Mass dependence of the damping parameter . . . . . . . . . . . . . . . . . . . . . .Wave packet widths dependencies . . . . . . . . . . . . . . . . . . . . . . . . . . . .Decoherence effects on R ll . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Decoherence effects on prompt searches for lepton number violation . . . . . . . .

Figure :
Figure : Partition of the parameter space into longitudinal dispersion regime (LDR), transverse dispersion regime (TDR), and no dispersion regime (NDR) through the long-and short-time thresholds as a function of the heavy neutrino mass m.Panel (a) shows the partition in the lab frame distance ct as defined in ( . ) and ( . ) and panel (b) shows the partition in the proper time τ frame as defined in ( .).
JS(b) Validity of approximations.

Figure :
Figure : Damping parameter λ as a function of Γ and ∆m.The contour lines for the damping parameter are shown for two different masses m as a function of the decay width Γ and the mass splitting ∆m.The numerical results for the damping parameter λ averaged over 500 events per parameter point are given four fixed values of exp(−λ) and in panel (a).The comparison between the numerical results and the analytical approximations for the no dispersion regime (NDR) and the transverse dispersion regime (TDR) are given for λ = ln 2 in panel (b).The mass dependence is given in more detail in figure .Very small decay widths of Γ 10 peV are governed by the longitudinal dispersion regime (LDR), which is not calculated in this work.Therefore, the predictions for these events are simulated using the same techniques as for the events falling into the TDR and hence are not reliable.Note that in the m = 500 GeV case, the analytical TDR line and the numerical line coincide.

Figure :
Figure : Simulation results for the damping parameter λ for four values of exp(−λ) as a function of the mass m and the mass splitting ∆m using the numerical analysis from appendix D in panel (a), the time-dependent results from appendix B in panel (b), and the distance-dependent results from [ ] in panel (c).Panel (d) compares these techniques.For the analytic results in panels (b) to (d), the results for the no dispersion regime (NDR) and transverse dispersion regime (TDR) are shown separately.

Figure :
Figure : Dependence of the time thresholds and the damping parameter on the scaling of each external wave packet width.The values for the widths of the baseline scenarios are given in table .The baseline scenario for the time thresholds shown in panel (a) is given in figure , and the baseline scenario for the damping parameter in panel (b) is given in figure a.The short-and long-time thresholds in panel (a) are most sensitive to the scaling of the lepton σ l and the proton σ p widths, respectively.The damping parameter in panel (b) is most sensitive to the scaling of the proton width.
Figure :Decoherence effects on the R ll ratio for values in R ll ∈ [0.1, 0.9].Panel (a) depicts the functional dependence given in ( . ) as a function of the naive value ( .), i.e. assuming λ = 0 and the damping parameter.The simulation results for the impact of decoherence effects in the (Γ, ∆m) parameter space are presented in panel (b).Results for the damping parameter λ are used, as discussed in section . .
Normal ∆m ν = (41.46± 0.29) meV Inverted ∆m ν = (749 ± 21) μeV Inverse m ν |θ| −2 m ν = 0.5 meV m ν = 5 meV m ν = 50 meV Table : Summary of the five benchmark models (BMs) considered in this publication, resulting in various mass splittings of the heavy neutrinos ∆m.Since the SPSS with one pseudo-Dirac pair captures the minimal linear seesaw, it suffices to define one BM for each light neutrino mass hierarchy that reproduces the observed mass splitting between two massive light neutrinos ∆m ν [ ].In the case of the inverse seesaw, the single pseudo-Dirac SPSS represents an incomplete theory since the model generates only one of the two light neutrino masses.Therefore, we define the mass splitting as a function of the generated neutrino mass m ν , for which three BM values are defined.

Figure :
Figure :Bands of the LNV over LNC event ratio R ll ∈ [0.1, 0.9] for the five BMs defined in table .The genuine results, taking damping into account as discussed in section ., are depicted as coloured bands, while the naive results that neglect damping are shown as grey bands.The mass splittings of the two BMs of the linear seesaw are too small to be affected by decoherence, as seen at the lowest two bands.However, the three BMs for the inverse seesaw, forming the uppermost three bands, deviate from the grey bands beneath them.Hence, for a fixed prediction of R ll , the value of the squared active-sterile mixing parameters between the genuine and the naive results can vary by up to four orders of magnitude.The grey area indicates displaced searches by CMS, and ATLAS[ , ].The shaded grey areas indicate searches for LNV signals by CMS and ATLAS [ , ].The two black lines at low masses indicate the reach of the HL-LHC and the FCC-ee [ -].The two black lines at high masses indicate the reach of the Large Hadron Electron Collider and FCC-eh [ ].
Figure :Bands of the LNV over LNC event ratio R ll ∈ [0.1, 0.9] for the five BMs defined in table .The genuine results, taking damping into account as discussed in section ., are depicted as coloured bands, while the naive results that neglect damping are shown as grey bands.The mass splittings of the two BMs of the linear seesaw are too small to be affected by decoherence, as seen at the lowest two bands.However, the three BMs for the inverse seesaw, forming the uppermost three bands, deviate from the grey bands beneath them.Hence, for a fixed prediction of R ll , the value of the squared active-sterile mixing parameters between the genuine and the naive results can vary by up to four orders of magnitude.The grey area indicates displaced searches by CMS, and ATLAS[ , ].The shaded grey areas indicate searches for LNV signals by CMS and ATLAS [ , ].The two black lines at low masses indicate the reach of the HL-LHC and the FCC-ee [ -].The two black lines at high masses indicate the reach of the Large Hadron Electron Collider and FCC-eh [ ].

Figure :
Figure :Impact of scaling the proton σ on the ll ratio.The naive R ll bands are depicted as grey lines, and the coloured lines represent the R ll taking decoherence into account.Dividing the proton width by a factor of 100 reduces damping, as shown in panel (a).Therefore, the effects in the inverse seesaw are only significant for smaller values of the active-sterile mixing parameter.Contrary, multiplying the proton width by a factor of 100 results in enhanced damping, see panel (b).The effects in the inverse seesaw are then already significant for larger values of the active-sterile mixing parameter.With enhanced damping, decoherence becomes relevant also for the linear seesaw with normal ordered light neutrino masses.Since the mass splitting of the heavy neutrinos in the linear seesaw is fixed, the damping parameter only varies as a function of the mass.Therefore, damping is relevant for masses in the range 10 GeV m 120 GeV, where the contours in figure a show a local minimum.

Figure :
Figure : The absolute value of the phase shift φ ij in panel (a) and the absolute value of the relative phase shift φ rel ij , defined in (C. ), in panel (b) as functions of the decay width and mass splitting.