Numerical estimate of minimal active-sterile neutrino mixing for sterile neutrinos at GeV scale

Seesaw mechanism constrains from below mixing between active and sterile neutrinos for fixed sterile neutrino masses. Signal events associated with sterile neutrino decays inside a detector at fixed target experiment are suppressed by the mixing angle to the power of four. Therefore sensitivity of experiments such as SHiP and DUNE should take into account minimal possible values of the mixing angles. We extend the previous study of this subject arXiv:1312.2887 to a more general case of non-zero CP-violating phases in the neutrino sector. Namely, we provide numerical estimate of minimal value of mixing angles between active neutrinos and two sterile neutrinos with the third sterile neutrino playing no noticeable role in the mixing. Thus we obtain a sensitivity needed to fully explore the seesaw type I mechanism for sterile neutrinos with masses below 2 GeV, and one undetectable sterile neutrino that is relevant for the fixed-target experiments. Remarkably, we observe a strong dependence of this result on the lightest active neutrino mass and the neutrino mass hierarchy, not only on the values of CP-violating phases themselves. All these effects sum up to push the limit of experimental confirmation of sterile-active neutrino mixing by several orders of magnitude below the results of arXiv:1312.2887 from $10^{-10}$ - $10^{-11}$ down to $10^{-12}$ and even to $10^{-20}$ in parts of parameter space; nonzero CP-violating phases are responsible for that.


Introduction
Neutrino oscillations clearly call for an extension of the Standard Model (SM) of particle physics. Sterile neutrino models can provide a simple theoretical framea e-mail: iv.krasnov@physics.msu.ru b e-mail: timagr615@gmail.com work explaining this phenomenon, which makes them popular among many candidates for physics beyond SM. In this framework one commonly introduces three Majorana fermions N I , I = 1, 2, 3, sterile with respect to SM gauge interactions SU (3) c × SU (2) W × U (1) Y .
One can write down the most general renormalizable sterile neutrino Lagrangian as: where M I are the Majorana masses, and Y αI stand for the Yukawa couplings with lepton doublets L α , α = e, µ, τ and SM Higgs doublet (H a = ǫ ab H * b ). When the Higgs field gains vacuum expectation value v = 246 GeV, the Yukawa couplings in (1) yield mixing between sterile N I and active ν α neutrino states. Diagonalization of the neutral fermion mass matrix provides active neutrinos with masses m i and mixing which are responsible for neutrino oscillation phenomena. If active-sterile mixing angles are small, then active neutrino masses are double suppressed, m ∼ U 2 M . This is the standard seesaw type I mechanism, for more details one can address [2].
The authors of Ref. [9] suggest that in upcoming particle physics experiments sterile neutrinos of GeV scale may appear in heavy hadron decays and can be detected as they decay into light SM particles. In proposed fixed target experiments such as SHiP [10] or 2 DUNE [11] the main source of sterile neutrinos are Dmesons decays, so mostly only sterile neutrinos with masses below 2 GeV are produced. We note that in a part of parameter space such sterile neutrinos can be responsible for leptogenesis in the early Universe which allows for direct laboratory tests of early time cosmology [12,13,14].
The minimal values of mixing between sterile and active neutrinos, consistent with the type I seesaw mechanism, have been estimated in [1] for the case when some of sterile neutrinos are lighter than 2 GeV and CPviolating phases are set to zero. Our estimates (and that of Ref. [1] as well) are performed for mixing with only two sterile neutrinos: the third sterile neutrino is considered to be unobservable in the discussed experiments. One scenario is that the third sterile neutrino can be too heavy to be produced in the mentioned experiments. It can also be too light to be kinematically recognizable there. Or it can be of interesting mass range, but very feebly interacting, i.e. practically decoupled. Although we don't specifically restrict this mixing in such a way, it would include also the special case where this sterile neutrino can play a role of warm dark matter 1 [15], from cosmological constraints its mass is restricted to be in keV range. Another widely studied special case is when the two heavier sterile neutrinos are strongly degenerate in mass: it is shown [16] that such sterile neutrinos can be used to provide explanation for leptogenesis (this model is called neutrino minimal extension of the SM or just νMSM). There are other sterile neutrino models that explain leptogenesis, and, consequently, baryon asymmetry of the Universe, for example, through Higgs doublet decay [14]. In this paper we don't introduce any special restrictions, but for numerical analysis we choose M 1 = 500 MeV, M 2 = 600 MeV, and then study how the results change with masses.
For the sterile neutrino masses about 500 MeV current upper limit on active-sterile mixing angles come from CHARM experiment [17] at |U e | 2 , |U µ | 2 ∼ 10 −6 . It can be improved in the near future by SHiP experiment [10], which plan is to reach |U e | 2 , |U µ | 2 ∼ 10 −9 . In more detail present bounds and expected sensitivities of some future experiments one can found, for example, in Ref. [10].
The lower limit on mixing is usually associated with seesaw mechanism, Big Bang Nucleosynthesis (BBN) and some limits that are specific to concrete model. The seesaw limit is a plain mathematical limitations 1 We call that case "warm" dark matter as opposed to "cold" dark matter (as in WIMPs case) due to the fact that keV scale sterile neutrinos have significantly non-zero velocities at the equality epoch. Cosmic structure formation gives an upper limit on those velocities, but at the level of 10 −4 .
imbued on model (1) to make it consistent with current central values for active neutrino parameters. This limit is the primary object of this study. Independent limit comes from cosmological role of sterile neutrinos. If at least one sterile neutrino decays during BBN, products of decay would change light element abundances. Observation of these abundances gives us limit on contribution of non-standard BBN scenarios.
The goal of this paper is to study the dependence of minimal mixing, consistent with the seesaw mechanism, on CP-violating phases, unaccounted before in Ref. [1]. During this research it also became obvious that the role of lightest active neutrino mass can be essential. The obtained results can be used to estimate the sensitivity of future experiments required to fully explore the parameter space of type I seesaw models with sterile neutrinos in the interesting mass range.

Active sector
Before going into details of concrete sterile neutrino model we provide the parametrization of active neutrino sector used in this paper.
For active neutrino masses, we use diagonal matrix m ν ≡ diag{m 1 , m 2 , m 3 }. Experiments provide us with two related parameters and 3σ allowed ranges (see p. 248 of [18] = ∆m 2 . It is usually defined that m 2 > m 1 (just for convenience), but we don't know which of m 1 , m 3 is smaller. Hence, the sign of ∆m 2 is unknown, and two different cases have to be considered: the normal hierarchy case m 1 < m 2 < m 3 and the inverted hierarchy case m 3 < m 1 < m 2 . Hereafter the values (values in brackets) correspond to normal (inverted) hierarchy of the active neutrino masses. If there is no difference, we don't use brackets at all. Absolute value of the lightest mass m lightest differs greatly from model to model and is not specified by present experiments. So in this paper it is treated as one of the free parameters.
For the normal hierarchy we have: 3 and for the inverted hierarchy: Cosmology constrains sum of the masses from above as [19]: We choose m lightest < 0.2 eV as the approximate constraint equivalent to (3). The transformation from flavour basis to massive basis is provided by hermitian conjugate of Pontecorvo-Maki-Nakagawa-Sakata unitary mixing matrix U P MN S : which in turn, can be parametrized as follows (see p. 248 of [18]): Here c ij and s ij stand for cos θ ij and sin θ ij , with i, j = 1, 2, 3, i < j. The angles entering (5) have been experimentally determined. We take best fit values and 3σ allowed ranges for sin 2 θ ij (see p. 248 of [18]): As it is more convenient to use angles θ ij themselves, rather than the values of sin 2 θ ij , we list them here (corresponding to the best fit value): CP-violating phases δ, α 1 , α 2 entering (5) are still not specified by experiments as strictly as angles θ ij and are one of the main subjects of study in this paper. In the most general case we have δ ∈ [0, 2π), α 1 ∈ [−π, π), α 2 ∈ [−π, π). For the Dirac phase δ we have the best fit value (see p. 248 of [18]): At 3σ no physical values of δ are disfavoured (see p. 248 of [18]). Basically we treat δ as a free parameter, but always provide graph for best fit value (7) if possible. Majorana phases α 1 , α 2 haven't yet been observed in any experiment. Consequently, they are considered as free parameters, for more detail see Sec. 5.

Sterile sector
Next we move on to the subject of sterile neutrino sector parametrization.
It is convenient to adopt the bottom-up parametrization for the 3 × 3 Yukawa coupling matrix Y [20]: where Sterile neutrino mass scale is not fixed in the seesaw mechanism, and in this paper we use for the numerical simulations values from the mass range M I < 2 GeV, most relevant for the upcoming fixed target experiments. We discuss different cases of sterile neutrino mass spectrum in Sec. 3. Matrix R is a complex orthogonal matrix, R T R = 1 3 . We use the following parametrization of R: where c i = cos z i , s i = sin z i , and z i ⊂ C are not restricted in any way.

Matrix of mixing angles
The main subject of this work is the matrix of mixing angles between active and sterile neutrinos: It depends on three complex angles z i entering (9) and three yet unknown CP-violating phases of U P MN S matrix (5). Also we know to a certain extent three U P MN S matrix angles (6) and two differences in active neutrino masses squared (2), but neither mass of the lightest neutrino m lightest nor the hierarchy of masses. We consider sterile neutrino production in a fixedtarget experiment due to mixing (10) in weak decays of hadrons. As D-meson decays are the main source of sterile neutrinos in the mentioned fixed target experiments, sterile neutrinos with M I > 2 GeV are too heavy to be produced. Sterile neutrino main signature is a weak decay into SM particles due to the same mixing. The number of signal events depends on the values of |U Iα | 2 , I = 1, 2, 3; α = e, µ, τ . Obviously, the sign matrix in (9) can be omitted during calculation of |U Iα | 2 . We should note that due to kinematics, mixing |U Iτ | 2 doesn't play any role in the decays 2 of sterile neutrino emerged in decays of charmed hadrons. Hence we study the minimal values of |U Ie | 2 and |U Iµ | 2 in order to determine what maximal sensitivity the coming experiments should achieve to fully explore the type I seesaw model.
The case of one sterile neutrino being lighter than 2 GeV (e.g. M 1 > 2 GeV, M 2 > 2 GeV, M 3 < 2 GeV) is of little interest as sterile neutrinos N 1 , N 2 can't be observed in the discussed experiments in this case. During the scan of possible values of unrestricted parameters z 1 , z 2 , z 3 , α 1 , α 2 , there always can be found such a set of these parameters, that |U 3e | 2 = 0, |U 3µ | 2 = 0. Because of that we can't rule out this model even if a fixed-target experiment doesn't detect sterile neutrinos with ultimately high precision. In this paper we only consider the case when one of sterile neutrinos doesn't have significant contribution to the mixing with active sector. In case when M 1 , M 2 , M 3 2 GeV it has to be specifically implied, as all three sterile neutrinos kinematically can be produced in D-meson decays. This noninteracting sterile neutrino can serve as a dark matter candidate, as it is decoupled from others [15]. If DM sterile neutrino are produced by oscillation in primordial plasma, its mass scale is in keV-range [15]. From experimental point of view such sterile neutrinos can't be kinematically recognizable in fixed-target experiments. If one wants to stay in the boundaries of νMSM to simultaneously provide dark matter candidate and leptogenesis in the early Universe, one needs two heavy sterile neutrinos masses to be degenerate. It was suggested, though, that leptogenesis could be successful in a much wider range of masses of sterile neutrinos, given that they are at the same scale, including GeV region, and mix to the active neutrino with comparable strength [13]. It is estimated in [6] for M I < 5 GeV case that mixing U 2 µI 10 −10 is con- 2 Here we count only observable decay modes; |U Iτ | 2 governs decay into unrecognisable τ -neutrino sistent with the leptogenesis scenario. The lower limits on mixing for M 1 , M 2 , M 3 < 2 GeV considered in [6] are |U µI | 2 ∼ 10 −13 for m lightest = 0.23 eV and |U µI | 2 ∼ 10 −12 for m lightest = 0. Study of the case of all three masses being below 2 GeV scale, M 1 , M 2 , M 3 < 2 GeV and none of the sterile neutrinos being decoupled from the active fermions (and so potentially discoverable in a beam-dump experiments), can be a subject for further research.
Lastly, the third neutrino can be heavy M 1 < 2 GeV, M 2 < 2 GeV, M 3 > 2 GeV. Naturally, in this case N 3 kinematically can not be produced in D-meson decays and has no effect in these experiments.
So for our setup the relevant observables are mixing angels between N 1 , N 2 and ν e , ν µ . Our aim is to find the lowest sensitivity enough to rule out the seesaw mechanism. It implies the absence of any signal of either of sterile neutrinos. Hence the relevant combinations to constrain are: Thus we search for minimal values of U e , U µ , which, at a given M 1 , M 2 guarantee full exploration of the seesaw mechanism for such case. It can be seen that U e and U µ don't depend on M 3 in this particular case. In our numerical studies, unless stated otherwise, we set M 1 = 500 MeV, M 2 = 600 MeV. We discuss what happens for other spectra in Sec. 5.3. We point out that from eq. (10) one can see that our results can be rescaled to other mass scales. If one simultaneously changes sterile neutrino masses by factor X: M I → XM I , than mixing also simply changes by that factor: U e → 1 X U e , U µ → 1 X U µ . We choose mass scale that can be tested in proposed fixed target experiments [10,11], but our results can be simply scaled for the case of heavier sterile neutrinos. Note that for heavier sterile neutrinos the main source of production is not the meson decays, but the decays of heavier SM particles, e.g. weak gauge bosons produced in colliders such as LHC, FCC.

BBN constraint
We should note, that for sterile neutrinos at GeV scale we have constraints from the Big Bang Nucleosynthesis. They follow from the fact that sterile neutrino decay products would change light element abundances originating from BBN. Sterile neutrinos can be born in the early Universe, although we don't consider any specific mechanism in this paper. They are not stable due to mixings with active neutrino, and may decay during BBN. SM products of sterile neutrino decays are very energetic and can destroy atoms that has already been produced, thus changing chemical composition of the Universe. Direct observations imply limits on how much new physics can affect these abundances.
These limitations are mainly independent from the seesaw constraint, and can significantly change with the introduction of some new physics affecting active-sterile neutrino mixing in the early Universe. One such example is the inflation theories which introduce coupling of the inflaton to sterile neutrinos, such as Ref. [21].
One can consider two realistic scenarios with small mixing. The first is that mixing is significant enough for sterile neutrinos to come to equilibrium and depart from it in the early Universe before BBN, and then, by the time of BBN, decay in SM particles. In this case we can obtain lower limit on the mixing. Second scenario is when mixing is greatly suppressed and sterile neutrinos never equilibrate. In such a case if mixing with active sector is lowered even further, BBN can no longer restrict mixing in this area.
First of all, the production rate of sterile neutrinos can be expressed as [22,23]: where T is plasma temperature, G F is the Fermi constant, c(T ) ∼ 1 is a numerical parameter varying slightly with temperature and sin 2 2θ is mixing parameter for the case where only one sterile neutrino mixes with only one active neutrino. We neglect actual numerical coefficients, differing for mixing with different active neutrinos at different temperatures, because that is of little importance for our estimate. From this point on we use |U | 2 instead of sin 2 2θ which corresponds to our case, then in total there are three sterile neutrinos. Equilibrium is achieved at H = Γ , where H is the Hubble parameter, with M P l being the Planck mass and g * standing for the effective number of degrees of freedom in plasma.
One can obtain numerically that for M I = 500 MeV the equilibrium can be achieved for |U | 2 |U b | 2 ≈ 2 × 10 −11 . The boundary value |U b | 2 corresponds to the situation when sterile neutrinos come into thermal equilibrium and exit it immediately at temperature T eq ∼ 10 GeV; with smaller mixing the sterile neutrinos would never be in equilibrium. One can get values of |U b | 2 , T eq by equating (12) and (13), expressing |U | 2 as function of T and finding its minima and T corresponding to it. We take c(T ) = 0.76 in accordance with [23], g * (10 GeV) = 86 1 4 . Physically T eq is the temperature of maximal production in (12).
After decoupling, the sterile neutrino concentration is: The smallest mixing we obtain for the seesaw model (see Figs. 5 -5) are typically smaller than |U b | 2 . The relevant case is then the second scenario: the long living sterile neutrino that never was in equilibrium. The BBN constraints one can obtain from Ref. [24]. Naturally, if sterile neutrino is never abundant enough for equilibrium, it's concentration is less than that in equation (14), and since n I ∝ Γ ∝ |U | 2 it can be estimated as: This rough estimate is enough for our purposes.
In [24] the BBN constraints for decays of new particle X are introduced for the variable: where s = g * 4π 2 90 T 3 is entropy density. In our case (16) matches to: where g * (10 GeV) = 86 1 4 , g * (1 MeV) = 10 3 4 and ǫ is the part of energy coming to concrete decay channel.
Combining (19) with Fig. 10 in Ref. [24] one can see that BBN excludes only mixing |U | 2 above a certain, although significantly small, value. We show in this paper that such straightforward limitation can be a stronger constraint than seesaw mechanism constraint. BBN constraint is independent from seesaw mechanism constraint and can be avoided with introduction of the new physics. On side note, BBN doesn't constrain too weak mixings, which we show might also be allowed by seesaw constraint. As we discuss in Sec. 5.2.3, this scenario can't be tested in any experiments in the near future, so we don't study it in details. We state that our numerical estimate can't recognise values of mixing below |U | 2 ∼ 10 −20 , that lays in the zone still excluded by BBN (19). This qualitative estimate is correct for the interesting sterile neutrino mass range 100 MeV < M I < 2 GeV.
Note in passing, as we explained at the end of Sec. 3, if we rescale the sterile neutrino mass range from 100 MeV < M I < 2 GeV to 1 GeV < M I < 20 GeV, we should change the values of |U | 2 in (19) by the factor of 0.1. From the definition, |U b | 2 also changes by the factor of 0.1, while T eq changes by the factor of (10) 1 3 ≈ 2.15. That corresponds to ζ I changing by the factor of 10 in (19). On the other hand, sterile neutrino lifetime changes by the factor of 10 −5 as seen from (18). To be more sincere, the lifetime takes even lesser values as at higher masses of the sterile neutrinos new decay modes become available. From Fig. 10 in Ref. [24] one can see that such changes should lessen the constraints coming from BBN for sterile neutrino with greater mass.
After obtaining values of parameters which correspond to the minimum, we write them in a file, and restore the value of U e , U µ using eq. (10), (11). Such approach guarantees that our results correspond to the actually possible vales of U e and U µ . Experimentally the mixing U e , U µ can be constrained from heavy meson decays as we discussed in Sec. 3, and we obtain the minimal mixing, which is still consistent with the seesaw mechanism.
Concrete formulas used for calculation are laid out in Appendix Appendix A. Note that possible values δ ∈ [0, 2π) can be restricted to δ ∈ [0, π), as (11) can't distinguish between some points of this area, as also shown in Appendix Appendix A.
To make it easier for applications, in Appendix Appendix B we provide a semianalytical approximation for curves presented in Figs. 3 -9. For calculated parameters of these fits see Appendix Appendix B.

Zero phases
Firstly we study the dependence of minimal U µ on U e with zero CP-violating phases for a set of values of m lightest and two hierarchies of masses. This work was done in [1]. To make comparison possible, firstly we perform calculations using parameters, presented in [1]. Obtained graphs should have been identical to those of Ref. [1], but differences between them have been significant enough to require an explanation.
Through discussion with authors of [1] we came to a conclusion, that declared values of experimentally observed variables didn't match ones used in calculation (the adopted values were provided by authors of [1]): Our graphs, constructed with these parameters, perfectly match ones from Ref. [1]. In Figs. 5 -3 we present the results calculated with presently accepted central values of neutrino parameters (2), (6) for both hierarchies and we use only these values from now on as well.
We plot dependence of minimal U µ on minimal U e for the normal hierarchy and δ = α 1 = α 2 = 0 and different values of m lightest on Fig. 5. Fig. 2 is a zoom in a small scale area of Fig. 5, studied in [1]. On Fig. 3 we plot the same dependence for the inverted hierarchy.
To make further descriptions more tangible, we take for each specific curve the value of U µ when U e = 0 and the value of U e when U µ = 0 and call them characteristic values of U µ and U e for this curve correspondingly. From the form of the graphs one can see that they represent the maximal values of U µ and U e . Usually we present only the smallest characteristic values. Thus for               the normal hierarchy defined in that way the characteristic values of the seesaw mixing are U µ ≈ 1.5 × 10 −11 , U e ≈ 2 × 10 −12 . For the inverted hierarchy they are U µ ≈ 3.5 × 10 −11 , U e ≈ 7.5 × 10 −11 . Basically, curves have the behaviour of "the greater the mass the higher the curves lay". One can see that, for zero CP-violating phases, "m lightest = 0" curve corresponds to the lower limit of the values of mixing. To fully explore type I seesaw model with corresponding sterile neutrino masses for zero CP-violating phases one just has to reach the sensitivity corresponding to that lower limit.

Non-zero phases
In this section we study the dependence of minimal U µ on U e , but with non-zero phases. We only lay out here some of the more characteristic graphs.

Inverted hierarchy
We plot dependence of minimal U µ on minimal U e for the inverted hierarchy and different values of m lightest on Fig. 4. The difference between graphs calculated using different values of δ can't be observed with naked eye. As dependence on δ doesn't play much role for these graphs, we only include graph for minimization on δ, α 1 , α 2 .
For the inverted hierarchy and minimization on δ, α 1 , α 2 ( Fig. 4) a significant difference can be seen as compared with zero phases case (Fig. 3). First of all while the curve corresponding to m lightest = 0 practically doesn't change its position, other curves change their behaviour of "the greater the mass the higher the curves lay" to the opposite one of "the greater the mass the lower the curves lay". In this way, in case of nonzero CP-violating phases the lower limit corresponds to the curve with the highest possible mass. From graphs in Fig. 4 one can see that with growth of m lightest characteristic values of U µ , U e lose as much as several orders of magnitude. Here we define the characteristic values in the same way as we did it in 5.1. m lightest = 0 curve keeps these values at U µ ≈ 3.5×10 −11 , U e ≈ 7.5×10 −11 , not differing much from δ = α 1 = α 2 = 0 case. For "m lightest = 0.05 eV" curve characteristic values are estimated to be U µ ≈ 1.4 × 10 −12 , U e ≈ 2.7 × 10 −11 . Characteristic values diminish even more rapidly with further growth of m lightest , and for "m lightest = 0.1 eV" curve these values are indistinguishable from the point of origin, U µ ∼ U e ∼ 10 −20 . This behaviour exceeds our expectation and we study dependence on the m lightest in more detail in Sec. 5.2.3. Nevertheless, this result shows that the estimation of the upper limit of m lightest can become the leading factor in determining the theoretical lower limit on the mixing angles. Non-zero values of minimized α 1 , α 2 are responsible for the difference with the results of Sec. 5.1.

Normal hierarchy
For the normal hierarchy the difference with zero phases case takes more complex shape.
For zero phases (Figs. 5, 2) we have the "the greater the mass the higher the curves lay" behaviour. It changes completely, as curves start to cross each other, behave differently in areas with big values of m lightest compared to the areas of small values. For some masses the mixing angle takes such minuscule values, that the corresponding curves become indistinguishable from the point of origin (they lose several orders of magnitude up to U µ ∼ U e ∼ 10 −20 ). Moreover, for different values of δ the graphs differ significantly from each other.
We study more closely the dependence of graphs on m lightest and δ to understand such behaviour in the next Sec. 5.2.3.

Dependence on m lightest
As each curve in Figs. 5 -9 monotonously declines with growth of U e , for convenience we choose the case of U e = 0, representing the highest value of U µ , as the characteristic point for each curve and study the dependence of U µ on m lightest .
Firstly, we plot dependence of minimal U µ on m lightest for the normal mass hierarchy and U e set to zero for a set of δ-phases (Figs. 5, 11). Starting from the m lightest = 0 point on the graph, where values of all curves lay at the same magnitude of 10 −11 , the curves start to decline as m lightest grows. The only exception is δ = π curve, which exhibits a short growth before it reaches a local maximum and starts to decline like every other curve. At this point all curves on this plot show a common behaviour: they decline swiftly, their values lose several orders of magnitude over minuscule increase in m lightest . For linear scale it seems as if values of U µ swiftly decline to zero values 3 and stay this way for a wide range of values of m lightest , before they start to increase just as swiftly as they have declined earlier. Due to its distinguished form we call this problematic area the "plateau". Another point of interest is that the curves no longer depend on delta after m lightest becomes great enough to exit "plateau" area, uniting into one curve, as can be seen on Fig. 5.
We study the "plateau" proximity area more closely to understand what values our function can actually reach in the "plateau". We plot dependence of minimal U µ on m lightest at U e = 0 in the "plateau" proximity for the normal hierarchy in logarithmic scale for δ = π 2 on Fig. 5. We obtain all active and sterile neutrino parameters listed in 2.1, 2.2 and by putting their numerical values in (10), we obtain our resulting values for U e , U µ in accordance with formulas in Appendix Appendix A.2. Even if we set U e = 0 analytically, numerically it can be preserved only to a certain degree. After reaching the aforementioned values of U µ ∼ 10 −20 for U µ with U e = 0 our numerical calculation has reached the limit of it's accuracy. It is visible that the closer we get to "plateau", the closer is the reconstructed value of "zero" U e to the value of U µ itself. Studying this area of parameters any deeper requires more refined procedure and is beyond the scope of this paper given the fact that mixing |U Iα | 2 ∼ 10 −20 can't be tested by experiments in the foreseeable future.
The reason why a small increment in the values of m lightest brings such drastic drop in the values of U e , U µ turns out to be a mutual subtraction. As we minimize α 1 , α 2 in plateau area they can be chosen in such a way that U µ1 , U µ2 can loose their leading orders, their absolute values dropping drastically. On a side note, the mixing with third sterile neutrino in such case can be more intense than mixing with other two, but it won't be observable in considered experiments if M 3 > 2 GeV.
Although we say that we have reached the limit of our calculation's accuracy, it doesn't mean that any of results presented here are inaccurate. What we mean is that the values are no higher than the ones we provide, but in the "plateau" area they can be even lower than U µ ∼ 10 −20 .
We plot dependence of minimal U µ on m lightest for the inverted mass hierarchy and U e set to zero for δ = 0 on Fig. 13. For inverted hierarchy there is no significant dependence on δ for all values of m lightest and so we don't include graphs with other values of δ. Starting from the m lightest = 0 point on the graph, where values of the curve lay at the magnitude of 10 −11 , the curves start to decline as m lightest grows. If one looks in linear scale, curve simply reaches "zero" value and after that stays in the "plateau". At m lightest = 0.1 eV curve still doesn't leave "plateau".

Dependence on sterile neutrinos masses
In this Section we study the dependence of the minimal mixing on the values of sterile neutrino masses. First of all we notice that explicit dependence of mixing values U αi can be taken from (10): M , U µ ∼ 1 M (as already stated in [1]). In general case M 1 = M 2 and mixing may be distributed between sterile neutrinos in uneven manner. In fact our numerical estimation shows that, for the minimal values of mixing we are interested in, the lighter sterile neutrino is practically decoupled. We plot (M 1 , M 2 ) space for the normal hierarchy on Fig. 5. Likewise, we plot (M 1 , M 2 ) space for the inverted hierarchy on Fig. 15. The regions below the lines will be ruled out by experiments with sensitivity U c = 5×10 −11 for the normal hierarchy. The different lines correspond to different values of m lightest . Parameter U c represents the sensitivity of the experiment needed to rule out the seesaw model in a specified sterile neutrino mass region. From these Figs. one can see that the lines corresponding to the minimal possible values of M 1 , M 2 follow the equation M 1 + M 2 = f (m lightest )/U c . This behaviour is the same as one found in [1], although concrete dependence is modified with the introduction of parameters δ, α 1 , α 2 .
For the normal hierarchy one can see that lines start to lay lower with the growth of m lightest until they reach zero in "plateau" area, and start to lay higher with the growth of m lightest after m lightest exits plateau values. For the inverted hierarchy lines lay lower with the growth of m lightest until they reach the minimum value at m lightest = 0.065 eV and start to lay higher with the growth of m lightest after that.
We plot dependence of M 1 + M 2 on m lightest for U c = 5 × 10 −11 for the normal hierarchy on Fig. 16 and inverted hierarchy on Fig. 17. The mentioned above behaviour can be seen on this graphs in more detail.
We plot dependence of parameter U c on the heavier sterile neutrino mass for the normal hierarchy needed to rule out seesaw model for fixed M 1 = 400 MeV, m lightest = 0.005 eV. Dependence of parameter U c on the heavier sterile neutrino mass for the inverted hierarchy for fixed M 1 = 400 MeV, m lightest = 0.005 eV. This way M 2 is the heavier mass in region 400 MeV < M 2 < 1 GeV. One can see that the minimal value of U c monotonically decreases with the growth of the heavier sterile neutrino mass. Minimal mixing values practically don't depend on the lighter sterile neutrino mass.

Conclusion
In this paper we study minimal possible mixing angles between sterile and active neutrinos |U Iα | 2 for the specific case of two sterile neutrinos with masses less than 2 GeV. These angles provide us with information on sensitivity which experiments such as SHiP or DUNE or their successors should achieve to fully explore type I seesaw model with two sterile neutrinos with masses below 2 GeV and one undetectable sterile neutrino. To that end we study the dependence of mixing matrix on model parameters (δ, α 1 , α 2 ), that hasn't been considered in work [1]. Characteristic values for zero phases are |U Iα | 2 ∼ 10 −11 . Introducing the dependence on CPviolating phases, we observe strong dependence on the lightest neutrino mass m lightest and these phases. For both hierarchies minimal mixing |U Iα | 2 could be lowered depending on m lightest and (δ, α 1 , α 2 ) to the values of 10 −20 at least. These results can be rescaled to other values of sterile neutrinos masses: if we simultaneously change M I → XM I (for all three sterile neutrinos), than mixing also simply changes by that factor: U e → 1 X U e , U µ → 1 X U µ . Such sterile neutrinos can be produced in the decays of weak gauge bosons and other heavy SM particles, e.g. in LHC, FCC. We conclude that still unknown parameters of active neutrino m lightest , δ, α 1 , α 2 may significantly change the mixing pattern and should be taken into account in future experiments.
Thus we can express z 2 , z 3 in terms of z 1 , α 1 , α 2 . We can transform (A.2) using (A.21, A.22) into a more simple form: Directly solving (A.21) one can express z 2 as: Assuming that we have already determined z 2 , the expression for z 3 can be obtained in the same way by solving (A.22): We use these equations to reconstruct U e , U µ by means of definition (A.1), (A.2) and compare the results with 0 and (A.23), respectively.

Appendix B: Semianalytic approximation of graphs
For convenience, we approximated the numerical results by the function Tables 1, 2 below show the approximation coefficients for each graph in the normal and inverted hierarchies, as well as the coefficient χ 2 /n and and the number of independent points n.