Heavy Neutrino-Antineutrino Oscillations in Quantum Field Theory

It has been proposed that the coherent propagation of long-lived heavy neutrino mass eigenstates can lead to an oscillating rate of lepton number conserving (LNC) and violating (LNV) events, as a function of the distance between the production and displaced decay vertices. We discuss this phenomenon, which we refer to as heavy neutrino-antineutrino oscillations, in the framework of quantum field theory (QFT), using the formalism of external wave packets. General formulae for the oscillation probabilities and the number of expected events are derived and the coherence and localisation conditions that have to be satisfied in order for neutrino-antineutrino oscillations to be observable are discussed. The formulae are then applied to a low scale seesaw scenario, which features two nearly mass degenerate heavy neutrinos that can be sufficiently long lived to produce a displaced vertex when their masses are below the $W$ boson mass. The leading and next-to-leading order oscillation formulae for this scenario are derived. For an example parameter point used in previous studies, the kinematics of the considered LNC/LNV processes are simulated, to check that the coherence and localisation conditions are satisfied. Our results show that the phenomenon of heavy neutrino-antineutrino oscillations can indeed occur in low scale seesaw scenarios and that the previously used leading order formulae, derived with a plane wave approach, provide a good approximation for the considered example parameter point.


Introduction
After neutrino oscillations have been proposed by Pontecorvo in the late 50's [1], they have lead to great insight into the (light) neutrino parameters. First introduced as neutrino-antineutrino oscillations, Maki, Nakagawa and Sakata considered oscillations into different flavours in 1962 [2]. Pontecorvo proposed in 1967 the possibility of solar neutrino oscillation, after which the solar neutrino problem was discovered [3]. In 1969, Pontecorvo and Gribov proposed neutrino flavour oscillations as a possible solution to this problem [4]. Since then, precise measurements of the light neutrino oscillations resulted in experimental values for the two mass squared differences of the light neutrinos as well as of the three mixing angles of the (PMNS) lepton mixing matrix and have provided a first indication of the range of the Dirac CP phase.
Despite this great success regarding the light neutrino parameters, the origin of the light neutrino masses, which requires an extension of the present Standard Model (SM) of elementary particles, is still unknown. One way to generate them consists in introducing right-chiral neutrinos as SM singlets. They can have Majorana mass terms as well as Yukawa couplings to the left-chiral SM neutrinos and the Higgs doublet. After electroweak symmetry breaking, the particle spectrum contains the three light neutrinos plus additional heavy neutrino mass eigenstates. If such heavy neutrinos are within reach of collider experiments, and if they are sufficiently long-lived, oscillations among the heavy neutrino and antineutrino interaction eigenstates 1 could lead to great insight into the heavy neutrino parameters and thus into the neutrino mass generation mechanism.
In parallel to the applications of neutrino oscillations also the framework in which they are described has evolved. The first approach to neutrino oscillations has been a quantum mechanical description in which the neutrinos are treated as plane waves. For neutral K meson oscillations, a quantum field model including wave packets was proposed in 1963 in Ref. [5]. In 1981, the authors of [6] pointed out conceptual problems of the plane wave approach for light neutrino oscillations and suggested a wave packet treatment as solution. To resolve the remaining difficulties, [7] introduced a quantum field theoretical model similar to [6], in which the propagating particle (assumed stable) is treated as an internal line in a Feynman diagram and the external particles are described by wave packets. A review of the existing quantum field theoretical approaches is given in [8], where also a framework employing external wave packets is discussed that can be used to describe the oscillations of unstable particles.
To describe heavy neutrino-antineutrino oscillations, a density matrix formalism (based on [10]) has been used in [11] for heavy neutrinos produced from meson decays. Using this formalism, formulae for the oscillation of the LNC and LNV decay rates have been calculated. The authors of [12] used a formalism for meson oscillations and plane wave arguments to derive formulae for heavy neutrino-antineutrino oscillations. These formulae are then applied to heavy neutrinos produced from W R decays in a left-right symmetric extension of the SM, for heavy neutrino parameters testable at the LHC. There it was shown that without resolving the heavy neutrino-antineutrino oscillations, the integrated effect can induce a non-trivial ratio between LNV and LNC processes at colliders (for other early papers calculating non-trivial LNV/LNC ratios, without specifying oscillation formulae, see e.g. [13,14]). The parameter region in which such non-trivial ratios are expected have been discussed e.g. in [15,16]. In [15] it has been shown, using the oscillation formulae from plane wave arguments (cf. [12]), that the signature of heavy neutrino-antineutrino oscillations can be resolved at collider experiments, considering minimal low scale type I seesaw extensions of the SM.
It has been shown in [15] that in the case of the low scale minimal linear seesaw model and inverse light neutrino mass hierarchy, the light neutrino mass splittings predict the heavy neutrino mass splitting, resulting in an oscillation length of order O(10 cm) that could be resolved e.g. at the (high-luminosity phase of) LHCb. Furthermore, it has been pointed out that in a realistic experimental setting, where the momenta of the heavy neutrinos are given by a distribution, reconstructing these momenta and considering the oscillations as a function of the heavy (anti)neutrino proper time is required to resolve the oscillation patterns. The method has been demonstrated for an example parameter point, consistent with the present searches and non-collider constraints [15].
Currently there is an increasing interest in exploring heavy neutrino-antineutrino oscillations. For example, oscillations for rare W boson decays at LHC are further studied in [17,18], and from tau decays in [19,20] (using the formulae of [11]). The authors of [21] use plane wave arguments (together with a discussion of coherence conditions) to discuss the oscillations for a rather general right-handed neutrino mass matrix in a left-right symmetric extension of the SM. Discussions of how the insight into the low scale seesaw parameters from heavy neutrino-antineutrino oscillations can help to test whether the baryon asymmetry of the universe can be produced by the leptogenesis mechanism, are given e.g. in [22,23]. So far, no QFT treatment using external wave packets has been performed for heavy neutrino-antineutrino oscillations yet.
The goal of this paper is to put the discussion of heavy neutrino-antineutrino oscillations on a more solid theoretical ground by treating them in the framework of QFT. We use the formalism of external wave packets (cf. [8]) to obtain a more fundamental derivation of the formulae for heavy neutrino-antineutrino oscillations, as well as a more fundamental discussion of the coherence and localisation conditions that have to be satisfied such that oscillations are observable. In addition, we apply our formulae to a specific realisation of the symmetry protected seesaw scenario (SPSS) [33,34] and expand them in small lepton number violating parameters. The paper is organised as follows: In section 2 the formalism of heavy neutrino-antineutrino oscillations is described, the oscillation formulae as well as formulae for the expected number of events are derived and the coherence and localisation conditions are discussed. Section 3 contains the applications to the specific SPSS model and section 4 the approximations of the oscillation formulae and the discussion of the LO and NLO effects. In section 5 we conclude.

QFT Formalism for Heavy Neutrino-Antineutrino Oscillations
As mentioned above, we define neutrinos (antineutrinos) as those particles produced from the decay of a W boson, together with a charged antilepton (lepton), respectively. As a specific example we consider in the following the dilepton-dijet signature at pp colliders, which can be LNC as well as LNV. The relevant Feynman diagrams, in which the W boson is produced from pp collisions, are shown in figures Fig.1 a and Fig.1 b. With a W boson decaying into an antilepton, the produced neutrino state is a superposition of neutrino mass eigenstates, where we only consider the heavy mass eigenstates. We define this projection of the (anti)neutrino interaction eigenstate onto the heavy neutrino mass eigenstates as "heavy neutrino interaction eigenstate". The neutrinos then propagate over a macroscopic distance after which they decay into either a lepton or an antilepton and an off-shell W boson, which in turn decays into two jets. If the neutrino superposition decays into a lepton and an off-shell W boson, the process is lepton number conserving (LNC) and if it decays into an antilepton it is lepton number violating (LNV). Although we focus on this example process, our results can be readily adapted to other processes as well 2 . Fig.1 a) Feynman diagram for the LNC process Fig.1

b) Feynman diagram for the LNV process
Starting with the general formula for the connected amplitude where H I is the interaction Hamiltonian andT the time ordering operator, we follow the procedure described in [8] for the LNC and LNV processes separately.
Using standard QFT methods in the canonical quantisation formalism, the process 2 The process in which the initial W + is replaced by a W − differs only in the leptonic mixing factors (later denoted by V ), which are complex conjugate to the ones appearing in the process with an initial W + . described by the Feynman diagrams figures Fig.1 a and Fig.1 b is obtained by expanding the exponential up to second order in the electroweak coupling constant. The external states are considered to be wave packets that are centred at the space-time points of production or detection as in [8]. This introduces integrals over the momenta of the external particles together with functions describing the shape of the wave packets, which we denote as Ψ, as well as space-time translation operators in the form of exponentials accompanying each wave packet. The propagator of the intermediate particle is obtained by contraction of the relevant fields in configuration space. Using the Fourier transformation it is written down in momentum space, introducing the integral over the neutrino momentum p. The relevant Feynman rules, which are applicable also to the LNV process, can be found in [24]. The amplitudes for the processes described by the Feynman diagrams figures Fig.1 a and Fig.1 b can be written as where α and β denote the flavour indices of the charged leptons at production and detection and LNC and LNV refers to the lepton number violating or lepton number conserving process, respectively. The sum runs over the propagating mass eigenstates, which are denoted by the indices i and j in this paper. The relevant part of the lepton mixing matrix, rotating the heavy neutrino mass eigenstates into the active neutrino interaction eigenstates is denoted by V and the so-called partial amplitude A LN X j is defined as where Ψ(k, K) describes a wave packet centred at three momentum K. The integration measure for the three momenta is written in a short hand notation where with E k being the energy of the respective particle. The interaction amplitude M LN X j is defined as the matrix element of the LNX process without the lepton mixing matrix elements and without the denominator of the propagator, where LNX can refer either to LNC or LNV. Note that we have suppressed the spin and polarisation labels of the external particles to simplify notation. The production (x 0 ) and detection (x 1 ) points in spacetime are defined with respect to the laboratory frame. Therefore the propagation distance, which is defined as L := x 1 − x 0 , is also to be understood in the laboratory frame. The propagation time is given by T : The authors of [8] proceed with the partial amplitude by evaluating the integrals over the three momenta of the external particles, which can be done analytically if the wave packets are assumed to have a Gaussian shape and if the interaction amplitude is approximated at the mean momenta of the external particles. The approximated interaction amplitude is written as M LN X j (p, Q, P 1 , X N ) := M LN X j (p, Q, P 1 , P 2 , K 1 , K 2 ), where X N denotes the mean momenta of the decay products of the heavy neutrinos. This is followed by the integration over x and x corresponding to the production and detection vertices.
The last step in the computation of the partial amplitude consists in the integration over the four momentum of the intermediate particle, which shows up in the propagator. The energy integral is done by the use of the Jacob-Sachs theorem [9], which basically puts the intermediate particle on-shell. The integration over the three momentum is done in three regimes separately depending on the propagation distance of the intermediate particle. The following steps are valid in the longitudinal dispersion regime, which is applicable when the propagation distance is larger than the dispersion length (cf. [8]). In section 2.1 it is argued that this regime is indeed the relevant one when we estimate the widths of the wave packets from measurement uncertainties. In this regime the three-momentum of the j-th neutrino mass eigenstate is approximated around the point of stationary phase p cl,j , which is the "classical" momentum given by where the classical velocity is defined as v cl = L/T with the corresponding gamma factor γ cl = (1 − |v cl | 2 ) −1/2 . The classical energy is defined accordingly as In the LNC case, the interaction amplitude M LN X j is dependent on p, which we approximate at the point p cl,j as well.
The "oscillation probability densities" are proportional to the absolute value squared of the respective amplitudes, averaged over the macroscopic propagation time To proceed with this integration Laplace's method is used in which the macroscopic propagation time is expanded aroundT 0 (cf. [8]) where p 0 := p 0 · L/|L|. The mean energy is defined asẼ 0 := m 2 0 + |p 0 | 2 , with the arithmetic mean of the heavy neutrino massesm 0 . The four momentum p 0 is defined via the mean momenta of the external particles using energy-momentum conservation at the production and/or detection vertex, which yields A related mass can be defined as Note that these so called "reconstructed" quantities, which are labelled with an subscript ( ) 0 , represent experimentally accessible quantities. As described below (see section 2.1), the reconstructed mass m 0 is related to the physical masses of the heavy neutrinos. In particular we note that if the heavy neutrinos are almost mass degenerate such that m i ≈m 0 , eq. (50) can be used to argue thatm 0 ≈ m 0 within the momentum uncertainty given by the wave packets widths. Using eq. (8) yields Due to the wave packet nature of the intermediate particle there is a non-zero probability to measure the decay vertex in a direction L from the production vertex not parallel to p 0 . However, those orthogonal directions are exponentially suppressed and negligible if the momentum uncertainties are small, i.e. if σ p (see appendix A) is smaller than the orthogonal momentum p 0 × L/|L|. We therefore introduce an integration over the direction of L, which is evaluated by approximating L/|L| ≈ p 0 /|p 0 |. This approximation allows to identify p 0 = |p 0 |. Also, with m j ≈ m 0 , which holds for nearly mass degenerate heavy neutrinos as mentioned above, the four momenta p cl,i and p 0 are approximately equal.
The time integration together with the integration over the direction of L leads to where L := |L| has been defined. We note that in the case of the no-dispersion regime, following the steps in [8], one obtains the same equation eq. (12) and therefore the next considerations (until section 2.1) hold for both regimes. The proportionality constant N 2 g can be obtained by the following normalisation condition, which has to be computed ∀{L, α} where dX N denotes an integral over the whole phase space of the decay products of the heavy neutrino PS N and spins represents the sum over all outgoing spins/polarisations and the average over all incoming spins/polarisations. Lepton mixing matrix factors are contained in V LN X αβ ij , which is defined as The oscillation length is given by where p 0 is defined in eq. (9). Additional terms which can be neglected, given the adequate kinematic and experimental conditions, are discussed in section 2.1.
Together with the normalisation condition, the oscillation probability densities P LN X αβ (L, Q, P 1 , X N ) are defined as densities with respect to the mean momenta of the decay products of the heavy neutrino. If this density is integrated over the considered phase space of the heavy neutrino decay products P S N p ⊂ P S N an oscillation probability is obtained as The results describe the probabilities that the superposition of heavy neutrino mass eigenstates, produced by the decay of a W boson together with an antilepton of flavour α, produces an (anti)lepton of flavour β if it decays after a distance L in the direction of p 0 via an LNC (LNV) process 3 . When summing these probabilities over the flavour of the final (anti)lepton, the resulting quantity β P LN V αβ (L, Q, P 1 ) can be interpreted as the probability that the produced heavy neutrino interaction eigenstate has oscillated into a heavy antineutrino interaction eigenstate. The quantity β P LN C αβ (L, Q, P 1 ) can be interpreted as the probability that the heavy neutrino interaction eigenstate has "survived", i.e. has not oscillated into a heavy antineutrino interaction eigenstate.
At this point eq. (16) together with the additional terms, to be discussed in section 2.1, can be regarded as the most general result. We now proceed simplifying it, in order to gain further insight. To this end, we first show under which conditions the interaction amplitudes can be factored out of the sum over mass eigenstates, and subsequently be absorbed into the normalisation constant. 4 After that we commit to a specific model within the SPSS, which features two almost degenerate heavy neutrinos.
The interaction amplitude is dependent on the masses of the propagating neutrinos through the numerator of the propagator, which reads ( / p cl,j + m j ). The neutrino masses are expressed as a deviation from the mean neutrino massm 0 . In the case of just two heavy neutrinos 5 , with masses m 4 and m 5 , where w.l.o.g. m 4 < m 5 can be assumed, one finds that where the dimensionless mass splitting parameter λ m is defined as and the mean mass is just given bỹ Using the definition eq. (5), the classical momentum can also be expressed using the mass splitting parameter, which yields where the mean momentum is defined asp cl,0 =m 0 γ cl v cl . Note that using eqs. (17) and (20) it is easy to reparameterize the four momentum of an on-shell particle. This makes it possible to factor out the mass dependence from the interaction amplitudes, yielding where Λ ij contains the factors describing the mass and momentum splitting and is given by The interaction amplitudes of the processes figures Fig.1 a and Fig.1 b can be written down by using the Feynman rule conventions described in [24], which are applicable also to the lepton number violating diagram. This yields and where g 2 is the coupling of the SU (2) gauge bosons and G F is the Fermi constant. The vertices Γ µ and Γ µ are given by γ µ P L and −γ µ P R respectively, with the left and right chirality projection operators defined as Note that, as above we have suppressed spinor and color indices as well as indices denoting the spin and polarisation of the external particles. For the neutrino we chose the fermion flow from left to right in figures Fig.1 a and Fig.1 Further simplifications are possible if the spin correlation between the production and detection vertex are neglected, i.e. if the numerator of the propagator can be written as This approximation is also done in the narrow width approximation and makes it possible to factorize the interaction amplitude into a production interaction amplitude and a detection interaction amplitude. Using this approximation, the interaction amplitudes for the LNC and LNV process are identical.
If the spin correlation is not neglected, the interaction amplitudes of the LNC and LNV process differ due to the chirality structure. For a given process it could be possible that the interaction amplitudes depend on the orientation of the momenta of the external particles in such a way that a probabilistic classification into LNC or LNV becomes possible, which could be interesting, e.g. for the SHiP experiment [25].
In order to simplify the expression eq. (12), following the above discussion, the spin correlation between the production and detection vertices are neglected. This allows to absorb the mass splitting independent parts of the interaction amplitudes in eq. (12) into the normalisation constant.
This leads to the following oscillation probability, which is independent of the mean momenta of the decay products of the heavy neutrino, spins and polarisations of the external particles Due to these simplifications the oscillation probability only depends on Q and P 1 in the combination |p 0 | = |Q − P 1 |, which yields The normalisation condition for the simplified oscillation probability is given by In section 4 the normalisation constant is evaluated explicitly for a specific example model of interest.
To proceed further we assume the experimental conditions and model parameters to be such that the heavy neutrinos travel a macroscopic distance before they decay, forming a displaced vertex. The number of expected events that feature such a displaced vertex can be expressed as a formula similar to the one described in [26]. This formula has to be modified in order to cover the circumstances of this paper. In particular an expression for the probability that the heavy neutrino decays in an LNC (LNV) manner involving a specific lepton flavour between a minimum and maximum distance is needed. Following the discussion in [8] regarding unstable oscillating particles, it is shown that the only relevant modification to the oscillation formula is given by the exponential discussed in section 2.1. If the masses and decay widths of the heavy neutrinos are nearly identical, it is possible to define a common decay length asL where the common decay width is defined as The exponential describes the probability that a particle is still present at distance L. Therefore the probability density that the particle decays at distance L is given by the derivative With this the probability that a particle decays in an LNX process into flavour β between x min (ϑ) and x max (ϑ) is given by where the subscript dv stands for displaced vertex and ϑ denotes the angle of the heavy neutrino with respect to the beam axis. Usually the interval [x min (ϑ), x max (ϑ)] will be chosen to lie inside the detector, such that the decay products can be measured. The detector geometry can be taken into account by the dependence on ϑ. The number of expected LNX events in which an antilepton of flavour α is measured at the production vertex and a lepton (LNC) or antilepton (LNV) of flavour β is measured at the detection vertex is given by whereσ N,0 is the mean production cross section of the heavy neutrinos, which depends on model parameters such as the masses of the heavy neutrinos and the details of the lepton mixing matrix.Br ljj,0 is the mean branching ratio for the decay of a heavy neutrino into a lepton and two jets and L is the time integrated luminosity. Therefore, the factorσ N,0B r ljj,0 L describes the total number of events in which a heavy neutrino is produced and decays into a lepton and two jets. It is assumed that the branching ratios as well as the production cross sections for different mass eigenstates are nearly identical in order for such an approximation to be appropriate. The remaining factor D N (ϑ, |p 0 |) accounts for the probability density that the reconstructed heavy neutrino has momentum of modulus |p 0 | and is produced with an angle ϑ with respect to the beam axis. As discussed above, P LN X dv αβ (x min (ϑ), x max (ϑ), |p 0 |) gives the probability that the heavy neutrino decays in an LNX manner into flavour β inside the interval [x min (ϑ), x max (ϑ)]. Note that if one wanted to consider the spin correlation, one has to use eq. (16) in the definition of eq. (34).

Observability Conditions and Dispersion Length
In order for the oscillations to be observable there are several conditions which have to be satisfied. This subsection describes those conditions and estimates their viability for typical parameters of long-lived heavy neutrinos detectable at, e.g. HL-LHCb. We consider as an explicit example a parameter point for the minimal low scale linear seesaw model that has also been used in Ref. [15].
To compute all relevant parameters it is necessary to know the kinematics of the process and the widths of the external wave packets. One approach to estimate these widths is based on the following considerations: Since the final particles at production and detection are reconstructed by measurements at a detector, the uncertainty of the measurement should be reflected by the widths of the respective wave packets. For charged leptons a relative momentum uncertainty in the range (∆p/p) lepton = [0.5%, 1%] holds for particles with a long enough track, cf. [27], which we therefore use for the widths of the charged lepton wave packets. The momentum resolution for the quarks, which are reconstructed from displaced jets, is much harder to determine. For a conservative estimate we therefore use a large range for their relative momentum uncertainty (∆p/p) quark = [5%, 30%] (and thus for the possible widths of the quark wave packets). Finally, the width of the wave packet of the initial W boson is taken to be its decay width (Γ W ≈ 2GeV).
Alternatively, one can also try to estimate the widths of the wave packets in position space based on the consideration that the uncertainty is determined by interactions with detector/beam particles and their respective widths. For the W boson one could use the proton-proton distance in the proton bunches, which at the LHC is about 5 × 10 −6 cm, whereas for the leptons and quarks one might take a wave packet width of the order of an atom radius, i.e. about 10 −8 cm. This would lead to significantly smaller widths in momentum space compared to the estimate using the measurement uncertainty, such that the appropriate regime is the no-dispersion regime (cf. [8]). We have checked that the relevant observability conditions in this case are all satisfied for our example parameter point, and our results from eqs. (16), (35), (73) and (74) can also be used for these estimates of the wave packet widths. From now on we will focus on the estimates for the wave packet widths from the momentum space considerations. Table 1: Parameters used for simulating the kinematics of the processes in figures Fig.1 a and Fig.1 b, in order to evaluate the observability conditions for the example minimal low scale seesaw parameter point used in [15]. γ denotes the gamma factor of the heavy neutrinos,m 0 their mean mass (which in the simulation coincides with m 0 ), and the squared mass splitting δm 2 45 = m 2 4 − m 2 5 in the scenario of [15] is predicted by the measured values of the light neutrino mass splittings. The used range for the uncertainties in the measurement of the momenta of the external charged leptons and jets are denoted as (∆p/p) lepton and (∆p/p) quark .
The kinematics of the processes described by figures Fig.1 a and Fig.1 b have been simulated, assuming two nearly degenerate heavy neutrinos. Since the purpose of the simulation is to compute the parameters necessary to check the observability conditions and not to simulate the oscillation process itself, all particles can be treated as plane waves, where their momentum represents the momentum of the peak of the wave packet. Simulating enough events, where the momentum of the W boson is taken in the range 340 GeV to 2 TeV. It is then possible to check if the observability conditions are fulfilled. For the simulation the parameter values in table 1 have been used. It has also been assumed for simplicity that there are two heavy neutrino mass eigenstates with masses m 4 and m 5 .
The observability conditions are given as exponential suppression factors. If it is not clear that they are satisfied those exponential factors have to be included into the probability eq. (12) or eq. (28), respectively. It is worth mentioning that including exponential terms into the probability changes the normalisation constant, which can be computed using eq. (13) or eq. (29). The quantities used in the computation of the observability conditions are defined in appendix A.
The oscillation length sets the scale of the experiment, since it is the length at which the measurements should be taken in order to observe oscillations. As stated above it is given by With the parameters in table 1 (with γ = 50) it can be computed to be where the two heavy neutrino mass eigenstates are labelled with subscript numbers 4 and 5.
The effective width σ pef f (see [8]) can be interpreted as the width of the wave packets of the heavy neutrinos. A rough estimate of the effective width can be obtained by neglecting the detection process and by approximating the lepton in the production process as a plane wave. Due to energy-momentum conservation, the shape of the effective wave packets of the heavy neutrinos is then given by the one of the initial W boson, which width is approximated by its decay width. However, simulating the kinematics of the process and computing the width numerically shows that it is in the range showing that the estimate σ pef f ≈ 2 GeV is indeed very rough.
The dispersion length is the threshold at which the spread of the wave packets of the heavy neutrinos becomes significant in all directions. At distances larger than the dispersion length the methods of the longitudinal dispersion regime have to be used. The dispersion length is given by [8] L where v 0 = |p 0 |/E 0 with p 0 and E 0 being the momentum and the energy given by the mean momenta of the external particles and energy-momentum conservation at either vertex as defined above. The simulation using the parameter values from table 1 shows that L osc > 1000 L disp . The assumption that the longitudinal dispersion regime is the relevant one, in the case in which the wave packet widths are estimated from measurement uncertainties, is therefore well justified.
The observability conditions can be divided into two groups. The so-called coherent effects, which are taken into account at the level of the wave function and the so-called incoherent effects, which are included at the level of the probability.

Coherent Effects
The so-called coherence length describes the decoherence of the wave packets, which can have two origins. The oscillations either vanish if the wave packets become separated due to different group velocities, or if the wave packets spread beyond the oscillation length, in which case the oscillations are averaged to zero. In momentum space both of these effects are taken into account by the exponential where (see [8]) These terms can be neglected if the momentum of the intermediate particle is much larger than the width of its effective wave packet. Using the parameters in table 1, all simulated events satisfy at least L coh 45 > 10 L osc 45 .
It is therefore justified to neglect the effects of decoherence in the first oscillation cycles for the parameter values in table 1.
We remark that the oscillation length and therefore the coherence length for oscillations including both the light and heavy neutrino mass eigenstates is smaller than 10 −12 m. Therefore it is appropriate to neglect the light neutrino mass eigenstates in the oscillations.
Localisation conditions determine whether there is decoherence from the start. The relevant exponential suppressing the oscillations reads (see [8]) where δm 2 ij := m 2 i − m 2 j . The effective width of the propagating neutrino is given in momentum space by σ pef f . The parameters σ m and ρ are determined by the widths of the external particles and their velocities (see [8]). Using the parameters from table 1 it follows that is satisfied in all events. This allows to neglect the effects from localisation in the oscillation formula for parameter values as in table 1.
In the process considered in this paper the heavy neutrinos are unstable intermediate particles. Therefore the full propagator should be used in eq. (4). As discussed in [8] this leads to the exponential decrease factor where the decay length is given by that has to be included in the sum over the mass eigenstates in eq. (28). With the definitionsm and the decay exponential can be written as where for i = j = 4 the plus sign and for i = j = 5 the minus sign is used. The first exponential can be absorbed into the normalisation constant using the condition eq. (29), since it does not depend on the mass indices i and j. If e.g. the decay widths of the mass eigenstates are too different, the second exponential can lead to a suppression of the oscillation pattern. On the other hand, the exponential is negligible if the mass eigenstates are nearly degenerate and if the decay widths are nearly equal. For the example point considered in this paper one can verify that L(δmΓ 0 + δΓm 0 ) p 0 such that this conditions is satisfied.
Since the detection of the decay products allows to reconstruct the invariant mass of the propagating particle, a condition that relates the mass of the propagating particle with the detection uncertainty is expected. This condition stems from the exponential that appears in the derivation of the amplitude (see [8]). Here where m i and m j are the masses of the heavy neutrinos and m 0 is as given above 6 The parameter σ m (see [8]) is related to the widths and velocities of the external particles. Since the detection process is described by the interaction with those external particles, their widths are in turn related to the precision of the momentum measurement. In conclusion this condition enforces the neutrinos to be either nearly degenerate in mass or highly relativistic, such that the propagating mass eigenstates are within the uncertainty of the momentum measurement. The simulation shows that σ m ∈ [0.0035 GeV, 32 GeV], such that the above condition is satisfied for the parameters in table 1.
In the processes considered in this paper (figures Fig.1 a and Fig.1 b) a W boson is decaying in flight and is therefore an unstable source. In principle one would have to treat the W boson as a propagator connecting the diagrams figures Fig.1 a and Fig.1 b and the particles producing the W boson. This extra propagator would however result in technical difficulties. The authors of [28][29][30] have used perturbation theory in a quantum mechanical model to describe neutrino oscillations and found an additional localisation condition, which suppresses oscillations if the unstable source moves a distance greater than the oscillation length during its lifetime. The unstable source has been assumed to have a mean momentum at rest in those derivations. A QFT approach to light neutrino oscillations has been used in [31], where a similar localisation condition to the one above has been derived. The authors of [32] considered a moving unstable source, i.e. a pion decaying in flight, in a QFT treatment. They obtained the constraint where v π , |p π |, m π , Γ π are the velocity, momentum, mass and decay width of the decaying pion. In our case the initial W boson takes the place of the pion as the unstable source. The heavy neutrino is highly boosted such that the velocities v W and v 0 are almost parallel. Furthermore it holds that v W v 0 for the parameter space considered in this paper, which implies that Even putting the momentum of the W boson as 2 TeV, which is the maximum of the range considered in this paper, we find that This constraint is therefore negligible for the process and parameter space considered.
As an additional remark, note that the above mentioned localisation conditions describe the suppression stemming from the fact that the production vertex is not known due to the finite lifetime and movement of the unstable source. At current and considered future colliders the W boson decays promptly, due to its large decay width, which makes the possible decay region much smaller than any macroscopic oscillation length. Already this argument shows that the additional localisation condition, due to the instability of the W, should be satisfied.
As discussed in [30], the unstable source can also lead to a loss of coherence for distances larger than the coherence length For the parameters in table 1, this additional coherence length can be neglected compared to the one given in eq. (40).

Incoherent effects
If the propagation distance is not precisely known, which is the case if the production or detection points are measured with some uncertainty, neutrinos that have travelled different distances overlap and wash out the oscillation pattern. This effect can be described by the following exponential (see [8]) where we have assumed that the propagation distance of the neutrinos is given by a Gaussian with width ∆L. Oscillations vanish if ∆L ≥ L osc ij . In our case the oscillation length is around 8cm and therefore much bigger than the uncertainty in the resolution of the position of the primary and secondary vertex.
In a real experiment there is a distribution of mean momenta of the external particles, such that the reconstructed momentum of the intermediate particle, denoted by p 0 , follows a distribution as well. The effect is already included in eq. (35), where the factor D N (ϑ, |p 0 |) describes the distribution of |p 0 |. A distribution of |p 0 | leads to a washout of the oscillation pattern, since different oscillation lengths superimpose. In order to resolve the oscillation patterns it is therefore helpful [15] to plot the oscillation probability as a function of the reconstructed proper time using the following relations where the reconstructed gamma factor is defined as γ 0 = ( 1 − |v 0 | 2 ) −1 , the reconstructed velocity is as before v 0 = p 0 /E 0 , the reconstructed time is given by T 0 = L/|v 0 | and the reconstructed proper time is given by τ 0 = T 0 /γ 0 . This leads to the following oscillation exponential where eqs. (19) and (47) have been used. Using an oscillation probability based on plane wave arguments, the above method has been demonstrated for an example parameter point (assumingm 0 = m 0 ) in [15]. Note that, as mentioned above, the quantities denoted by a subscript ( ) 0 are the ones which are reconstructed by experimental measurements of the external particles in the process.

Low Scale Seesaw with Symmetry Protection
To further develop and apply the above results, we consider SPSS models (cf. [33,34]) i.e. low scale seesaw models where the smallness of the light neutrino masses is protected by a slightly broken "lepton number"-like symmetry. As a particular example we focus on the "minimal low scale linear seesaw" model with only two righthanded (sterile) neutrinos, that has also been discussed as an example in [15]. The Lagrangian of this model takes the following form where L SM is the Standard Model (SM) Lagrangian, α = (e, µ, τ ) is a family index andφ := φ * with the Levi Civita symbol and the SM Higgs doublet φ. The Yukawa couplings to the sterile neutrinos are denoted by Y α and Y α . In the symmetry limit of the model, the Yukawa couplings Y α are zero, and the "lepton number"-like symmetry is only broken slightly by the Yukawa coupling Y α , for which we assume that Y α Y β for all entries. Possibilities to realise a low scale linear seesaw in the context of SO(10) Grand Unified Theories have been discussed e.g. in [35,36].
After electroweak symmetry breaking the part of the Lagrangian responsible for neutrino masses and mixing can be written as where n = (ν e L , ν µ L , ν τ L , The symbols m and m denote column vectors given by Y v EW / √ 2 and Y v EW / √ 2, respectively. v EW ≈ 246 GeV is the electroweak vacuum expectation value. The mass matrix can be diagonalised using a Takagi decomposition This can be achieved following the steps in [37] where first a block diagonalisation with an exponential ansatz is performed, followed by a diagonalisation of the active neutrino 3x3 and sterile neutrino 2x2 block. Expanding the exponential to second order yields where θ = m * (M * ν h ) −1 . It can be easily checked that the mass matrix is indeed block diagonalized to second order in θ.

Furthermore, one finds that
diagonalizes the heavy neutrino 2x2 block. Regarding the heavy neutrino-antineutrino oscillations, the interesting part of the mixing matrix is the upper right 3 × 2 block, which describes the mixing between the active neutrino interaction eigenstates and the heavy neutrino mass eigenstates. In the oscillation formula we called this part V . From eqs. (63) and (64) we find that where the phase is defined as 2φ = arg ( θ ) · ( θ) * and the expansion parameters are defined as Note that θ can be seen as a function of θ and θ . From the hierarchy between the Yukawa couplings of the sterile neutrinos it follows that θ θ. We restrict ourselves to maximally first order in the elements of θ in the following.
As a remark one can see that higher order terms in θ can be absorbed in a slight rescaling of the Yukawa couplings Y , since in the symmetry limit (where θ = 0) V has the exact form Therefore those terms do not qualitatively change the oscillation formulae that are obtained in section 4.
Using the expansion eq. (65) it is also possible to obtain expressions for the masses of the heavy neutrinos and thus also for the mass splitting parameter λ m , introduced in eq. (18), which can be expressed as Another important parameter entering the oscillation formula is the quadratic mass splitting δm 2 45 , which enters the formula for the oscillation length. It can be expressed as where the mean mass of the heavy neutrinosm 0 is defined in eq. (19) and can be expressed asm

Approximations of the Oscillation Formula
Using the results from section 3 the oscillation probability eq. (28) can be expanded in the small parameters θ appearing in the lepton mixing matrix and λ m which has been introduced in eq. (18) to account for the mass splitting of the heavy neutrinos. Also, as we discussed, the mass splitting parameter λ m can be neglected at leading order as it is much smaller than O( θ ), see eq. (69). The following definitions are used To derive the approximate oscillation formulae, the normalisation constant is computed using the condition eq. (29), yielding where the ellipses contain orders higher then O( θ ). With this the probability eq. (28) can be expanded in θ . Note that we keep the exponential exp (iΦ ij L) exact. Up to first order in θ , the probability eq. (28) in the LNC case yields In the LNV case the expansion of eq. (28) up to first order in θ yields where in both cases the LO terms are written in the first line and the NLO terms in the second line. Note that if the initial W boson is replaced by a W − , the leptonic mixing matrix factors are complex conjugate to the ones in the process where the initial boson is a W + . This results in a sign change of the NLO contributions. The leading order term in those expansions describe the oscillations from neutrinos into antineutrinos, whereas the first order term describes flavour oscillations. This can be seen by adding up the LNC and LNV probabilities, which means that the sign of the outgoing lepton is ignored. This should make the oscillations of neutrinos into antineutrinos vanish, and indeed the oscillatory part of the leading order terms cancel each other. We are left with Summing over the outgoing flavours makes the oscillatory part of the above equation vanish. This happens because β θ * β θ β exp(−2iΦ) ∈ R, and therefore β I β = 0.
With the "lepton number"-like symmetry being broken, one also expects lepton number violation in the limit where the distance L goes to zero. As mentioned above, the no-dispersion regime, which is the relevant one in this limit, results in the same formulae for the oscillation probabilities if the observability conditions are met. From eq. (74) it can therefore be seen that there is no lepton number violation at zero distance. This leads to the conclusion that this effect has to be introduced at a higher order and is therefore much smaller than the lepton number violation due to oscillations. That this effect is indeed present at higher orders can be confirmed by numerically diagonalizing the lepton mass matrix and using eq. (28) to compute the oscillation probability.
Taking only the leading order into account, the only relevant model parameters are the Yukawa couplings Y α , or equivalently the mixing parameters θ α , and the quadratic mass splitting δm 2 45 appearing in φ 45 . The leading order effects can therefore be described by the symmetry limit of the SPSS [33,34] plus the quadratic mass splitting as an additional parameter. This holds even more generally, for all realisations of the SPSS.
If the mechanism of light neutrino masses is given by the minimal linear seesaw described in eq. (59), one can reparameterize the model in terms of active neutrino parameters according to [38]. Assuming an inverse ordering of the light neutrino masses m ν i yields where U denotes the unitary PMNS matrix, and Note that we have absorbed the parameter appearing in [38] into the definitions of Y and y . An interesting observation is that the heavy neutrino mass splitting is given by the light neutrino mass splitting [15] and therefore the squared mass splitting is given by Taking the values of the active neutrino mixing angles and mass squared differences from [39,40], the only undetermined parameter is the Majorana phase α. Analyzing the parameterization of the Yukawa couplings one finds that the products  Fig.2 a, Fig.2 b, Fig.3 a and Fig.3 b. The values of the oscillation probabilities makes it clear that any higher order effects on the oscillation patters will be extremely challenging to observe under realistic conditions for the chosen example parameters from table 1.

Summary and Conclusions
In this paper we have applied the framework of quantum field theory (QFT) with external wave packets (cf. [8]) to derive the probabilities for the oscillations between long-lived heavy neutrino and antineutrino interaction eigenstates (cf. footnote 1), where we define a neutrino (antineutrino) as the neutral lepton that is produced together with a charged antilepton (lepton) and a W boson. These heavy neutrinoantineutrino oscillations can lead to an oscillating rate of lepton number conserving (LNC) and violating (LNV) events at colliders, as a function of the distance between the (anti)neutrino production and displaced decay vertices.
Our most general formula for the oscillation probability is eq. (16) together with the additional terms discussed in section 2.1. The latter can be neglected given the adequate kinematic and experimental conditions, and are referred to as observability conditions. The oscillation probability densities can be further simplified to eq. (28) by neglecting the spin correlation between the production and detection vertex. Including the decay probabilities of the heavy neutrinos, formulae for the expected number of LNC/LNV events with a certain displacement between primary and secondary vertex have been given in eq. (35).
The simplified formulae for the oscillation probabilities have been applied to low scale seesaw models where the smallness of the light neutrino masses is protected by a slightly broken "lepton number"-like symmetry, i.e. to the SPSS (cf. [33,34]). As a particular example we have focused on the "minimal low scale linear seesaw" model with only two nearly mass-degenerate heavy neutrinos, that has also been discussed as an example in [15]. Within this class of models, an expansion of the probabilities in the small "lepton number"-like symmetry breaking parameters θ has been performed, yielding the LO and NLO contributions (cf. eqs. (73) and (74)).
For the example parameter point used in Ref. [15], we have discussed the observability conditions (cf. section 2.1) and found that they are all satisfied. However, if the momentum p 0 is given by a distribution, which is the case if the mean momenta of the external particles follow a distribution as discussed in section 2.1, the oscillation pattern can be washed out as has already been pointed out in [15]. The proposed solution to this is to reconstruct the four momentum p 0 (from the measurements) and to consider the oscillations as a function of the heavy neutrino reconstructed proper time. In [15] it has been demonstrated, using estimated uncertainties for the HL-LHCb and the above-mentioned example parameter point, that the proposed solution is indeed feasible.
Comparing our simplified LO results with the existing literature, we found that we agree with the results from [11] (when we set their parameters θ LN V 21 = −π and θ LN C 21 = 0 to match the considered low scale seesaw scenario). Our LO formulae also agree with the ones derived using the formalism for meson oscillations and plane wave arguments and used e.g. in [12,15]. Our results in the most general form, i.e. eq. (16) together with the additional terms discussed in section 2.1, allow to discuss effects beyond LO and to check the observability conditions (or include them explicitly in the calculations).
Our NLO results showed that beyond the LO heavy neutrino-antineutrino oscillations, the probablities are also modulated by "flavour oscillations", as discussed in section 4. On the other hand, for the case of the "minimal low scale linear seesaw" model (with parameters around the considered example point), it has turned out that the NLO effects are very small, with a suppression which makes them undetectable at the currently considered future collider experiments. While this does not necessarily have to be the case for other choices of parameters, it indicates that there is a parameter region of interest for the LHC (and future colliders) where the LO formulae are sufficient. In this region the only model parameters relevant for heavy neutrino-antineutrino oscillations, in terms of the SPSS parameters, are the three flavour-dependent active-sterile mixing angles θ α and the mass squared difference δm 2 45 between the two heavy neutrinos.
In summary, our results show that the phenomenon of heavy neutrino-antineutrino oscillations can indeed occur in low scale seesaw scenarios and that the previously used leading order formulae, derived with a plane wave approach, provide a good approximation for (at least) the considered example parameter point. Our results help to put existing studies based on LO formulae on a more solid theoretical ground (by providing the observability conditions which have to be checked in the QFT framework) and can be used in future studies to explore the phenomenon in other parameter regions and for different types of low scale seesaw models.

A Formulas for the Observability Conditions
In the following, the formulas to compute the observability conditions of section 2.1 are given. For more details we refer to [8].
The particles at production are labeled P i and their wave packets are assumed to have a Gaussian form of width σ xP i in configuration space. The wave packet peaks, in momentum space, at momentum P i and since they are assumed to be on-shell, their peak energy is given by E P i = m 2 p i + |P i | 2 . The peak velocity is then defined as v P i = P i /E P i . For the particles at detection the letter P is simply replaced by D. The velocity v 0 is defined using energy-momentum conservation at the production and/or detection vertex, which yields v 0 := p 0 /E 0 .
Labelling the incoming particles at production P i,in and the outgoing ones P i,out yields and The relevant parameters are defined as follows.
The following symbols are defined in the longitudinal dispersion regime, on which we focused in this paper. We denote a velocity (v) projected onto the direction in which the oscillation distance is measured 7 (L) as ν. (93) 7 If one does only measure the scalar distance of oscillation, one might replace the direction L/|L| with the direction p 0 /|p 0 |, which can be interpreted as the main direction in which the heavy neutrinos travel.