Lepton universality in a model with three generations of sterile Majorana neutrinos

The extension of the Standard Model lepton sector by three right-handed Majorana neutrinos (heavy neutral leptons, HNL) with masses up to GeV scale is considered. While the lightest HNL is the dark matter particle with mass of the order of 5 keV, the remaining two HNLs ensure standard (active) neutrino mass generation by means of the see-saw type I mechanism. Two heavy sterile neutrinos with quasi-degenerate masses up to 5 GeV can induce the deviation of lepton universality violation parameter in the decays of $\pi^\pm$ and $K^\pm$ mesons from the the Standard Model value. Contours are obtained for the permissible values of this parameter within the framework of two mixing scenarios, taking into account the lifetime boundary for heavy neutral lepton from Big Bang nucleosynthesis in the Universe. When calculating the HNL decay width in the framework of the model with six Majorana neutrinos, three active and three heavy, both two-particle and three-particle lepton decays, essential for masses below the mass of the pion, were taken into account. When calculating the decay widths, the limiting case known as the"Dirac limit"is not used. The results based on the explicit form of mixing matrices for three HNL generations and the diagram technique for Majorana neutrinos, which explicitly take into account the interference terms for diagrams with identical mass states, can lead to some differences in lifetime from the results using the"Dirac limit"and the displacement of the corresponding experimental exclusion contours of the"mass-mixing"type. For the second mixing scenario, a mass region of $460~\text{MeV}

1 Introduction.The main features of the model.
Extension of the Standard Model (SM) leptonic sector by heavy neutral leptons of right helicity (HNL, also referred to as sterile Majorana neutrinos) known for a long time [1,2] has been analyzed multilaterally recently due to the attractive general features of such an extension, within the framework of which the symmetry between right and left neutrinos is restored, new large energy scales are not necessarily introduced, neutrino oscillations and their masses generation by means of the see-saw mechanism are successfully explained and a number of important cosmological applications of the model, such as baryon asymmetry of the Universe generation, description of the inflationary stage of the early Universe and its accelerated expansion at the present time are successfully interpreted [3,4].
New useful features are realized in the construction of the so-called minimal neutrino standard model νMSM [5,6] which is a minimal extension of the SM.In νMSM framework the HNL masses do not exceed the electroweak scale and there are no other new particles up to the Planck scale.Cosmological observations impose significant limitations on the model parameter space, which lead to at least three HNL and establish a strict upper limit on the mass of the lightest active neutrino min(ν e,µ,τ ).The lightest heavy lepton N 1 with mass of the order of 10 keV, the lifetime more than τ Universe ∼ 10 17 sec and mixing parameter of the order of 10 −13 -10 −7 , plays a role of the dark matter (DM) particle in such an extension [7,8].The direct method of N 1 DM detection is due to the possibility of observation of the one-loop decay process of N 1 → γν in galactic media [9].Two remaining heavy leptons ensure the mechanism of mass generation of standard (or active) neutrinos, their masses can vary in a wide range of values up to multiGeV scale, however, enough baryon asymmetry through oscillation-induced leptogenesis can be generated even if they are of the order of MeV and their mass splitting is rather small [10].
In this paper we estimate the possible departure of the lepton universality parameter from zero value due to HNL contributions in the νMSM-like model where an explicit form of mixing for the three lepton generations is used.The Lagrangian of extension has the form where l L = (ν L , e L ) T is the left lepton doublet, ν R are HNL flavor states, (ν R ) c ≡ Cν T R (C = iγ 2 γ 0 ), ν R ≡ ν † R γ 0 , H is the Higgs doublet ( H = iτ 2 H † ), F is the Yukawa matrix and M M is a Majorana mass matrix.After spontaneous symmetry breaking M D = F H = F v (v = 174 GeV) is the matrix of Yukawa term.The full 6×6 mass matrix where the flavor states (ν L ) α , (ν R ) I and the mass states ν k , N I (α = e, µ, τ , k, I = 1, 2, 3) are connected by the unitary transformation where P L is the left projector and U ν , U N are unitary 3 × 3 matrices.Block-diagonal form of the mass matrix (3) looks as where m = diag(m 1 , m 2 , m 3 ), M = diag(M 1 , M 2 , M 3 ), W † MW = diag(m ν , M N ), M is defined by Eq.( 3).In the following diagonalization procedure [17] the unitary W-matrix is represented as an exponent of an antihermitian matrix and decomposed The flavor states are related to the mass states in the following form For the left neutrino the main contribution in the flavor basis is given by the first term in (8) which corresponds to the the well-known phenomenological relation ν Lα = α (U PMNS ) αj P L ν j [16].Deviation from unitarity for the PMNS matrix in the approximations W ∼ O(θ 2 ) (and also W ∼ O(θ 3 )) is given by U PMNS ≃ 1 − 1 2 θθ † U ν and defined by the η = − 1 2 θθ † -matrix.The lagrangian terms for HNL currents interaction with W ± , Z bosons have the form In the following consideration we are keeping only the first and the second order terms in θ, as it is customary to do in the available literature.The HNL mixing parameter is defined in the approximation W ∼ O(θ 2 ) as Θ ≡ θU * N .The standard set of active neutrino masses is defined in the framework of the O(θ 2 ) scenario as a solution of see-saw type I equation with ambiguous definition of M D by means of the U N mass matrix in the HNL sector [17,18] where Ω is an arbitrary orthogonal matrix, ΩΩ T = I.
In conclusion of this Section, simplifying for greater clarity the model to one generation of neutrino ν, we recall the terminology used in the literature (see [4]) for the limiting cases of the mass term parameters in the Lagrangian, Eq.( 2).Case M M =0 is called the "Dirac limit", since the Weyl spinors ν L and ν R represent the left-and right-chiral components of the Dirac neutrino with the mass term M D νν = M D (ν L ν R + νR ν L ).The lepton number is preserved.Case M M << M D is called the pseudo-Dirac limit, since it is possible to divide the four components of the Dirac neutrino into two Majorana neutrinos with left-chiral components (ν Case M M ≥ M D corresponds to the case of seesaw mechanism (M M >> M D is the "seesaw limit"), since we have two Majorana mass states, one of which has a mass m 1 of the order of M 2 D /M M , the other m 2 of the order of M M , and the mixing parameter of the two flavor states is of the order of M D /M M .To simplify calculations, the Dirac limit is used in the available literature when Feynman rules for processes involving active neutrinos and HNL are an analogue of Standard Model rules.In this paper, the Dirac limit is not used, the corresponding diagram technique for processes involving Majorana fermions was developed in [19] and [20] for calculations within supersymmetric models and can be directly applied to HNL production and decays.
For further analysis of the lepton universality within the framework of two characteristic mixing scenarios (Section 3) compatible with cosmoligical limitations (Section 2), calculations are reproduced for two-particle semileptonic HNL decays (Section 4) and calculations are made for three-particle leptonic HNL decays (Section 5) in the model with all six Majorana neutrinos.They are used for the lifetime restrictions on the mixing parameter space (Section 6).Bounds on the characteristic lepton universality parameter, Eq.( 1), in the decays of K ± and π ± are considered in the framework of characteristic mixing scenarios in Section 7.

The lightest HNL as a candidate for the role of a Dark Matter particle
In the following it is assumed that heavy neutral leptons N 1,2,3 are ordered by mass and N 1 is the lightest one.For a mass M 1 of the order of keV the main decay channel is N 1 → ννν.The decay width corresponding to the four-fermion effective Lagrangian defined by Eq.( 10) has the form where α = e, µ, τ .Details of calculation can be found in the Appendix A. Heavy lepton N 1 must not decay at time scale of order of the age of the Universe, which means τ N 1 ≥ 4 × 10 17 sec.This limitation is significantly strengthened when taking into account the one-loop induced decay N → γ, ν, which can give a distinctive signal with photon energy Although the increase of the width is small Γ N →ννν /Γ N →γν ≡ k rad = 8π 27α EM ≈ 128, the limitation on the lifetime can be increased by the eight orders of magnitude due to specifics of the gamma-astronomical observations, see [21,22], providing τ N 1 > 10 25 seconds.It is convenient to introduce the effective mass parameter allowing to associate the masses of active neutrinos with the mixing matrix Θ.The connection of the effective mass parameter with the phenomenological value of mixing Then the N 1 lifetime in seconds can be expressed as and the gamma-astronomical constraint can be rewritten as where we used an estimate for the lifetime τ X = 10 25 seconds.It is shown by solid blue line in Fig. 1.
A known direct constraint from below on the HNL mass is M 1 > 0.4 keV, since the distribution of HNL as fermionic dark matter in the phase space of the galaxy is limited by the distribution for a degenerate Fermi gas (Tremaine-Gunn bound, see [23]).
The cosmological restriction in for the density of N 1 dark matter in the Universe appears in the scenario where the mixing of active and sterile neutrinos Θ is quite small, Figure 1: Bounds for the effective mass parameter Σ α (m D ) α1 , see (15), summed by the flavour index, as a function of dark matter particle mass.See Eq.( 24) for the fine-tuning of mixing.Blue region is excluded at the level τ N 1 > 10 25 seconds by gamma-astronomical non-observation of N 1 → γν decay.Grey areas are more detailed exclusions from HEAO-1, XMM and Chandra experiments which are recalculated for m D using the contours in [27].Violet area is excluded by Tremaine-Gunn bound [23].There is an overproduction of dark matter by means of Dodelson-Widrow mechanism [24] in the excluded light-red area above the red line.and the sterile neutrino has never been in thermal equilibrium.The dominant mechanism of the formation of sterile neutrinos (Dodelson-Widrow mechanism, see [24] ) arises from the active -sterile neutrino oscillations.The energy fraction of sterile neutrinos in the Universe in the case of non-resonant production [25,26] is given by In particular, the density of N 1 particle expressed using the effective mass parameter defined by Eq.( 15) is It leads to a restriction from above on m D summed by flavours The excluded area where Ω N > Ω DM = 0.12 is shown in Fig. 1 by light-red color.
Combining the obtained constraints (18) and (21) for the effective mass parameter we get 3 On the classification of mixing scenarios In the following we consider three possibilities of Ω matrix parametrization most appropriate to the constraint given by Eq.( 22), • "Fine-tuning" of mixing for normal and inverted hierarchies where Ω 2×2 is a 2 × 2 orthogonal matrix.In this form of mixing the constraint is imposed directly on the mass of the lightest active neutrino m lightest (m 1 for NH, m 3 for IH) Here we used the unitarity condition for PMNS matrix assuming that U ≃ U PMNS up to O(θ 2 ).With such form of mixing matrix there is a "fine-tuning" of mixing that explicitly highlights the non-zero mass of the lightest active neutrino, unlike the scenarios considered in the following, where the small finite numerical value of mass is not so significant.The effective mass parameter (m D ) α1 summed by the flavor index α gives a counterpart of the parameter U 2 α which is used in experimental reconstructions, see Section 6.
• Mixing expressed by the real orthogonal rotation matrix Ω ∈ SO(3, R) with the following parameterization where c j = cos α j and s j = sin α j , α j ∈ R are Euler angles.Note that the condition (22) for matrix (25) restricts only α 1 and α 2 angles.Moreover, in this scenario one can take m lightest = 0.
• Mixing expressed by the complex special orthogonal matrix Ω ∈ SO(3, C) with the same parametrization as given by Eq.( 25) but replacement of α j → ω j = α j + iβ j , β j = 0.The same as in the previous case, here m lightest = 0.
Possible deviations of Ω matrix from the form of "fine-tuning" above were analyzed in [28].However, in the following we focus mainly on the form of "fine-tuning" which is consistent with the cosmological constraints in a wide range of HNL dark matter masses, demonstrating also flexibility of the mixing factor in the HNL decays ( 13) and (14), which can vary due to changes both of M 1 and the lightest neutrino mass.
Discussion of ambiguity of the choice of the type of matrix Ω in the general case can be found in [18], the most significant of them are related to the processes of lepton flavor violation [29] in different sectors of the model.The minimal parametric choice Ω =I corresponds to a special case of "fine tuning" with normal hierarchy, when redundant parameters are not introduced.A similar form for the inverse hierarchy occurs when Ω is anti-diagonal matrix.In these two cases where U αI are the elements of U PMNS .In further consideration this case of mixing, Eq.( 26) or Eq.( 27), is designated as mixing scenario 1.
For parametric scenarios in the seesaw type I models which are more interesting for collider phenomenology it is needed to combine very small active neutrino masses of the order of F 2 v 2 /M M with moderately heavy HNL, providing observable signals within the LHC and next colliders energy reach, and enhance at the same time small mixing factors of the order of m ν /M HNL , providing observable rates at the luminosity frontier.This is achieved either by fine-tuning of the mixing matrices in a specific scenarios with additional symmetries [30], or in the framework of Casas-Ibarra diagonalisation with complex-valued parameters.First sort of models gives quasi-Dirac neutrinos processed by the standard calculation technique, which are not fully consistent with the second sort of models beyond the "Dirac limit", where evaluations are performed with Majorana fermions.In the latter case Ω 2×2 , Eq.( 23), is chosen as an element of SO( 2, C) In further consideration, this choice is designated as mixing scenario 2.
Three new parameters are introduced in (28), ξ = ±1, Re(ω) and Im(ω).The Dirac limit of scenario 2 has been considered in detail in the literature.Significant enhancements of the collider signals appear with complex-valued ω parameter which leads to the factors X ω = e Im (ω) in the mixing matrix Θ. Detailed phenomenological analyses of active and sterile neutrino mixing in [31] showed that a phenomenologically consistent hierarchy of mixings Θ eI , Θ µI and Θ τ I with suppressed Θ eI relative to other matrix elements can be achieved in a wide interval of X ω independently on the values of HNL masses.Translating the experimental upper bounds on Θ αI from the shortest possible lifetimes of N 2,3 from π ± and K ± meson decays into the upper bound on X ω one obtains at the HNL mass scale 10 2 MeV Im(ω) = 4.5 for the lifetime of the order of 1 sec and Im(ω) ∼ 7 for the lifetime of the order of 0.01 sec.Values of Im(ω) > 6 − 7 lead to a large mixing parameters of the W ± , Z-neutrino interactions not consistent with the data.Complexvalued parametrization of Ω was also used for the study of HNL properties at the TeV scale [30,35], see also [36].
Extensive literature is devoted to the study of the question of the number of HNL generations.For νMSM model the case of only two HNL generations in comparison with the case of three generations has been analysed within the cosmological framework in [5] for arbitrary Ω and diagonal M M with the result that the number of HNL generations equal to three is preferred.Neutrino phenomenology for the case of two right-handed neutrinos has been analysed in [32] where the decoupling limit of the three right-handed neutrino model has been constructed using a specific form form of Ω-matrix in the basis where M M matrix and the mass matrix of charged standard leptons are diagonal and real with an underlying symmetry for the Yukawa couplings or the elements of M M (texture zeroes).Constraints on thermal leptogenesis and LFV processes have been found for such a case.In the presence of a sufficiently large number of acceptable cosmological scenarios, we will adhere to the framework of the νMSM model with dark matter production through activesterile neutrino mixing, where the mass difference of N 2 and N 3 is small in comparison with the known mass splittings of the light left-handed neutrino mass states [33].
Significant recent reconsideration for the case of two HNL generations in the region of the parameter space corresponding to the mass less than the mass of K-meson and performed taking into account the available set of modern data, see [34], is discussed in Section 4 below in connection with the comparison for the case of three HNL generations considered in this paper.The analysis of lepton universality within scenario 2 in the Dirac limit assuming the νMSM framework in the approximation of the two-particle decays was performed in [15].
Since in the following the decomposition of the anti-Hermitian matrix W, Eq.( 6), by powers of θ to the second order terms is used for the transition to the mass basis of leptons, the question naturally arises about the scope of applicability of such a decomposition and taking into account the O(θ 3 ) terms of the decomposition (and higher).Within the framework of a non-minimal O(θ 3 ) decomposition, it is necessary to take into account the terms of the order of O(θM D ) when [37] M whereas, within the framework of the standard mininmal approximation for the see-saw mechanism, it is assumed that M N = M M .For non-minimal decomposition of the W matrix, the condition must be met which is a condition for the self-consistency of the diagonalization procedure, taking into account the O(θM D ) terms.For scenario 2, the mixing matrix, in addition to a small parameter of the order of m/M , contains a potentially large factor of the order of Ω −1 , the limited contribution of which must be checked.This issue is discussed in Section 6.

Semileptonic HNL decays
HNL decays with meson in the final state can be divided in four groups [38]: • pseudoscalar neutral meson M 0 ps = π 0 , η, η ′ , K 0 , D 0 , B 0 in the final state with the decay width • pseudoscalar charged meson • vector charged meson M ± v = ρ ± in the final state with the decay width where G F is the Fermi constant, V qq ′ is the Cabibbo-Kobayashi-Maskawa (CKM) matrix element, f M and g M are the corresponding meson decay constants [38], HNL , m α is the mass of charged lepton l α , λ is two-particle kinematic function, κ M is an additional dimensionless correction factor for hadronic matrix element 0|J Z µ |M 0 v (see [38] for details).Values of κ M are given in Table 1.These two-particle widths, see Values for correction factor κ M in Eq.( 33) from [38], s W = sin θ W sine of the Weinberg angle Fig. 2a and Fig. 3, introduce essential contributions to the total N 2 width starting from the π 0 threshold.Black lines -decay with π 0 , red -K 0 , green -η, blue -η ′ , orange -D 0 , purple -B 0 , dashed black -ρ 0 , dashed gray -ϕ.This type of decay is not sensitive to flavor and mixing scenarios 1 and 2 for both hierarchies do not have significant differences.

Leptonic HNL decays
In this section we calculate squared amplitudes and HNL decay widths for the four-fermion effective interaction terms when M 1,2,3 ≪ M W , M Z .Three different decay amplitudes with respect to the mixing factors can be distinguished, ) keeping all masses of leptons nonzero by integrating in invariant variables over the Dalitz plot.The details of these calculations are given in Appendix B. In the model under consideration, all six neutrinos are Majorana fermions, which requires careful determination of the signs of the interference terms.See in this connection [19] for the fermion flow technique for diagrams with Majorana fermions, or [20] for the case of generic basis of γ-matrices.Three cases for amplitudes take place: Case 1.Three active neutrino mass states appear in the final state: N I (p) → ν i (k 1 ), ν j (k 2 ), ν j (k 3 ), i, j = 1, 2, 3. Insofar as the approximation O(θ 2 ) is used, the PMNS matrix can be operated as a unitary one (U † U) = δ ij and there are no Zν j ν k =j vertices.Neglecting the active neutrino masses for the sake of clarity we get the squared amplitude of the following structure Details of the calculations can be found in Appendix A.
Diagram for the case 1. HNL decay into three active Majorana neutrinos.
The mixing factor |(U † Θ) iI | 2 is reduced taking the sum over the active neutrino mass states and the decay width Case 2. Two leptons of different flavors and a neutrino mass state appear in the final state: Neglecting the interference term between diagrams with intermediate W + and W − , Diagrams for the case 2. HNL decay into Dirac charged lepton and antilepton of different flavors associated with one active Majorana neutrino After summing by the active neutrino mass states and using UU † =I due to approximate unitarity of PMNS matrix, the decay width takes the form where In the limiting case m e ≪ m µ , m µ ≪ m τ Eq.( 35) is reduced to where This result is the same as the one obtained in [38].Note that if one assumes that the active neutrinos are Dirac fermions, then there is no interference between diagrams with W + and W − because they correspond to different final states νl + l − and νl + l − .Discussion of the approaches to evaluations for Dirac and Majorana fermions can be found in [40], where the Dirac limit is used for the neutrinos.Beyond the Dirac limit, the interference term between diagrams with Majorana neutrinos vanishes if the active neutrino mass is taken to be zero.To illustrate the suitability of approximations of this kind it is useful to calculate the decay width with interference term using a simplified amplitude neglecting the m M I power terms, and integrating over the triangle Dalitz plot One can observe that the interference term, suppressed by the mass ratio, may be not small in the case of scenario 2 with huge mixing factors of the order of e Im (ω) and for HNL masses at the eV scale (see the discussion of such scales in [36]).The sign of interference term depends on the sign of Re(ω).
Case 3. Two leptons of the same flavor and a neutrino mass state appear in the final state: . The amplitude contains both charged and neutral currents and includes in this case three interfering diagrams Diagrams for the case 3. HNL decay into Dirac charged lepton and antilepton with the same flavor and one active Majorana neutrino.
The mixing factor appearing in the squared amplitude is reduced by summing over the states of active neutrinos Integration of the amplitude gives the decay width where (a) Partial widths of 3-body leptonic HNL decays N → l α , l α , ν evaluated in the Dirac limit (dashed line) and using Eq.(39) (solid line).
(b) Ratios of partial widths for 3-body leptonic HNL decays N → l α , l α , ν evaluated in the Dirac limit and using Eq.( 39) for Im(ω) = 0 (solid line) and Im(ω) = 10 (dased line).In the charged lepton massless limit r = 0, F 1 (r) =1, F 2 (r) =0 and the factor −2C 1 in front of the interference term is positive.Comparison of the decay width calculated using Eq.( 39) and the decay width in the Dirac limit is shown in Fig. 4a,b.The difference in three-particle widths, caused by the opposite signs of the interference terms in Eq.( 39) compared with the Dirac limit, can be several times, however, the main contribution to the total HNL width made by two-particle modes is significantly greater than the threeparticle modes.
Case of HNL in final state.Decay channels with HNL in the final state N 2,3 → N 1 +... are suppressed by the factor of ∼ O(Θ 2 ) in comparison with the channels described above.They are insignificant for the following analysis.
6 Constraints on the mixing parameters for N 2 and N 3 In both scenarios defined in Section 3 the contribution of N 1 , the DM particle, to the lepton universality parameter is small and the degree of lepton universality violation (LUV) depends on N 2 and N 3 .µ for various HNL masses.Green lines are the lower bounds derived from the BBN lifetime limit, left panels -NH, right panels -IH; thin green lines are for the simplified model with two HNL's, thick green lines -for the full model with three HNLs and the lightest active neutrino mass m 1 /M DM or m 3 /M DM of the order of 10 −8 .Excluded domains in gray color correspond to the upper limits from two types of accelerator experiments: experiments with the missing energy reconstruction (TRIUMPH [42], PIENU [43] (π decay), NA62 [44], E949 [45], KEK [46] (K decay)) and experiments with the displaced vertices (DELPHI [47], PS-191 [48], CHARM [49], NuTeV [50] ), the division of the gray region into subdomains corresponding to different accelerator experiments can be found in [34].Seesaw bound contours of U 2 e and U 2 µ in the simplified case of two HNLs and/or the case of three HNLs with m 1(3) = 0 are shown by red and dark red lines.
Taking into account the valuable remark in review [36] regarding the use of the socalled "model-independent approach" in the analysis of data from various experiments3 , we note the need for careful translation when bringing the results to a common denominator.In the general case when the mass and the mixing parameter are not independent variables, the exclusion contours are subject to dependence on the field-theoretic model of the expansion of the lepton sector.The partial probabilities of HNL decays in the model under consideration with six Majorana neutrinos differ from the corresponding probabilities in the model with active Dirac neutrinos and Majorana sterile neutrinos.These deviations can be significant and may lead to some quantitative displacement of the exclusion contours, although a qualitative correspondence will take place.
Significant development beyond the model-independent approach has been performed in [34] where two HNL generations with degenerate masses have been introduced with the 3×2 Ω matrix analogous to Eq. (28).In this case m 1 =0 in the active neutrino mass matrix for NH and m 3 =0 for IH.

Lower bound for BBN
Cosmological considerations imposing restrictions on the lifetime of N 2 and N 3 on the level of τ N 2,3 < 0.1 -1 sec [53] were recently improved in [41] giving the minimal level of 0.02 sec.These bounds are model dependent and obtained in the framework of a rather specific scenarios of the Big Bang nucleosynthesis.A simplified estimate for BBN lifetime limit is used in the following evaluations τ HNL < τ BBN = 0.02 sec, M HNL > 140 MeV 0.1 sec, M HNL < 140 MeV (43) In the scenario 1 framework the minimal mixing matrix (26) does not contain redundant parameters, so the lifetime dependence on the HNL mass is unambiguous, see Fig. 6.
For scenario 2, it is necessary to take into account the constraint for Ω matrix elements following from the self-consistency condition of the model with Casas-Ibarra diagonalization extended to the cubic terms in the decomposition of the W matrix, see Eq. (7).The exponential factor e 2Im(ω) may give a huge increase of the mixing parameters Using the constraint of Eq.( 29) for the orthogonal matrix Ω, one arrives to and for m 2,3 ≃ O(0.1 eV) the inequality must hold In this scenario first the allowed lifetime domain on the Im(ω) − M plane is found, see Fig. 7, which is then translated to the allowed domain for the lepton universality parameter ∆r M , see Fig. 8 and Fig. 10a.In the phenomenological analysis, in the following we use the commonly accepted denomination for the model parameter space where m 1(3) is the mass of lightest active neutrino for normal (inverted) hierarchy.The case of simplified model with two HNL generations, see [34], can be reproduced in the limiting case of the model under consideration when m 1(3) → 0 is taken, which is valid for the range of HNL masses where the lower bound of U 2 α ≫ |Θ α1 | 2 at fixed non-zero m 1 or m 3 .

Lower bound for seesaw
In the limiting case of Im(ω) → 0 which is the transition to the real-valued Ω matrix, it is necessary to take into account all terms, since they become comparable in the order of magnitude lim The lifetimes of N 2 (green lines) and N 3 (dark red lines) as a function of their masses M 2 = M 3 in the case of normal (solid lines) and inverted (dashed lines) hierarchies, scenario 1. Red line correspond to BBN constraints (43), which is taken from [41].N 2 contour for the case of normal active neutrino hierarchy practically coincides with N 3 contours for both hierarchies.lim where φ = Re(ω).If zero mass of the lightest active neutrino is taken, m 1 = 0 (NH) or m 3 = 0 (IH), then the values of mixing matrix elements form the "seesaw bound" as it is called in the existing literature, see for example [34].Note that in the case of nearly degenerate M 2 and M 3 inherent to νMSM-like model which we adhere, there is no φ-dependence of the mixing parameter U 2 α in the case of Im(ω) = 0 when Θ-matrix has the form ( 49) or ( 50) (due to sin 2 φ + cos 2 φ =1).Seesaw bound is inroduced to mark a minimal possible value of phenomenological parameters4 , in particular for U 2 α we can write It is appropriate to call this parameter at m 1 = 0 or m 3 = 0 as the "absolute seesaw bound" for NH or IH, and if we choose a non-zero mass of the lightest active neutrino, then to call (51) as just a "seesaw bound".These two bounds coincide when MeV such inequatity is not respected, since so the mixing component of Dark Matter HNL N 1 becomes the dominant term of U 2 α , or at least the same order of magnitude term as the other terms for M 1 ∼ 1 keV.The difference in bounds is illustrated in Fig. 5.

Restrictions on the lepton universality parameter
in the decays of π ± and K ± mesons Limitations on the lepton universality parameter are imposed by the restrictions on the N 2 and N 3 lifetime from the Big Bang scenario and experimental restrictions from meson decays.
Although the calculations were performed taking into account all the abovementioned decay modes, it is instructive to write out simple formulas for the case of leptonic twoparticle decays.In the effective four-fermion approximation the width of the scalar meson where G F is Fermi constant and f M is meson formfactor, α = e, µ, τ (only channels allowed by energy-momentum conservation are admitted), I = 1, 2, 3 number of HNL generation, j = 1, 2, 3 active neutrino mass states; In the case of only one active neutrino in the final state the mass corrections are neglected and Convenient variable to analyse BSM deviations from lepton universality is the ratio [11] or its derivative demonstrating the deviation of the ratio from zero where masses of active neutrinos are neglected and only decays which are allowed kinematically M I < m M − m α are accounted for.If M I > m M − m α then G M αI = 0.The SM parameter and the BSM parameter are Using the unitarity condition for the full 6 × 6 mixing matrix one can rewrite (56) in the form

Numerical analysis
Taking into account the combined restrictions (see Fig. 5), the allowed values of the mass are M > 430 MeV for U 2 e -bounds and M > 300 MeV for U 2 µ -bounds.Consequently, in the case of π-decay, the lepton universality is violated due to the non-unitarity of the PMNS matrix and For K-meson decay the function G K eJ is nonzero in an allowed range 430 MeV < M < 497 MeV.Moreover, if M 2 and M 3 are quasidegenerate, then G K α2 ≡ G K α3 and For the dark matter fermion N 1 we take M 1 = 5 keV.The values for the contribution of terms with m 1 = 10 −5 eV are given in Table 2. Components of the mixing matrix are can be neglected in comparison with other terms.As demonstrated by Fig. 9, the kinematic factor G K µ2 does not exceed 5, so the value of (G Using this approximation, the parameter of lepton universality violation (LUV) can be written in a simple form In the following the upper and lower bounds are denoted by u and u for U 2 e and by v and v for U 2 µ .Then the maximum and minimum values of the LUV parameter for the pion are and for the kaon (M > 430 MeV) If the mass of the lightest active neutrino m 1(3) = 0 (or, equivalently, if there are only two generations of HNL) then δ dm = 0.In the considered mass range M > 410 MeV, one can assume that u, v > 10 7 , so u ≫ v, δ dm and v ≫ u, δ dm .Therefore, when the HNL mass is greater than the kinematic threshold M K − m e , it follows that (∆r non−unitary ) max ≃ v and (a) Lepton universality violation (LUV) parameter in K + decays with the mixing matrix corresponds to Eq.( 28), scenario 2. Orange lines correspond to the minimum value of LUV in the case of three HNL generations when , where m 1(3) ∼ 10 −5 eV.Solid line denotes the case of NH, dashed line is for IH.Blue line shows a maximum value of LUV parameter for both hierarchies, see Eq.(66) for details.Horizontal dotted black line is the experimental value of ∆r K ≃ 10 −3 , see [56].Vertical green line is a kinematic threshold M = M K − m e for kaon decay K + → e + , N 2(3) .
(b) Lepton universality violation parameter in π + and K + decays when the mass of HNL is greater than all kinematic thresholds for considered decay channels.In such a case LUV occurs due to PMNS non-unitarity: The maximum value of ∆r non−unitary is estimated from the upper experimental bound on U 2 µ , see Eq.(63) for details.
Figure 10: Lepton universality violation (LUV) parameters ∆r π and ∆r K as a function of the HNL masses in the case of quasidegenerate M 2 ≃ M 3 ≡ M. The BBN bound is respected for the depicted mass interval according to the allowed domains in Fig. 5.The choice of M 1 and the active neutrino masses for both normal and inverted hierarchies are the same as in Fig. 8.
(∆r non−unitary ) min ≃ −u.The case 450 MeV < M 493 MeV is separately considered since LUV occurs due to the decay channel K + → e + , N 2,3 with a kinematic factor up to 10 5 (see Fig. 9) and in such case These boundaries are shown in Fig. 10a.Maximum value of LUV parameter for HNL mass greater than the kinematic threshold (M > 493 MeV) is shown in Fig. 10b.

Summary
The most general case of extending the lepton sector of the SM by three sterile Majorana neutrinos in order to generate the masses of standard neutrinos using the see-saw mecha-nism allows one to build a structured hierarchy of mixing parameters within a well-defined basis for mass states.Cosmological limitations on the lifetime and energy fraction of the lightest mass state of neutral heavy leptons, considered as a dark matter particle, restrict the mass to vary in an interval of 0.4-40 keV within the sensitivity of modern experiments, which allows the use of preferred forms of active and heavy neutrino mixing matrices for the analysis of data from experiments with extracted beams and colliders.Limitations of the possible form of the mixing matrix in the see-saw type I models, which are of undoubted interest, can be obtained by combining constraints on the matrix elements imposed from above by the absence of signals on colliders and beam dump experiments, while constraints on the mixing from below are provided by baryogenesis scenarios within the framework of the Big Bang concept, as well as restrictions on flavor oscillations.However, the numerical values for the boundaries strongly depend on the mixing scenario within a particular model which, from a technical point of view, is ambiguously implemented by a certain choice of the Ω-matrix in the M D − U P M N S − U N connection, Eq.( 12).The lepton universality violation parameter in K ± and π ± decays is sensitive to the neutral heavy leptons due to their additional contributions to meson decay widths, dependent on the mixing matrix.
Two mixing scenarios analysed above are rather different, in the "fine-tuned" scenario 1 the mixing matrix depends on the mass ratio (m i /M i ) 1/2 of active neutrino and HNL, which suppresses any HNL production or decay channel and imposes rather strict restrictions on N 1 from the cosmological lifetime and DM energy fraction requirements, not affecting N 2 and N 3 masses.In the scenario 2 with three additional parameters for the N 2 and N 3 mixing, intensively discussed in the literature, a sort of tuning is needed to generate small masses of active neutrinos by means of see-saw type I and the baryon asymmetry of the Universe by means of N 2 − N 3 flavor oscillation mechanism.HNL production is enhanced by the mixing parameter e Im(ω) and stronger restrictions on the N 2 and N 3 masses (N 2 ∼ N 3 ) from below are imposed not affecting to large extent N 1 which has a negligible mixing.Strong hierarchy of Θ eI and Θ µI , Θ τ I elements of the mixing matrix naturally appears in scenario 1 and can be easily configured for the NH case in scenario 2. The cosmological upper bound on the HNL lifetime is critical for determining the bound on the masses of N 2 and N 3 , so their lifetime was carefully evaluated taking into account the three-particle leptonic decay channels which are dominant below the thresholds of the two-particle channels.Our main results can be summarized as • In the range of HNL masses 120 MeV < M < 140 MeV (close to the mass of πmeson) a small window for the U 2 e consistent with the data on the missing energy reconstruction and search for displaced vertices was found in the case of a normal hierarchy (NH): 3 • There is a possible range of parameters consistent with the experimental upper bound for U 2 e , U 2 µ and the limitations of Big Bang Nucleosynthesis (BBN) with the following dependencies on the mass hierarchy: Using the theoretical SM value for R π and R K from [54,55] and the experimental value from Particle Data Group [56], the lepton universality violation parameter for π ± is found to be ∆r π = (−4 ± 3) × 10 −4 and ∆r K = (4 ± 4) × 10 −3 for K ± .A maximal value of the lepton universality violation parameter ∆r K is comparable or exceeds the experimental value in the HNL mass range between 460 -485 MeV.It follows that the currently achieved level of experimental accuracy is quite moderate, as a result of which future experiments of high precision have a great potential for the discovery of BSM physics lepton mixing scenarios.
Some differences between the BBN exclusion contours obtained in the literature and in evaluations above can take place due to approximations for the explicit form of mixing matrices for three HNL generations in the scenarios under consideration, and calculations of the widths of three-particle HNL decays beyond the Dirac limit.The HNL lifetime bound is essential for manipulations with displaced vertices.If a lifetime bound is denoted by τ 0 and a decay width is factorized in the form of a mixing factor times matrix element squared, Γ 0 = |U α | 2 M 2 , then a bound for the mixing factor is |U α | 2 > 1/(τ 0 M 2 ).This qualitative estimate gives a shift of the |U α | 2 contour upwards (downwards) when either the lifetime bound or the matrix element squared decreases (increases).For HNL masses less than the threshold of around 0.14 GeV, see [41], the BBN contour in Fig. 5 can be shifted upwards, which is explained by a step-like approximation of the upper limit for the lifetime τ 0 increase to approximately 0.1 sec.At masses exceeding the threshold of 0.14 GeV the contour in Fig. 5 can be shifted downwards, which is explained by different ways of taking into account the permissible areas for Im ω and M, see Fig. 7.When generating a contour in this paper, the lower permissible boundary on the Im ω − M plane exactly corresponds to a displacement along the curve in Fig. 7, whereas for the case of two HNL generations either asymptotic behavior is used, see Eq.( 44), or the limiting case of the so-called "dominant mixing" is taken, when the ratios have the form of a ratio of some numerical constants, which is equivalent to a stepfunction approximation.Departures in the case of inverse neutrino mass hierarchy are a consequence of completely different structure of mixing matrices for NH and IH in our case.Some displacement of the contours appears also due to contributions to the width by the interference terms or different signs of these terms in the three-particle decays beyond the Dirac limit, however, against the background of two-particle contributions above the thresholds, this displacement is not very significant.As a result of evaluations, in our case the open windows for HNL masses are slightly smaller and for larger values of M HNL , the lower permissible limit is changed.Significant changes occur with the threegeneration seesaw bound curves compared to the two-generation absolute seesaw bound curves due to an additional term in the parameters U 2 e and U 2 µ which includes non-zero active neutrino masses.The abovementioned changes do not affect the mass range, which is most interesting for observing the violation of lepton universality.

Z
Table 4: Feynman rules for Majorana fermions used in calculations of pure leptonic decay amplitudes.The direction of the fermion flow is chosen from the bottom up, as indicated by an arrow near the vertex, P R,L = (1 ± γ 5 )/2, s W = sin θ W . Nonstandard pseudovector structure of ν νZ vertex, for example, follows from the equality of vector current to zero for Majorana fermion.
Diagram for the case 1 with with an explicitly specified choice of fermion flow direction indicated by arrows.
Denoting the terms of the full decay amplitude as one can observe that there is no interference in the approximation of zero masses of active neutrinos in the final state For the amplitude terms we get and finally for the full squared amplitude Identical particles in the final state give an additional factor 1 2 for decay width.Note that for Dirac active neutrino terms LR and RR do not arise and the factor in front of right-hand side will be 32.It is also necessary to take into account the charge conjugate mode by multiplying by two.There is no charge conjugated mode for Majorana neutrinos.In the Dirac limit, see [38], Eq.(3.5) for Γ(N → νν ν), summation over flavor indices and multiplication by a factor of 2 for charge conjugated final states gives the same result as Eq.( 34) above.
Case 2. Acting by the same rules in combination with Fiertz transformations, see Table 2, the amplitudes M the squared terms are The decay width for all lepton masses nonzero is given by Eq. (35).
Case 3.For the three diagrams in this case, the notation M 3 XY , XY = LL, LR, RL, RR is used for the neutral current diagram and the notation M 1 , M 2 is used for W + and W − exchange diagrams, respectively.Then the full amplitude |M WZ | contains 21 terms: 6 terms for the squared diagrams and 36−6 2 = 15 terms for the interferences.Moreover, we need to change the fermion flow depending on each specific interference term to build a trace of gamma matrices.For instance, we need to change the fermion flow for Zνl-vertex in the diagram with W + in order to calculate M 1 M † 3 XY .Full interference term looks as Symbols "′" and "′′" mark those diagrams in which it was necessary to change the fermion flow for calculation, where and we use notation (s W = sin θ W where θ W is the Weinberg angle) The square of the full amplitude is simplified after summing by the neutrino mass states, Eq.( 25), the decay width for lepton nonzero masses is given by Eq. (39).Note that for Dirac neutrinos, opposite signs of interference terms above appear.The equations for the boundaries of the physical region in terms of invariants are obtained requiring G(s 1 , s 2 , 1, r 2 , r 1 , r 3 ) = 0, so for s 1 they are (a) Branching ratios for semileptonic twobody decays of N 2 with charged pseudoscalar or vector mesons and charged lepton in the final state.Solid lines correspond to e − in final state, long-dashed lines -to µ − , dotted lines -to τ − .Blue lines -decay with π + , green -K + , orangeρ + , purple -D + , red -B + .(b)Branching ratios for pure leptonic threebody decays of N 2 .Gray linesννν decay mode, green linesνe − e + , purple lines νµ − µ + , dashed purple linese − µ + , dashed blueνe − τ + , dashed orangeνµ − τ + , red linesντ − τ + .mode.There is a very weak sensitivity of partial widths in relation to the choice of the mixing scenario.

Figure 2 :
Figure 2: Branching ratios for HNL decays.Upper plot -scenario 1 for both the normal and the inverted hierarchy.Middle plot -scenario 2 with Ω = Ω N H (ξ, ω) where ξ = +1, Re(ω) = 0 and Im(ω) = 7.This parameter set comes from the Big Bang nucleosynthesis (BBN) limits, see Fig.5.Bottom plot -the same as the middle plot but for the case of IH.

Figure 3 :
Figure3: Branching ratios for semileptonic two-body decays of N 2 with neutral pseudoscalar or vector mesons and active neutrino in the final state.Curves for NH and IH practically coincide, they are evaluated for the same parametric scenarios as in Fig.1.Black lines -decay with π 0 , red -K 0 , green -η, blue -η ′ , orange -D 0 , purple -B 0 , dashed black -ρ 0 , dashed gray -ϕ.This type of decay is not sensitive to flavor and mixing scenarios 1 and 2 for both hierarchies do not have significant differences.

Figure 4 :
Figure 4: Comparison of the HNL leptonic decay widths calculating in the Dirac limit and in model with 3 Majorana active neutrinos

Figure 5 :
Figure 5: Constraints on U 2e and U 2 µ for various HNL masses.Green lines are the lower bounds derived from the BBN lifetime limit, left panels -NH, right panels -IH; thin green lines are for the simplified model with two HNL's, thick green lines -for the full model with three HNLs and the lightest active neutrino mass m 1 /M DM or m 3 /M DM of the order of 10 −8 .Excluded domains in gray color correspond to the upper limits from two types of accelerator experiments: experiments with the missing energy reconstruction (TRIUMPH[42], PIENU[43] (π decay), NA62[44], E949[45], KEK[46] (K decay)) and experiments with the displaced vertices (DELPHI[47], PS-191[48], CHARM[49], NuTeV[50] ), the division of the gray region into subdomains corresponding to different accelerator experiments can be found in[34].Seesaw bound contours of U 2 e and U 2 µ in the simplified case of two HNLs and/or the case of three HNLs with m 1(3) = 0 are shown by red and dark red lines.

Figure 7 :
Figure 7: The BBN constraints for HNL lifetime (green line for N 2 and dark red line for N 3 ) as a function of its mass in scenario 2 with following parameterization of mixing Ω = Ω NH (ξ = +1, ω) for NH (a) and Ω = Ω IH (ξ = +1, ω) for IH (b) with Re(ω) = 0. Blue lines correspond to the restriction from above on the Im(ω) parameter (47) with ǫ = 1 (solid line) and ǫ = 10 −3 (dotted line).The allowed parameter domains for corresponding HNL are above the lines marked by N 2 (green line) and N 3 (dark red line) and below the blue solid line.

Figure 8 :
Figure 8: Lepton universality parameters ∆r π , red lines, and ∆r K , blue lines, as a function of the HNL mass in the case of quasidegenerate M 2 ≃ M 3 ≡ M for scenario 1.First generation HNL N 1 with the mass M 1 ≃ 5 keV is the dark matter particle.Vertical green lines are the lower bounds from BBN lifetime constraints for IH (left line) and NH (right line).Horizontal dashed lines corrspond to experimental values of LUV parameters: ∆r π = 5 • 10 −4 and ∆r K = 10 −3 , see [56].

M > 430 ••
MeV for U 2 e with NH, M > 430 MeV for U 2 e with IH; M > 290 MeV for U 2 µ with NH, M > 300 MeV for U 2 µ with IH.Combining this constraints, we can conclude that M > 430 MeV in addition to the allowed window of HNL mass mentioned above.In the model with three generations of sterile Majorana neutrinos where the lightest active neutrino mass is m 1 ∼ 10 −5 eV and M DM ∼ 1 keV the lower limit for HNL mass coming from the U 2 e BBN bound M > 430 MeV (NH) is raised in comparison with the simplified model with two generations of HNL where m 1 = 0 and M > 170 MeV.The lepton universality violation parameter ∆r M , M = π ± , K ± does not exceed O(10 −4 ) at M 2 and M 3 masses of the order of 10 2 MeV in the scenario 1, demonstrating very high sensitivity to the BBN lifetime restrictions (M > 480 MeV for IH and M > 830 MeV for NH).In allowed BBN region LUV parameter does not exceed O(10 −10 ) for NH and O(10 −7 ) for IH.The situation is more involved in the threeparametric scenario 2 where the lifetime bound is given by an exclusion contours on the mixing -mass, Im(ω) − M plane.

←
1 and M 2 for the diagram with intermediate W + and the diagram with intermediate W − can be written as ν k (k) Diagrams A2.Diagrams for case 2 with an explicitly specified choice of fermion flow direction (thin arrows).