Probing the Majorana neutrinos and their CP violation in decays of charged scalar mesons $\pi, K, D, D_s, B, B_c$

Some of the outstanding questions of particle physics today concern the neutrino sector, in particular whether there are more neutrinos than those already known and whether they are Dirac or Majorana particles.There are different ways to explore these issues. In this article we describe neutrino-mediated decays of charged pseudoscalar mesons such as $\pi^{\pm}$, $K^{\pm}$ and $B^{\pm}$, in scenarios where extra neutrinos are heavy and can be on their mass shell. We discuss semileptonic and leptonic decays of such kinds. We investigate possible ways of using these decays in order to distinguish between the Dirac and Majorana character of neutrinos. Further, we argue that there are significant possibilities of detecting CP violation in such decays when there are at least two almost degenerate Majorana neutrinos involved. This latter type of scenario fits well into the known neutrino minimal standard model ($\nu$MSM) which could simultaneously explain the Dark Matter and Baryon Asymmetry of the Universe.

A related issue in neutrinos physics is the absolute mass values of the known neutrinos. While the experimental evidence of neutrino oscillations within the known three flavor states [41][42][43][44][45][46] clearly shows that these particles cannot all be massless, the oscillations are only sensitive to mass differences, not to their absolute values. In contrast, 0νββ decays are sensitive to the absolute mass and may help in their determination, if neutrinos turn out to be Majorana particles. So far the best bounds on the absolute masses of the light neutrinos come from Cosmology m ν 0.23 eV [47].
A pending question is then why light neutrinos are so light, specifically so much lighter than all other Standard Model (SM) fermions. Interesting enough, the existence of such very light neutrinos can be explained via the seesaw mechanism [48][49][50][51][52] where more neutrinos are required and where all of them are, in general, Majorana particles. In the simplest form of this mechanism, the masses of the light neutrinos are ∼ M 2 D /M R ( 1 eV), where M D is an electroweak scale or lower. At the same time, additional neutrinos, usually much heavier (with masses M R 1 TeV) and sterile under electroweak interactions except through small mixing with the SM flavors, are required. This mixing is suppressed as ∼ M D /M R ( 1). Besides the simplest scenario, there are other seesaw scenarios in which the heavy neutrinos may have lower masses, namely near or below 1 TeV [53][54][55][56][57][58][59] and even near the 1 GeV scale or below [15,[60][61][62][63][64][65], and at the same time their mixing with the SM flavors may not be so extremely suppressed as in the original scenarios.
In the first part of this work we discuss lepton number violating (LNV) semileptonic decays of charged pseudoscalar mesons such as K ± and B ± , mediated by a heavy Majorana neutrino on its mass shell, cf. Ref. [26]. The pions, which are the lightest mesons, can have only leptonic decays; we discuss hypothetical leptonic decays that could be mediated by on-shell heavy neutrinos. Such decays could be either lepton number conserving (LNC) or lepton number violating (LNV) when the mediating neutrinos are Majorana particles, cf. Ref. [37], while only LNC decays occur when the neutrino is of Dirac type. We present ways of determining the nature of neutrinos using the differential decay rates of these processes.
Yet another interesting issue in neutrino physics is the possible existence of CP violation in the lepton sector, a phenomenon that could be measured, for example, in neutrino oscillations [66]. Alternatively, as we present in the second part of this work, leptonic CP violation may show in leptonic LNC and LNV decays of charged pions, cf. Ref. [38], or in semileptonic LNV decays of K ± and B ± , cf. Refs. [39,40]. It turns out that such CP violation becomes appreciable and possibly detectable in these decays if the scenario contains at least two on-shell heavy neutrinos that are almost degenerate. Interestingly, this scenario fits well into the neutrino minimal standard model (νMSM) [60,61,[67][68][69][70][71][72] which contains two almost degenerate Majorana neutrinos of mass ∼ 1 GeV and another lighter neutrino of mass ∼ 10 1 keV, besides the three light neutrinos of mass 1 eV. This model can explain simultaneously the existence of neutrino oscillations, dark matter and baryon asymmetry of the Universe. Furthermore, in more general frameworks of low-scale seesaw, baryon asymmetry (but not dark matter) is explained while keeping even larger values of the heavy-light mixing [73] than in the νMSM, and in such frameworks the case of almost degenerate Majorana neutrinos is preferred [74] since it allows larger mixings. CP violation effects in the neutrino sector in scenarios with nearly degenerate heavy neutrinos have also been investigated earlier [75,76] using a more detailed formalism, although it amounts to the same effect described here which is the interference of amplitudes with two slightly different dispersive and absorptive parts in the neutrino self energy.
In Section II we discuss the LNV semileptonic decays of mesons M ± → ± 1 ± 2 M ∓ , mediated by an on-shell Majorana neutrino (which henceforth we call N ) where M is a heavy pseudoscalar (M = K, D, D s , B, B c ), M is a lighter pseudoscalar, and j (j = 1, 2) are charged leptons, and we present there the corresponding branching ratios. In Section III we present the expressions and values of the branching ratios for the LNC and LNV leptonic decays of charged pions mediated by on-shell N sterile neutrino, π ± → e ± N → e ± e ± µ ∓ ν, as well as the differential branching ratio dBr/dE µ for these decays. We discuss the possibilities of detecting such branching ratios and to discern from them the Majorana or Dirac nature of neutrinos. In Section IV we then extend the analysis of the mentioned leptonic and semileptonic decays to a scenario where we have at least two heavy on-shell sterile neutrinos involved (N 1 , N 2 ), and we present an analysis of CP-violating asymmetries A CP ≡ [Γ(M − )−Γ(M + )]/[Γ(M − )+Γ(M + )] for such processes. In Appendices 1 and 2 we present explicit formulas for our LNV semileptonic decays, and in Appendices 4 and 5 explicit formulas needed for the analysis of our LNC and LNV leptonic decays of the charged pion. Appendix 3 contains formulas needed for evaluation of the decay width of the heavy neutrino N , and in Appendix 6 we derive an identity relevant for CP violation asymmetry. In Section V we discuss and summarize our results.

II. LEPTON NUMBER VIOLATING SEMILEPTONIC DECAYS OF SCALAR MESONS
If there is a Majorana sterile neutrino N with mass M N ∼ 1 GeV, its existence could be discerned by detecting semileptonic LNV decays of heavy mesons mediated by on-shell N . Here we will consider such LNV decays M ± → This is in part one of the reasons for the high suppression suffered by all lepton flavour violating processes involving heavy neutrinos.
Factor 1/2! above is the symmetry factor when the two produced leptons are equal; d 3 is the integration differential of the final three-particle phase space The resulting decay width can be written as where k is the mixing factor and Γ ± (XY * ) are the normalized (i.e., without the explicit mixing) contributions from the X channel and the complex-conjugate of the Y channel (X, Y = D, C, where D and C stand for the direct and crossed channels) The expressions for T ± (X)T ± (Y ) * (X, Y = D, C) are given in Appendix 1. T ± (X) is the relevant part of the amplitude in the X channel and forms part of the total decay amplitude T (M ± ), cf. Appendix 1. In Equation (5), notice that subscripts ± for the contributions Γ(DD * ) and Γ(CC * ) are unnecessary because |T + (D)| 2 = |T − (D)| 2 and |T + (C)| 2 = |T − (C)| 2 . P (X) (X = D, C) are the propagator functions of the intermediate neutrino N in the two channels The overall constant K 2 in Equation (7) is Here, f M and f M are the decay constants of M ± and M ∓ , and V QuQ d and V quq d are the corresponding CKM matrix elements. We denote the valence quark content of M + as Q uQd ; of M + as q uqd . When the intermediate neutrino N has such a mass that it is on mass shell, Equation (2), the squares of the propagators (8) are reduced to Dirac delta functions because where p k = p 1 , p 2 for X = D, C. In this on-shell case, the DD * and CC * contributions in Equation (5) are large, and the interference contributions DC * and CD * are negligible in comparison (cf. Ref. [39] for details on this point), leading to when 1 = 2 , we even have Γ(DD * ) = Γ(CC * ). The normalized decay width Γ(DD * ) can be calculated explicitly, and it turns out to be and Γ(CC * ) is obtained by the simple exchange y 1 ↔ y 2 The notations used in Equations (12) and (13) are λ(y 1 , y 2 , y 3 ) = y 2 1 + y 2 2 + y 2 3 − 2y 1 y 2 − 2y 2 y 3 − 2y 3 y 1 (14a) In this expression, N N (M N ) ≡ N N ( = e, µ, τ ) are the effective mixing coefficients; these are numbers ∼ 10 0 -10 1 which depend on the mass M N . In Appendix 3 we write down the relevant formulas for the evaluation of these coefficients. The results of these evaluations are presented in Figure 2, for the case of Majorana and Dirac neutrino N , in the entire neutrino mass interval 0.1 GeV < M N < 6.3 GeV which will be of interest in this work. For further clarifying remarks we refer to Appendix 3. Equations (11) and (12) imply that Γ(M ± → ± 1 ± 2 M ∓ ) is proportional to 1/ K (∝ 1/|B N | 2 ). Hence, we can define a canonical branching ratio Br, being the part of the branching ratio Br(DD * ) ≡ Γ(M ± → ± 1 2 ±M ∓ ; DD * )/Γ(M ± → all) with no explicit or implicit heavy-light mixing factors Use of the expressions (12) and Equations (16) and (17) then gives for the canonical branching ratio the following expression: where the notations (14) and (9) are used. In the limit of massless charged leptons (M 1 = M 2 = 0) this expression becomes simpler Br(y N ; 0, 0, y ) = 3π 16 In the mentioned branching ratios, an often important effect of suppression due to the decay (i.e., nonsurvival) probability was not included. Namely, if the detector for the considered decays has a certain length L, the produced (on-shell) massive neutrino N could survive during its flight through the detector, and would decay later outside it. Such decays are thus not detected and should be eliminated from the width and the branching ratio of the considered process M ± → ± 1 ± 2 M ∓ , by introducing a suppression factor (nonsurvival probability) where t ≈ L/β N is the time of flight of N through the detector (β N is the velocity of N in the lab frame), and γ N = (1 − β 2 N ) −1/2 is the Lorentz time dilation factor. Hence, the suppression factor, which should multiply the branching ratio, is In the last relation, we used β N ≈ 1 and τ N = 1/Γ N [≡ 1/Γ(N → all)], in the units used here (c = 1 = ). This decay-within-the-detector probability P N has been discussed and presented for the processes with intermediate onshell particle (such as N ) in Refs. [16,[37][38][39][80][81][82]. In this respect, here we follow mostly the notations of Ref. [39]. Usually, the quantity P N is small and is then written as which agrees with Equation (22) in the limit of small P N . The suppression factor (22) can be rewritten as Here, L N is the decay length, and L N is the canonical decay length (independent of the mixing parameters B N ) where K and Γ(M N ) are from Equations (16)- (18), cf. also Figure 2. Equation (24b) suggests that it is convenient to define a canonical (i.e., independent of mixing) probability P N for the decay of N within the detector as We present the inverse canonical decay length, L −1 N , for γ N = 2, in Figure 3 as a function of M N . We note that L −1 N increases very fast (as M 5 N ) when M N increases. Therefore, the supression due to the factor P N may not necessarily be strong (i.e., P N 1) for semileptonic LNV decays of heavier mesons M ± , such as B ± . If we use eyeball estimates for the coefficients N N of the left-hand Figure 2, approximate expressions for the factor K of Equation (18) for Majorana neutrinos can be written In order to estimate better the values of K Equations (27) and thus the suppression factor P N Equation (24), we need to know the present upper bounds for the squares |B N | 2 as a function of M 2 N . These upper bounds we take from compilation of values of Ref. [29], based in turn on upper bound values obtained in Refs. [83][84][85][86][87][88][89][90][91][92][93][94][95][96]. We present them in Table I, for specific chosen values of M N in the mass range of interest. We remark that the upper bounds have in some cases strong dependence on the precise values of M N , see Ref. [29] for further details. In order to use only rough estimates for the values of K, we present in Table II Tables I and II, the upper bounds for |B τ N | 2 are at present significantly less stringent and are expected to become more stringent in the future. When we combine Equations (24b) with (27) and Table II, we obtain for the decay-within-the-detector probability P N ≡ P N K the following estimates and upper bounds, relevant for the K decays (M N ≈ 0.25 GeV), D and D s decays (M N ≈ 1 GeV), and B and B c decays       (32), with L = 1 m and γN = 2, for some of the LNV decays M ± → ± ± π ∓ . The value of MN is chosen such that the maximal value of Br eff is obtained (the value of MN is given in parentheses, in GeV). For M ± = K ± , two different values are given, for = e and = µ. For all other cases, = µ is taken (when = e the values are similar). In order to have the analysis and the formulas simpler, in the rest of this Section we will assume that one mixing parameter, |B N |, dominates over the other two mixing parameters: For example, it may well be that = µ, i.e., that |B µN | |B eN |, |B τ N |. Then, of the branching ratios Br(M ± → ± 1 ± 2 M ∓ ) the largest will be M ± → ± ± M ∓ which, according to Equations (29) and (19) (note that DD * and CC * give the same contribution since 1 = 2 now), is: Multiplying this expression by the probability P N of the decay in the detector, Equation (26), we obtain the effective branching ratio Br eff We see that in the effective branching ratio, Br eff , the complicated dependence on mixing parameters encoded in K [cf. Equation (18)] disappeared because factors K cancel here. All the mixing effects in Br eff are in the simple factor |B N | 4 . Unfortunately, this factor represents a strong suppression, in comparison with Br of Equations (19) where In the last identity (31b) we introduced the canonical (i.e., without any mixing dependence) effective branching ratio Br eff where we recall that Br was defined in Equations (19) and (20). Here Br eff is half the value of Br eff in Ref. (39) because the latter expression referred to the exchange of two Majorana neutrinos instead of one.
Only when M ± = B ± or B ± c , i.e., when the mass of the on-shell N can be high (M N 1 GeV), would it be possible to have P N ∼ 1 [Equation (28c)]; and in such a case Equations (31) do not apply, but rather Equations (30), i.e., Br eff = Br in this case. c , respectively. For the meson decay constants and CKM matrix elements, needed for the evaluation of K 2 factor of Equation (9), and for the masses and lifetimes of the mesons, we used the values of Ref. [97]. The values of the decay constants f B and f Bc were taken from Ref. [98]: f B = 0.196 GeV, f Bc = 0.322 GeV. We note that the presented formulas for Br eff and Br can be evaluated also for the decays M ± → ± 1 ± 2 M ∓ when 1 = 2 . Furthermore, when the final leptons are τ leptons (and M ± = B ± or B ± c ), the values of branching ratios turn out to be similar to those in Figures 6 and 7, but the range of M N in this case is shorter: Table III displays values of Br eff for representative values of M N in the decays M ± → ± ± M ∓ .
As an illustrative example, let us consider the decays D ± s → µ ± µ ± π ∓ , which is one of the preferred decay modes proposed at CERN-SPS [80,81], and, in addition, let us assume that |B µN | 2 is the dominant mixing. In such a case, Equations (31) and Table III imply for the experimentally measurable branching fraction Br eff The present rough upper bound on the mixing for M N ≈ 1 GeV is |B µN | 2 10 −7 , cf. Table II. Equation (33) then implies that Br (eff) 10 −12 for such decays. The proposed experiment at CERN-SPS [80,81] could produce D and D s mesons in numbers by several orders of magnitude higher than 10 12 , which would open the possibility to explore whether there is a production of sterile Majorana neutrinos N in the mass range M N ∼ 1 GeV.
If the decays B ± c → µ ± µ ± π ∓ are considered (we do not use B ± decays as they are CKM-suppressed compared to B ± c ), the results of Figure 7a and Equation (31b) imply an effective branching ratio which is similar to the case of D s , Equation (33) in a detector of the same length (L = 1 m, in our example). Since D s is significantly lighter than B c , the relevant neutrinos that give a sizeable effect are also lighter and thus longer living, implying a smaller P N factor for the same detector length. On the other hand, D s decays have no CKM suppression compared to B c (|V cs | ≈ 1 while |V cb | ≈ 0.04) and the B c channels are more numerous, so that the true branching ratios of the latter are smaller. These two effects compensate in a given detector, as shown in Equations (33) and (34). However for a longer detector the observable D s branching ratio increases considerably [80,81].

III. CHARGED PION DECAYS MEDIATED BY ON-SHELL MASSIVE NEUTRINOS
In the previous Section we considered semileptonic decays of mesons heavier than the pion. For pion decays, there are purely leptonic modes only, since the pion is the lightest meson. In this Section we present and discuss the branching ratios for the LNC decay π ± → e ± N → e ± e ± µ ∓ ν e and the LNV decay π ± → e ± N → e ± e ± µ ∓ν µ . We also  (19) and (20), for these decays.
FIG. 7: The same as in Figure 6, but for the LNV decays of the charmed mesons B ± c .
consider the differential branching ratios dBr/dE µ , where E µ is the energy of the produced µ ∓ . In contrast to the previous Section where the intermediate N neutrino has to be Majorana, here N can be either Majorana or Dirac.
In the Dirac case, only the LNC mode is possible, while both LNC and LNV modes occur in the Majorana case. However, the experimental distinction between these two modes cannot be resolved by simply examining the final state, since the produced neutrino (ν e orν µ ) is not detectable. The major part of this Section refers to Ref. [37], and for certain details we use here results of Refs. [38,39]. The formalism is somewhat more complicated now because we have four final particles (in the previous Section there were three). Nonetheless, several features turn out to be similar as in the previous Section. The considered processes are presented in Figures 8 and 9. As in the previous Section, we consider a scenario with at least one heavy sterile neutrino N , which has suppressed heavy-light mixing coefficients B N with the first three neutrino flavors ν ( = e, µ, τ ), cf. Equation (1). The mentioned decay rates may be nonnegligible only if the intermediate neutrino N is on-shell, i.e., and the process is of the s-type, Figures 8 and 9. Specifically, this means 106 The lepton number violating (LNV) process π + → e + N → e + e + µ − νµ which can be mediated by a neutrino N only if N is a Majorana particle: (a) the direct (D) channel; (b) the crossed (C) channel.
A. Branching Ratios for π ± → e ± e ± µ ∓ ν The decay widths Γ (X) (π ± → e ± e ± µ ∓ ν) (X=LNC, LNV) can be written in terms of the corresponding reduced decay amplitudes T Here, 1/2! represents the symmetry factor from two final state electrons, and d 4 is the integration element of the phase space of the four final particles Here, p 1 and p 2 are the momenta of e; in the direct channel, the e momentum at the first (left-hand) vertex is p 1 , and in the crossed channel it is p 2 , cf. Figures 8 and 9. When we use the expressions for the amplitudes T (X) π,± of Appendix 4, for the specific considered case N 1 = N (and no N 2 ), the decay width (36) can be written as where X = LNC, LNV; k (X) ± are the corresponding mixing factors and Γ (X) π (Y Z * ) are the normalized (i.e., without explicit mixing dependence) decay widths The expressions for T π (CC * ) in Equation (38) have no subscripts ±. Further, in Equation (40), P (X) (Y ) represent the N propagator functions of the direct and crossed channels (Y = D, C) where Γ N ≡ Γ(N → all), and K 2 π is the following constant: It turns out that in the case when the intermediate N neutrino is on-shell, i.e., when its mass is in the interval of Equation (35), the squares of the propagators (41) reduce to simple delta functions due to the inequality Γ N M N where p k = p 1 , p 2 for X = D, C. In this on-shell case, the DD * and CC * contributions in Equation (38) are large and equal, and the interference contributions DC * and CD * are negligible in comparison; we refer to [38] for details on this point. Hence the decay width (38) can be written in the on-shell case as Here, the normalized decay width Γ turns out to be the same for X = LNC and X = LNV, where the following notations are used: and the function F(x µ , x e ) is given in Appendix 5. In the approximation M e = 0, the above result becomes simpler where the function f ( It can be checked that the decay rate (44), with N on shell, coincides with the factorized expression In Appendix 5 we also provide the differential decay rates dΓ (X) (π ± → e ± e ± µ ∓ ν)/dE µ with respect to the final muon energy in the rest frame of N neutrino. They turn out to have quite different forms for X = LNC and X = LNV cases (we will return to this later in this Section).
In order to obtain the branching fractions of the considered processes, we need to divide the decay width by the total decay width of the charged pion, Γ(π ± → all) where the expression (50b) represents the by far most dominant decay mode π ± → µ ± ν µ , and represents a (very small) relative correction coming from the π ± → e ± ν e decay (δg π ≈ 1.3 × 10 −4 ). As we can see in Equations (44) and (45), the decay width Γ(π ± → e ± e ± µ ∓ ν) and its normalized counterpart Γ(π ± → e ± e ± µ ∓ ν) are inversely proportional to the (very small) decay width Γ(N → all) ≡ Γ N , which in turn is proportional to the mixings K ∼ |B N | 2 ( = e, µ, τ ), cf. Equations (16)-(18) and Figure 2. As in Section II A, this effect represents the N -on-shell effect Equation (43) and it makes the considered width Γ(π ± → e ± e ± µ ∓ ν) by many orders of magnitude larger than it would be in the case of off-shell N . In the (narrow) mass interval (35) for on-shell N in the considered pion decays, the factor K in the width Γ N Equation (16) has the following approximate form (cf. Appendix 3 and Figure 2): This expression, within the precision given here, is valid equally for the Majorana and for the Dirac N and agrees with that given in Ref. [38] for the Majorana case. However, the affirmation in Ref. [38] that K for Dirac N is smaller by a factor of two is not correct.
In certain analogy with Section II A, we can define here the canonical branching ratio Br π as the part of the branching ratio Br (X) π ≡ Γ(π ± → e ± e ± µ ∓ ν)/Γ(π → all) which contains no explicit or implicit heavy-light mixing dependence (and is independent of X = LNC or LNV) where the notations (46) are used, and in Equation (53c) we used the expression (125) in Appendix 5. We note that however, the on-shellness (43) brings in factor 1/Γ N ∝ 1/ K ∼ 1/|B N | 2 , reducing the suppression to Br (X) π ∝ |B N | 2 . In Figure 10 we present the canonical branching fraction Br π as a function of M N in the entire interval of onshellness (35). In Figure 11 this canonical branching fraction is presented at the lower edge of the on-shellness interval, where the effects of M e = 0 turn out to be appreciable. We can see that the branching ratio is the largest when M N ≈ 0.130 GeV. The differential branching ratios dBr (X) /dE µ , where E µ is the final muon energy in the N -rest frame, are obtained directly from the differential branching ratios d Γ (X) π /dE µ , and the latter quantity for X = LNV is given explicitly in Appendix 5.
The canonical differential branching ratios, free of any mixing dependence, can be defined in analogy with the definition of the canonical total branching ratio (53a), and, in contrast to canonical branching ratios (53) they do depend on whether X = LNC or Explicit expressions for these quantities, when X = LNV and X = LNC, are given in Appendix 5 in Equations (126) and (128), and in the limit M e = 0 in Equation (129).
The differential (and full) branching ratios for the process π ± → e ± e ± µ ∓ ν differ in the cases when the on-shell N is Majorana or Dirac. When N is Dirac, only the X = LNC process contributes. When N is Majorana, both LNC and LNV processes contribute. We can write the differential and full branching ratios in these cases in terms of their canonical counterparts, by defining first the combined canonical differential branching ratios where 0 ≤ α ≤ 1. Then it is straightforward to check that in the cases of Dirac and Majorana N neutrino, the branching ratios can be expressed in terms of the above quantites (55) where we recall the definition (39) of the coefficients k (X) , and the "Majorana LNV admixture" parameter α M appearing in Equation (56b) is defined as Integration of the relations (56) over E µ leads to the full branching ratios where we used the fact that the E µ -integrated expression Br π (α) is independent of α, cf. Equations (53). In Figures 12 we present these canonical differential branching ratios as a function of the muon energy E µ in the rest frame of N , for four different values of mass M N in the on-shell interval (35), and for five different values of the admixture parameter (α M = 1, 0.8, 0.5, 0.2 and 0), The value α M = 0 corresponds to the case of Dirac N . From the curves of Figure 12 we can conclude: if N is Majorana neutrino with a significant value of the admixture parameter α M and with mass in the interval (35), then the measurement of such differential branching ratios may be able to confirm the Majorana nature of N . If the differential branching ratios are studied with respect to the muon energy E µ in the pion rest frame, the distinction between the Dirac and Majorana case is more difficult, cf. Ref. [37].
B. Effect of the Long Neutrino Lifetime on the Observability of π ± → e ± e ± µ ∓ ν The branching ratios presented in this Section so far, Equations (58) in conjunction with Equations (53), should be multiplied by the probability P N for the on-shell neutrino N to decay within the detector (of length L), as explained in Section II B, Equations (22)-(26) and Figure 3. As a consequence, the effective (true, measurable) branching ratios are not Br (∼ |B N | 2 ) of Equations (58), but Br eff = P N Br In Equations (59a) and (60a) we used the expressions (58) for Br and Equation (26) for P N . In Equations (59b) and (60b) we introduced canonical (i.e., without any mixing dependence) effective branching ratio Br π,eff Br π,eff ≡ P N Br π with Br π given in Equations (53), and the canonical nonsurvival probability P N presented in Figure 3 in Section II B for a wide range of N neutrino masses. In Figure 13 we present P N in the here relevant narrower mass interval (35). As in the case of semileptonic LNV decays of Section II B, we notice that also here in the effective branching ratios, Equations (59) and (60), the complicated dependence on the mixing parameters entailed by the factors K = N N |B N | 2 of Γ N [cf. Equation (18)], cancels out, and there remains only simple dependence on the mixing parameters, in the form |B eN | 2 |B µN | 2 or |B eN | 2 (|B eN | 2 + |B µN | 2 ). Further, in comparison with the branching ratios Br of Equations (58), which are ∼ |B N | 2 , the effective (true) branching ratios Br eff of Equations (59) and (60) are unfortunately significantly more suppressed by the mixing parameters, namely Br eff ∼ |B N | 4 .
The future pion factories, such as the Project X at Fermilab, will be designed to produce charged pions with lab energies E π of a few GeV and luminosities ∼ 10 22 cm −2 s −1 [106,107], and ∼ 10 29 charged pions could be expected per year.
On the other hand, if the larger among the mixing elements (B N | 2 , = e, µ) is |B eN | 2 ( 10 −8 ), the LNV processes dominate, and the effective branching ratios (63) have the following upper bounds: In such a case we have |B µN | 2 < |B eN | 2 10 −8 , and up to 10 7 events could be detected per year. The present upper bounds on |B N | 2 suggest that the first scenario, Equation (64), is more likely. The measurement of the effective branching ratios alone cannot distinguish between the Dirac and the Majorana character of intermediate neutrino N . However, as argued in Section III A and presented in Figure 12, the measurement of the differential branching ratios for the considered processes is a promising way to discern the character of the neutrinos. However, the differential branching ratios (56) must be multiplied by the nonsurvival probability P N , in order to obtain the effective (true, measurable) differential branching ratios dBr eff /dE µ . In analogy with Equations (59) and (60) we obtain, using Equation (56) dBr where the canonical differential effective branching ratios were introduced, in analogy with Equation (55) with dBr π (α)/dE µ given in Equation (55). The parameter α M appearing in Equation (66b) was defined in Equation (57). In analogy with Figure 12, and using the values of P N from Figure 13, we deduce that the values of the y-axes of Figure 12a-d, must be multiplied by P N ≈ 2.0 × 10 −3 , 2.8 × 10 −3 , 3.7 × 10 −3 , and 4.8 × 10 −3 , respectively (assuming L = 1 m and γ N = 2), in order to obtain the representation for the curves of the canonical differential effective ratios dBr π,eff (α)/dE µ .

IV. CP VIOLATION IN CHARGED MESON DECAYS MEDIATED BY MASSIVE STERILE NEUTRINOS
CP violation in the lepton sector could be measured by neutrino oscillations [66]. Here we consider the possibilities of measuring CP violation in meson decays mediated by sterile neutrinos N , such as the semileptonic LNV decays of charged heavy pseudoscalar mesons considered in Section II, or the leptonic (LNV and LNC) decays of charged pions considered in Section III.
It turns out that CP violation in all such decays is possible in scenarios with at least two massive sterile neutrinos N j (j = 1, 2). CP violation in the neutrino sector is expected whether neutrinos are Dirac or Majorana particles. However, in the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) mixing matrix [108][109][110], the number of possible CPviolating phases is larger when the neutrinos are Majorana particles. If n = 3 + N is the total number of neutrinos (N is the number of sterile neutrinos), the number of CP-violating phases is n(n − 1)/2 if neutrinos are Majorana, and (n − 1)(n − 2)/2 if neutrinos are Dirac, cf. Ref. [111]. CP violation in the decays π ± → e ± N j → e ± e ± µ ∓ ν was investigated in Ref. [38], and in the decays M ± → 1 N j → 1 2 M ∓ in Refs. [39,40]. In both cases, it turns out that, even though the rates are extremely small, the CP violation asymmetry in these charged decays may become appreciable, even close to order unity, when the two intermediate neutrinos can go on shell and are almost degenerate in mass M N1 ≈ M N2 . This is to be contrasted with CP violation in charged meson decays due to the standard CKM mechanism in the quarks sector, where the rates are larger but the asymmetries are much smaller, e.g., of order 10 −4 in K ± → 3π decays [112].
A. CP Violation in Semileptonic LNV Decays M ± → ± 1 ± 2 M ∓ As mentioned above, we will consider the scenario with at least two sterile neutrinos, N j (j = 1, 2), both of which can go on shell in the intermediate state, i.e., with masses M Nj satisfying the condition shown in Equation (2). The processes of interest are again those of Figure 1, except that now in both the direct (D) and crossed (C) channels there are two possible neutrinos exchanged: N 1 or N 2 . The relative measure of CP violation for these processes will be the asymmetry: The corresponding LNV decay widths Γ(M ± → ± 1 ± 2 M ∓ ) will be obtained now, in the scenario of two sterile neutrinos (N = 2), in close analogy with the calculation in Section II A which was performed for the case of one neutrino N (N = 1). The relations (3) and (4) and (9) there hold without change now, but the relations (5)-(8) obtain the following slightly more general form when N = 2: Here, k (±) j are the mixing coefficients and Γ ± (XY * ) ij are 2 × 2 matrices, and represent the normalized (i.e., without the explicit mixing dependence) contributions of N i exchange in the X channel and complex-conjugate of the N j exchange in the Y channel (X, Y = C, D) The expressions for T ± (X)T ± (Y ) * (where X, Y = D, C) are the same as in Equation (7) and are given in Appendix 1, and P j (X) (X = D, C) are the propagators of the exchanged neutrinos N j in the direct and crossed channels We will disregard effects due to non-diagonal neutrino widths in their mass basis. For these details we refer to Ref. [40]. The total decay width of N j , Γ Nj , is given by Equations (16)- (18), where each N j has its own mixing parameter K j Equation (18), i.e., K j = N Nj |B Nj | 2 , where the coefficients N Nj as a function of M Nj are given by the left-hand Figure 2 in Section II A.
As we will see below, the CP asymmetry parameter A CP (M ), Equation (68) First let us calculate the quantities appearing in the numerator and the denominator of the CP violation parameter A CP (M ) Equation (68). We introduce the following notations which will be needed below: For example, in the specific case when 1 = 2 = µ, we have θ 21 = 2(φ µ2 − φ µ1 ) = 2(arg(B µN2 ) − arg(B µN1 )). As in Section II A, when both N j are on-shell, it turns out that the interference contributions to the quantities S ∓ (M ) from the direct (D) and crossed (C) channels (DC * and CD * ) are suppressed by several orders of magnitude in comparison with the contributions from the direct (DD * ) and crossed (CC * ) channels, and we will neglect them (we refer to [39] for details on this point). Then it follows from the expression (69) and In the sum (76) M Nj Γ Nj λ 1/2 (1, y Nj , y 1 )λ 1/2 1, y y Nj , y 2 y Nj Q(y Nj ; y 1 , y 2 , y ) where j = 1 or j = 2; Γ(CC * ) jj is obtained from Γ(DD * ) jj by the simple exchange y 1 ↔ y 2 [cf. Equation (13)]. For evaluation of the CP-violating difference S − (M ), Equation (75), the quantity Im Γ(XX * ) 12 (X = D; C) is of central importance. In the integrand of ImΓ(XX * ) 12 appears as a factor the following combination of the propagators of N 1 and N 2 [cf. Equation (71)]: In Equation (79b) we used the narrow width approximation: and in Equation (79b) the parameter η was introduced which parametrizes any deviation from the naive expectation η = 1. We expect η ≈ 1 when ∆M 2 In Appendix 6 we argue that this parameter η, in the case of near degeneracy ∆M N M N1 , is a simple function of only one variable y ≡ ∆M/Γ N , where ∆M N ≡ M N2 − M N1 > 0 and This implies that in the case of two almost degenerate sterile neutrinos (∆M N M N1 ) we have for the factor in Equation (79b) the following identities: We note that the factor 4 in the denominator on the right-hand side of Equation (83) is nontrivial, because a somewhat different result would have been obtained by a more simple and direct consideration of the expression (79a) in the limit Γ Nj M Nj (j = 1, 2). The result (81) [or equivalently, Equation (83)] has been confirmed also by numerical evaluation of Im Γ(XX * ) 12 , Refs. [38,39] (see below). The mechanism (79) [with the identity (83) , i.e., when the two sterile neutrinos are almost degenerate. If ∆M N ∼ Γ N (i.e., y ∼ 1), then the quantity η(y)/y = y/(y 2 + 1) in Equation (83) is very small and CP violation effects disappear.
The mechanism (79) was used in Ref. [38] in the context of CP violation of leptonic decays of charged pions, and in Refs. [39,40] in the here presented context of CP violation of semileptonic LNV decays of heavy pseudoscalars.
We recall that the expression (79) has the same structure with Dirac delta functions as Equation (10) in Section II A; however, the factors in front of these Dirac delta functions are different now. Therefore, integration over the final particles' phase space can be performed now in the same way as in Section II A, i.e., analytically. This leads, in analogy with Equation (12), to the result where we recall the notations Equations (14) and (82), y Nj ≡ M 2 Nj /M 2 M , the function Q is presented in Appendix 2, and we denoted ∆M We note that in Equation (84) we have not yet assumed the near degeneracy of the two sterile neutrinos.
From here on in this Section, we will consider the case of near degeneracy of the two on-shell sterile neutrinos (81) and (82) hold and the quantities Equations (83) and (84) become appreciable and CP violation can thus become significant. Therefore, we have In this case the identity (81) holds, cf. Appendix 6, and the expression (84a) becomes simpler where in the last identity we used the expression (12), the identity (82), and the fact that Γ N2 /Γ N1 = K 2 / K 1 [cf. Equation (18) for N 1 and N 2 ]. The normalized decay matrix elements Γ ± (XY * ) ij , Equation (71), were evaluated in Ref. [39] also numerically, by Monte Carlo integration and using finite small widths Γ Nj in the propagators. The numerical calculations confirmed the presented formulas, among them the expressions (78), (86), and (81)- (83). The form (81) of η(y) was confirmed numerically with a precision better than a few per mille. The numerical evaluations also confirmed that the directcrossed interference terms (DC * and CD * ) are really negligible. We refer to [39] for details.
Further, the mentioned numerical evaluations gave us values of the N 1 -N 2 overlap parameter δ 1 as defined in Equation (77), i.e., the parameter which appears in the expression (76) and represents the N 1 -N 2 overlap effects. It turned out that the numerical values of the parameters δ j (j = 1, 2), as well as of η, are practically independent of: the channel contribution considered (DD * or CC * ), of the type of pseudoscalar mesons (M ± , M ∓ ), and of the light leptons ( 1 , 2 = e, µ) involved in the considered decays. The numerical results show that the parameter δ ≡ (1/2)(δ 1 + δ 2 ) is a function of only one variable, namely y ≡ ∆M N /Γ N (the same is true for η) The numerical values of the parameter δ are given in Table IV as a function of y. It is not clear whether there exists a simple analytic expression for δ as a function of y. Further, in Figure 14 we present the quantities η/y and δ as a function of y. The values for δ in Table IV are practically equal to the values of the corresponding δ parameter in the rare leptonic decays of the charged pions π ± → e ± N → e ± e ± µ ∓ ν, cf. next Section IV A and Ref. [38].
Branching ratios of experimental significance here can be defined by dividing the expressions S ∓ (M ), Equations (73) and (75) and (76), by the corresponding sum of total decay widths [Γ(M + → all) + Γ(M − → all)], which is practically equal to 2Γ(M ± → all) If employing the canonical (independent of mixing) branching ratio Br(y N ; y 1 , y 2 , y ) ≡ Br(DD * ) defined in Equation (20) of Section II A, and its CC * analog, Br(CC * ) ≡ Br(y N ; y 2 , y 1 , y ), then it turns out that the branching ratios Br ± (M ) of Equation (88) can be rewritten in terms of Br, of the heavy-light mixing parameters |B Nj | and K j = N Nj |B Nj | 2 , of function η(y)/y = y/(y 2 + 1) and the overlap function δ(y) tabulated in Table IV Br This leads to an expression for the CP violation parameter A CP (M ) defined in Equation (68), which now involves only the heavy-light mixing parameters |B Nj | and K j [cf. Equation (18)], the function η(y)/y = y/(y 2 + 1), where y ≡ ∆M N /Γ N , and the overlap function δ(y) tabulated in Table IV: In the usually considered case 1 = 2 (≡ ), i.e., when the considered decays are M ± → ± ± M ∓ , the formulas (89) and (90) get simpler, because in such a case Br(CC * ) = Br(DD * ) ≡ Br, and B 2Nj = B 1Nj ≡ B Nj , and κ 2 = κ 1 = κ A CP (M ) = 4 sin θ 21 These formulas become even simpler if the absolute values of the heavy-light mixings of N 1 and N 2 are equal (but not their phases), i.e., when In such a case, we have all κ ≈ 1, and K 2 ≈ K 1 ≡ K, and therefore the expressions (91) reduce to In these expressions, we assumed, in addition, that the N 1 -N 2 overlap terms are small (O(δ)).
In order to obtain the corresponding effective (true) branching ratios, we have to multiply the above branching ratios Br ± by the decay-within-the-detector probability P Nj ≡ P N K j (L/1m), as in Section II B (L is the length of the detector). Again, the complicated mixing dependence entailed in the parameters K j gets cancelled in this multiplication. When adopting the simplifying assumption Equation (92), i.e., the validity of Equations (93), we obtain the following effective branching ratio Br eff and the CP violation effective branching ratio A CP Br eff : In these formulas, we used the canonical effective branching ratio Br eff as defined via Equations (32) and (20) and depicted in Figures 4-7 as a function of M N (≡ M N1 ≈ M N2 ). The values of the Lorentz factors in the lab system are taken to be γ N = 2 for both N 1 and N 2 , keeping in mind that Br eff scales as 1/γ N . We recall that θ 21 = 2(φ 2 − φ 1 ) = 2(arg(B N2 ) − arg(B N1 )). We notice that on the right-hand side of Equation (94a) there is an additional factor two in comparison with Equation (31b) of Section II B; this factor two comes from the fact that we now have contributions of two intermediate neutrinos N 1 and N 2 , and we neglected the contributions from the N 1 -N 2 overlap (O(δ)).
As at the end of Section II B, let us consider now as an illustrative example the decays D ± s → µ ± µ ± π ∓ . In addition, let us assume that |B µN | 2 is the dominant mixing. In such a case, the estimate Equation (33) is still valid, and the CP-violating difference of the effective branching ratios, A CP Br eff , is obtained by comparison of Equations (94a) and (94b) Since for M N ≈ 1 GeV we have at present |B µN | 2 10 −7 , cf. Table II, Equation (95) means that A CP Br (eff) 10 −12 for such decays. As already mentioned at the end of Section II B, the proposed CERN-SPS experiment [80,81] could produce D and D s mesons in numbers by several orders of magnitude higher than 10 12 , and production of the sterile Majorana neutrinos N j could be explored. Further, if there exist two almost degenerate sterile neutrinos of mass M N ∼ 1 GeV (this is so in the νMSM model [60,61,[67][68][69][70][71][72]), such that y ≡ ∆M N /Γ N ∼ 1, then we would have η(y)/y ≡ y/(y 2 + 1) ∼ 1. In such a case the estimate (95) would imply that the CP-violating difference of effective branching ratios, A CP Br (eff) (D s ), is of the same order as the effective branching ratio Br (eff) (D s ) (if the phase difference |θ 21 | 1). Therefore, if experiments can discover the mentioned νMSM-type Majorana neutrinos, they will possibly detect also CP violation effects coming from the Majorana neutrinos.
The case of B ± c → µ ± µ ± π ∓ is similar to the case of D ± s → µ ± µ ± π ∓ described above, cf. Equations (33) and (34) at the end of Section II B. Therefore, Equation (95) is valid also for CP violation in such decays of B ± c . For the relative advantages and disadvantages of D ± s and B ± c decays, we refer to the comments at the end of Section II B.
B. CP Violation in Pion Decays π ± → e ± e ± µ ∓ ν In this Section, we will only briefly outline the calculation of the CP violation asymmetry in the (LNC and LNV) semileptonic decays π ± → e ± e ± µ ∓ ν as described in Section III. We will assume the presence of at least two nearly degenerate sterile neutrinos N j (j = 1, 2) that can go on shell in the intermediate state, as in Section IV A. The present Section is a similar extension of the analysis of the decays π ± → e ± e ± ν of Section III to two sterile neutrinos. The results of the present Section are largely based on Ref. [38]. Only few details will be presented here, for further details we refer to Ref. [38].
Similarly to the previous Section IV A, the quantities relevant for the CP violation will be where X = LNC, LNV. The total branching ratios are Br ± = Br when N j are Majorana neutrinos, and Br ± = Br (LNC) ± when N j are Dirac neutrinos. We adopt the same conventions and the same notations as in the previous Section IV A. In addition, since we have now LNV and LNC processes, we introduce the additional notations As in Section IV A, the requirement that the quantities sin θ (X) (here: X = LNV, LNC) be nonzero, and the requirement of the near degeneracy of the two neutrinos (∆M N M N1 ≡ M N ) in conjunction with the expressions Equations (79)-(83) for Im(P 1 (D)P 2 (D) * ), are needed in order that the CP violation parameters A (X) π,CP = 0 acquire nonnegligible values. Analysis similar to that of the previous Section IV A (but algebraically more complicated) leads then to the results for the quantities defined in Equations (96). More specifically, the results for the Dirac case, Br π,CP , are the following: The expression for the canonical quantity Br π , appearing in Equation (98a), is given in Equation (53) in conjunction with the notation (46) in Section III A The results for the Majorana case, Br π,CP are the following: The function η(y)/y = y/(y 2 + 1) is the same as in Section IV A (with: y ≡ ∆M N /Γ N ). Even more so, numerical evaluations give for the N 1 -N 2 overlap parameter δ(y) the same values as in the semileptonic decays of Section IV A, cf. Table IV and Figure 14 there.
If we assume that |B N2 | ≈ |B N1 | (for = e, µ, τ ), i.e., Equation (92), then we have K 1 ≈ K 2 ≡ K, and the expressions for A π,CP simplify significantly As in the case of semileptonic LNV decays of the previous Section IV A, we see that the CP asymmetry parameter A π,CP can become appreciable and even of order one if the following two conditions are fulfilled simultaneously: (a) at least one of the angles θ (X) (X=LNC,LNV), defined in Equation (97), is appreciable; (b) the quantity y ≡ ∆M N /Γ N is y ∼ 1 (near degeneracy). In such cases, the estimates for the effective (true) branching ratios Br

V. CONCLUSIONS
We have studied lepton number violating (LNV) semileptonic decays of charged pseudoscalar mesons, specifically π ± , K ± , D ± , D ± s , B ± and B ± c , mediated by heavy neutrinos that can go on their mass shell. We first presented the LNV semileptonic decays of charged Kaons and of the heavier mesons D ± , D ± s , B ± and B ± c , in processes of the form M ± → ± 1 ± 2 M ∓ , mediated by on-shell massive neutrinos, where M is the decaying meson and M a correspondingly lighter meson. We estimated the branching ratios as functions of the neutrino masses and mixing parameters, and found the scenarios where upper limits on the mixing parameters can be obtained. We also studied the effect on the observability of these decays due to the long neutrino lifetime, as the secondary decay vertex is likely to fall outside the detector for the range of neutrino masses that are relevant to these processes.
We then presented our corresponding study of charged pion decays, which in this case are purely leptonic since pions are the lightest mesons. Here we can have modes that conserve lepton number (LNC) as well as modes that violate lepton number (LNV), if the intermediate neutrinos are of Majorana type, while only the former modes occur if the intermediate neutrino is of Dirac type. However, these modes are not distinguished by the final state because the latter involves a standard neutrino, which is not experimentally observable. We find that it could be possible to discern the Majorana or Dirac nature of neutrinos if one is able to observe features in the final state distribution.
We finally explored the possibility of observing CP violation in the lepton sector using these meson decays mediated by massive neutrinos on shell. The CP signal in charged meson decays is the usual asymmetry between the decays of opposite charge mesons. We found that leptonic CP violation may show in semileptonic LNV decays of charged Kaons and charged B mesons, as well as in LNC and LNV decays of charged pions, depending on the mass of the intermediate neutrinos. It turns out that such CP violation becomes appreciable and possibly detectable if there are at least two heavy neutrinos almost degenerate in mass that can go on their mass shell. The neutrino mass splitting that gives maximal CP asymmetries is close to the neutrino decay width. This type scenario fits well into the so called neutrino minimal standard model (νMSM), which contains two almost degenerate Majorana neutrinos of mass near 1 GeV and another lighter neutrino of mass of order 10 1 keV, a model that can explain simultaneously neutrino oscillations, the dark matter and the baryon asymmetry of the Universe. and that elements of the D-C interference matrices Γ ± (CD * ) and Γ ± (DC * ) are simply related If the two final leptons are of the same flavor ( 1 = 2 ), one can use the property that the integration d 3 over the final particles is symmetric under exchange of p 1 and p 2 (because M 1 = M 2 ), and we have the following additional symmetries:

A.2. Explicit Expression for the Function Q
The expression in Equations (12) and (78) is arrived at by using in the integration over the phase space of three final particles [Equations (3) and (4)], for the contribution of the N neutrino, the identity The first identity can be used for the DD * contribution (where p N = p M − p 1 ) and the second for the CC * contribution (where p N = p M − p 2 ). When one uses the identity (10) in the DD * contribution, and the analogous identity for the CC * contribution, the integration over dp 2 N becomes trivial, and the d 2 -type of integrations can be performed. Notice that this is equivalent to the factorization approach Γ(M → 1 N )Br(N → 2 M ), which holds when N is on-shell. The obtained expression for Γ(DD * ) is then the expression Equation (12) when N = 1 [Equation (78) when N ≥ 2] with the notations (14), where the obtained function Q has the following form: In the limit of massless charged leptons (y 1 = y 2 = 0), this reduces to A.3. Calculation of the Total Decay Width of Neutrino N In this Appendix, for completeness, we summarize the formulas needed for evaluation of the total decay width of a massive sterile neutrino N , cf. Equations (16)- (18) and Figure 2.
In Ref. [29] (Appendix C there), the formulas for the leptonic decay and semimesonic decay widths of a sterile neutrino N have been obtained, for the masses M N 1 GeV. For higher values of the masses M N , the calculation of the semileptonic decay widths becomes difficult because not all the resonances are known. Hence, for such masses the authors of Refs. [30,113] proposed an inclusive approach, based on duality, for the calculation of the total contribution of the semileptonic decay width of N . It consists of representing the various (pseudoscalar and vector) meson channels by quark-antiquark channels. This approach was applied for M N ≥ M η ≈ 0.958 GeV. Here we summarize the formulas given in Ref. [30] for the decay width channels (cf. also: [29]). In some of these formulas, twice the decay width is given [2Γ(N → . . .)], signalling the fact that for each possible decay of Majorana neutrino in charged particles, there is an equally possible decay into charge conjugate channel (something not possible if N is Dirac particle).
In Equation (110a) factor 2 was included because for Majorana neutrino N both decays N → − + ν and N → + − ν contribute ( = ). When M N < M η ≈ 0.968 GeV, the following semimesonic decays contribute, which involve presudoscalar (P ) and vector (V ) mesons: where factor 2 in the charged meson channels appears because both decays N → − M + and N → + M − contribute (M = P, V ) if N is Majorana. The factors V P and V V are the CKM matrix elements involving the valence quarks of the mesons; and f P and f V are the corresponding decay constants, whose values are given, e.g., in Table 1 in Ref. [30]. The pseudoscalar mesons which contribute here are: P ± = π ± , K ± ; P 0 = π 0 , K 0 ,K 0 , η. The vector mesons which contribute are: V ± = ρ ± , K * ± ; V 0 = ρ 0 , ω, K * 0 ,K * 0 . If M N ≥ M η (=0.9578 GeV), due to duality the (many) semimesonic decay modes are represented by the following quark-antiquark decay modes [30]: R I2(0, xq, xq)+ (g (q) In the formulas (110) We note that in the evaluation of the total decay width Γ N , the expressions (112a) and (112b) should be added when N is Majorana; if N is Dirac, the same summation should be taken, but the expressions (112a) should be multiplied by 1/2. The same approach is valid also in the case of summation of expressions (110) and (111).
In Equations (110b) and (112b) there appear the following SM neutral current couplings: Further, the neutral current couplings κ V of the neutral vector mesons in Equation (111d) are The following kinematical expressions I 1 , I 2 , F P and F V were used: with λ function defined in Equation (46a). These formulas allow us to obtain the total decay width Γ(N → all) as a function of M N . Using these formulas, we evaluated the coefficients N N , appearing in Equation (18) at the mixing terms |B N | 2 , and presented them in Figure 2 as a function of M N for the cases of Majorana and Dirac neutrino N . One may notice a small kink in the curves of Figure 2 at M N = M η (=0.9578 GeV). This kink appears because at M N ≥ M η the use of duality is made (the replacement of the semileptonic decay channel contributions by the quark-antiquark channel contributions). We can see that the duality works quite well at M N ≥ M η , with the possible exception for the case = τ because τ lepton has a large mass.
A.4. Explicit Amplitudes for N j -Mediated Decays π ± → e ± e ± µ ∓ ν In this Appendix we summarize, for completeness, formulas needed in Sections III and IV B. These formulas were derived and presented in Ref. [38], for the case of exchange of two different neutrinos N 1 and N 2 . Here we summarize them in a slightly more general form, for the case of N different neutrinos N j (j = 1, . . . , N ). In Section III the simpler case N = 1 is taken, as it is sufficiently representative for the branching ratios considered there. In Section IV the case N = 2 (or: N ≥ 2) is considered, with two (on-shell) neutrinos N 1 and N 2 almost degenerate, as in such a case significant CP violation effects can occur in the neutrino sector.
The squared amplitude |T (X) π,± | 2 for the N j -mediated leptonic decays of neutrinos, appearing, for example, in Equation (36) (where X = LNC, LNV), is a combination of contributions from the two channels D (direct) and C (crossed) (cf. Figures 8 and 9), and, in general, of the contributions of N neutrinos N j The constant K 2 π is given in Equation (42), and the mixing factors k In the case of N = 1, these coefficients are in Equation (39). P π,± (CD * )], appearing in Equations (116), get simplified when summed over the helicities of all the final leptons. In the case of the X = LNV processes (cf. Figure 9) they acquire the following A.5. Explicit Expression for Γ (X) π and d Γ (X) π /dE µ for π ± → e ± e ± µ ∓ ν with On-Shell N Equation (45) refers to the expression obtained by performing the integration (36) over the phase space of the four final particles [cf. Equation (37)], of the integrand written explicitly in Appendix 4. In the integration, the on-shellness (43) is assumed, which makes the integration over p 2 N trivial. At the final stage of integration, the differential decay width d Γ (X) /dE µ over the muon energy E µ , in the rest frame of the N neutrino, is performed. The expressions for d Γ (X) /dE µ were written in Refs. [37,38], and we write them down here for completeness. In the case of X = LNV it is The integration of this expression over E µ can be performed explicitly (in Ref. [37] it was performed only in the limit M e = 0). The result is Equation (45) with notations (46), where the function F(xµ, x e ) was obtained in Ref. [38]. We write it down here again, for completeness.
We can obtain the LNV canonical differential decay width, according to Equation (54), from the normalized differential decay width Equation (123). For this, it turns out to be convenient to use the following identity: 2 K Γ(π ± → all) which is obtained by using Equations (42), (16) and (17) and (50). This then gives us where The LNC canonical differential decay width turns out to be It turns out that, upon integration of this expression over E µ , we obtain the same result as in the X = LNV case, i.e., Equations (53c) with (46) and (124), or equivalently, Equations (45) with (46) and (124). We must add that we found a typographical error in Equation (A.16) of Ref. [37], where E 2 must be replaced by 2E 2 , and in Equations (B.1c) and (B3) of Ref. [38], where (E µ ) max should read (M 2 Nj + M 2 µ − M 2 e )/(2M N ). In the limit M e = 0 (which is a good approximation), the canonical differential decay widths (126) and (128) get simplified The full (integrated) canonical branching ratio in the M e = 0 limit is obtained by taking the x e = 0 limit of Equation (53c) where the function f is written in Equation (48).

A.6. Delta Function Approximation for the Imaginary Part of the Propagator Product
In this Appendix we investigate the expression for the imaginary part of the propagator product, Im(P 1 (D)P 2 (D) * ), Equation (79a) of Section IV A. For convenience we introduce in this Appendix the following simplified notations x, M 2 , ∆ and ξ: We note that ∆ > 0 by convention; and 0 < ξ < 2. Further, Γ N1 + Γ N2 = 2Γ N , in accordance with the definition of Γ N Equation (82). Since we always have Γ Nj M Nj (the neutrinos N j are sterile), the relation (80) holds, i.e., We can write the right-hand side of Equation (79a) for Im(P 1 (D)P 2 (D) * ) as Im (P 1 (D)P 2 (D) * ) = R 1 + R 2 where R 1 and R 2 can be written, in our notation, as In Equations (134b) and (134d), the identity (132) was used, and we introduced two (dimensionless) parameters η j (j = 1, 2). We want to obtain these two parameters η j . They can be obtained by integrating analytically the explicit expressions (134a) and (134c) for R j (x) over x. For example, integration of R 1 (x) gives where Therefore, in the case of near degeneracy (∆ M 2 ) we have M 2 * = M 2 . If we now use in the integration over dx the expression (134b) instead, take into account M 2 * = M 2 in the case of near degeneracy, and compare with (135), we obtain the following expression for the parameter η 1 by comparison with (135): where in Equation (137b) we use the usual notation in this paper y ≡ (M N2 − M N1 )/Γ N = ∆/(2M Γ N ). Here we note that ∆ ≡ (M 2 N2 − M 2 N1 ) = (M N2 − M N1 )2M N1 in the case of near degeneracy ∆ M 2 ≡ M 2 N1 . Doing the same procedure with the quantity R 2 , we obtain for η 2 the very same result as for η 1