Combining Experimental and Cosmological Constraints on Heavy Neutrinos

We study experimental and cosmological constraints on the extension of the Standard Model by three right handed neutrinos with masses between those of the pion and W boson. We combine for the first time direct, indirect and cosmological constraints in this mass range. This includes experimental constraints from neutrino oscillation data, neutrinoless double $\beta$ decay, electroweak precision data, lepton universality, searches for rare lepton decays, tests of CKM unitarity and past direct searches at colliders or fixed target experiments. On the cosmological side, big bang nucleosynthesis has the most pronounced impact. Our results can be used to evaluate the discovery potential of searches for heavy neutrinos at LHCb, BELLE II, SHiP, ATLAS, CMS or a future lepton collider.


Introduction
The Standard Model (SM) of particle physics and the theory of general relativity together can explain almost all phenomena observed in nature, covering scales that range from the extension of the observable universe down to the inner structure of the proton [1]. In the SM, matter is composed of fermions. All fermions except neutrinos are known to exist with both, left-handed (LH) and right-handed (RH) chirality. One reason to suspect that massive RH neutrinos exist is that these offer an explanation for the observed neutrino masses through the seesaw mechanism [2][3][4][5], hence they could explain neutrino flavour oscillations, the only confirmed evidence for physics beyond the SM observed in the laboratory to date. In addition they could explain a number of phenomena observed in outer space as well as a few unconfirmed experimental anomalies, for a review see e.g. Ref. [6]. Most noticeable, RH neutrinos can generate the observed baryon asymmetry of the universe (BAU), see [7] for a discussion, via leptogenesis [8], and are a natural Dark Matter candidate [9,10]. For a specific mass pattern they can even explain all of these phenomena simultaneously [11,12], see [13,14] for a detailed summary.
As for any new postulated particle, the experimentally most interesting properties of RH neutrinos N I are their masses M I and the interaction with other particles. The latter is characterised by a matrix of Yukawa coupling constants F αI with ordinary leptons of flavour α or, alternatively, mixing angles Θ αI . Unfortunately, the magnitude of the M I and F αI is entirely unknown. The main constraints come from neutrino oscillation experiments, which are only sensitive to a particular combination of them that determines the physical neutrino masses at low energies (essentially the ratio between the square of the coupling and the mass). To be specific, these are given by the eigenvalues m 2 i of the flavour matrix m ν m † ν , where m ν is given in Eq. (5). If the N I -interactions are to be described by perturbative quantum field theory, then the M I should be at least 1-2 orders of magnitude below the Planck mass. This can be estimated by inserting the observed neutrino mass differences into (5). On the lower end they can have eV (or even sub-eV) [15] masses, and any value in between is experimentally allowed. Cosmological constraints can be used to push the lower bound up. If one requires that the observed BAU is generated by leptogenesis, then past studies suggest that M I has to be larger than a few MeV [14,16]. This bound can be pushed up to roughly 100 MeV if one combines the negative result of past direct searches [17] with the cosmological requirement that the N I -lifetime is shorter than about 0.1s [14], though there are some uncertainties described in section 4.6. For longer lifetimes the energy released during the N I decay would affect the process of big bang nucleosynthesis (BBN) [18] and lead to an observable change in the abundance of light elements in the intergalactic medium (unless the N I are very long lived, in which case leptogenesis would not work). This still leaves a window 100MeV < M I < 10 15 GeV of roughly sixteen orders of magnitude open.
Hence, theoretical prejudice is the only guideline when picking a mass scale. In grand unified theories the RH neutrino mass scale is usually placed near (slightly below) the scale of grand unification. This choice is highly popular because it allows to explain neutrino masses with F αI ∼ 1 and leads to standard thermal leptogenesis. The disadvantage of this scenario is that the N I will most likely never be observed directly, though indirect constraints can be derived, e.g. from neutrinoless double β-decay [19] and cosmological considerations [20].
A discovery of the N I in the laboratory in the near future is only possible if the M I are at the TeV scale or below. Apart from the usual "lamp-post approach" argument of experimental testability, an appeal of this choice is that more massive RHNs would destabilise the Higgs mass [21], leading to a severe fine tuning problem in the absence of mechanisms that cancel radiative corrections such as supersymmetry. More particular motivations for RHNs below the TeV scale can be found in models based on classical scale invariance [22], in the framework of the "inverse seesaw" [23,24] and other models with approximate lepton number conservation (e.g. [25][26][27][28][29][30][31][32]), "no new scale"-considerations [33] or by applying Ockham's razor to the number of new particles required to explain the known beyond the SM phenomena [11]. Such a setup based on a minimal particle content is what we consider in the present work. More specifically, we study existing experimental constraints on N I in the (type I) seesaw model with M I below the mass of the W-boson. For a similar study with TeV masses we refer the reader to Ref. [34,35]. Experimentally this mass range is interesting because it is within reach of the LHC, b-factories or future colliders [36][37][38][39][40][41][42][43][44]. Cosmologically it is interesting because it allows to explain the BAU via leptogenesis from flavour oscillations during the thermal production of the N I in the early universe [11,45].
In principle the number n of RH neutrinos is a free parameter, as they are SM gauge singlets and not subject to any anomaly cancellation requirement. If the seesaw mechanism is the sole origin of neutrino masses, this implies n ≥ 2 because two non-zero mass differences between the light SM neutrinos have been observed, and the seesaw mechanism requires one RH neutrino per observed non-zero light neutrino mass. 1 Leptogenesis also requires n ≥ 2. Previous studies of experimental constraints are mostly focused on the case n = 2 (or even n = 1), either for simplicity or because the third neutrino is assumed to be a Dark Matter candidate that interacts so feebly that its influence on neutrino masses and leptogenesis is negligible. The known constraints for n = 2 can be found in Ref. [6,13,47] and references therein, recent works include [12,14,17,18,39,40,44,[48][49][50][51]. In this work we extend these past studies and assume that there are three RH neutrinos, in analogy to the number of generations in the SM.
The motivation to extend the analysis from n = 2 to n = 3 is two-fold. First, with in the seesaw framework we know that n ≥ 3 if the lightest neutrino turns out to be massive. Second, the range of N I parameters for which leptogenesis can be realised with M I ∼ GeV turns out to be very different for n = 2 and n = 3. For n = 2 the two masses masses M 1 and M 2 have to be quasi-degenerate in order to explain the observed BAU, see [52,53] and [11,12,14,16,[54][55][56]. This is also a necessary condition to measure the CP-violation in the sterile sector [57], which is responsible for the BAU in most of the parameter space [14]. Moreover, leptogenesis then requires the Yukawa coupling constants F αI to be too tiny to give measurable branching ratios in most existing experiments to find the N I [12,14,34,58]; small improvements of the existing bounds for M I < m K have recently been achieved at the E949 experiment [59] and may also be possible with NA62 [60] and at T2K [40] or LBNE [61]. These could, however, be found in dedicated future experiments [41] based on present day technology [12,14]. With n = 3 both of these restrictions can be overcome [42,62,63]; in this case leptogenesis does not require a mass degeneracy, and the parameter region where the BAU can be explained is within reach of existing experiments [42]. This provides strong motivation to re-analyse the known bounds to understand the parameter space for n = 3.
The purpose of the present work is to combine all known direct and indirect experimental constraints as well as bounds from BBN on RH neutrinos with masses between the pion and W-boson mass. This allows to compare the allowed parameter space to the cosmologically motivated region where leptogenesis is possible and the sensitivity of future searches. Our analysis includes indirect constraints from neutrino oscillation data [64], lepton flavour violation [65], neutrinoless double β-decay [66] and big bang nucleosynthesis [18]. As far as direct searches are concerned, we include bounds from the experiments PS191 [67], NuTeV [68] (both re-analysed in [17]), NA3 [69], CHARMII [70], DELPHI [71] (as given in [47]), LHCb [44] and BELLE [39]. They are stronger than those from violation of lepton universality [72][73][74][75], though these may become relevant in the near future [51].
The present article is organised as follows. In section 2 we recapitulate the seesaw model and introduce our notation. In section 3 we define four benchmark scenarios, characterised by the range of M I , which we will use to explore the parameter space. In the following section we recapitulate different bounds on this parameter space, coming from neutrino oscillation data (section 4.1), lepton flavour oscillation (section 4.2), neutrinoless double β-decay (section 4.3) and big bang nucleosynthesis (section 4.6). In section we numerically obtain combined constraints on the RH neutrino's mass and mixing in our benchmark scenarios . We discuss the implications of these constraints for cosmology and future experiments in section 6 and conclude in section 7.

The seesaw model
The well-known (type-I) seesaw model is defined by adding n neutral fermions ν R with RH chirality to the SM. These can couple to the SM neutrinos ν L in the same way as the RH and LH components of the charged leptons are coupled together, hence the name RH neutrinos. The Lagrangian reads Here flavour and isospin indices are suppressed. L SM is the SM Lagrangian, l L = (ν L , e L ) T are the LH lepton doublets in the SM and Φ is the Higgs doublet withΦ = ǫΦ * , where ǫ is the antisymmetric SU(2) tensor, M M a Majorana mass term for ν R with ν c R = Cν R T being its own charge-conjugate particle, and F is a matrix of Yukawa couplings. The charge conjugation matrix is C = iγ 2 γ 0 . We work in a flavour basis where M M = diag(M 1 , M 2 , M 3 ).
For M I ∼ GeV (or indeed any M I ≫ 1 eV) there are two distinct sets of mass eigenstates, which we represent by flavour vectors of Majorana spinors ν and N. The elements ν i of the vector are mostly superpositions of the "active" SU(2) doublet states ν L , they have light masses Here c.c. stands for the charge conjugation defined above. The elements are mostly superpositions of the "sterile" singlet states ν R and have masses of the order of M I . The only way how these interact with the SM at low energies is via this mixing with active neutrinos, which is characterised by the mixing matrix Θ with elements Θ αI ≪ 1. V ν is the usual neutrino mixing matrix and U ν its unitary part, V N and U N are their equivalents in the sterile sector. It is given by where M D ≡ F v and v is the temperature dependent Higgs field expectation value (v = 174 GeV at temperature T = 0). The matrix θ is related to Θ as θ ≡ ΘU T N . The unitary matrices U ν and U N diagonalise the mass matrices respectively. Since the eigenvalues of M M and M N coincide in very good approximation, we will in the following not distinguish between these and use the notation M I for both.
Here m diag can be parametrised by complex angles ω ij as The non-zero elements of the R (ij) are We work in the flavour basis where the charged lepton Yukawa couplings are diagonal, then the matrix U ν can be parametrised as (12) diag(e iα 1 /2 , e iα 2 /2 , 1).
Here U ±δ = diag(e ∓iδ/2 , 1, e ±iδ/2 ) and the non-zero entries of the matrices V are given by Here θ ij are the mixing angles amongst the active neutrinos, and α 1 , α 2 and δ are CPviolating phases. This parametrisation allows to encode all constrains from neutrino oscillation experiments in U ν and m diag ν . Another advantage is that (7) nicely factorises into a part 1 v U ν m diag ν that only contains parameters associated with the light active neutrinos and a rest that contains the parameters in the sterile sector. However, it is only valid as long as radiative corrections to the neutrino mass s m i are small because the tree level masses are used to parametrise F .

The benchmark Scenarios
Only five of the 7n − 3 parameters contained in (1) in addition to the SM for n = 3 have been determined experimentally (two mass differences |m 2 i − m 2 j | and three mixing angles θ ij ), see [64] for the recent values used in our analysis. With the CP-violating phases, three more can at least in principle be determined in low energy experiments involving the ν i only (though this is very challenging for the Majorana phases α 1 and α 2 ). This leaves us with a large number of free parameters in the sterile sector, making it very difficult to understand the parameter space analytically. This is in contrast to the case n = 2, in which there are only four parameters in the sterile sector, and further simplifications arise if one assumes the mass degeneracy required for leptogenesis. Therefore most of our following analysis is numerical, and we define a number of representative benchmark scenarios.
The high dimensionality also raises the question how to project the parameter space, i.e. plot our results. Since most of the parameters defined in the equations following (7) do not directly correspond to physical observables, it is not straightforward to identify physically meaningful exclusion regions in terms of these parameters. We therefore mainly express our results in terms of the physical masses M I and the mixings These are the quantities that enter branching ratios and rates observed in experiments. To define our benchmark scenarios we first specify an experimentally accessible region between M exp min and M exp max . The most promising experiments are LHCb, BELLE II and the proposed SHiP experiment. For all of these, the accessible region lies in the GeV range. There is no near-future experiment that could access masses between the B-meson mass and W-boson mass, though the perspectives at future colliders have been studied [43].
To define our benchmark scenarios, we also have to specify the total mass range M min < M I < M max within which we allow the M I to lie and which in general is broader than the experimental range. The bounds on the mixings U 2 αI of those N I that lie inside the experimentally accessible region may depend on the properties of those N J that have masses outside this region. As pointed out in the introduction, there are essentially no experimental constraints on M min and M max unless one adds cosmological considerations. For our numerical analysis we use M min = 100 MeV, motivated by BBN and direct search bounds. Note that this is not strictly required by the BBN bounds, not even if one requires successful leptogenesis. The BBN bounds for masses M I > 140 MeV are not known in detail. Here we require that the N I life-time is shorter than 0.1s. This is based on an extrapolation of the results presented in Ref. [18], which suggest that 0.1s is a realistic criterion for see also our discussion in Section 4.6. Here m π is the charged pion mass. For M I < 140 MeV, the limits turn out to be weaker [18] than this requirement, such that the BBN bounds that we impose may be too conservative in this regime. In principle even much lighter masses could be realised. It could be that only two N I contribute to leptogenesis, while the third one is lighter, very feebly coupled and has a lifetime that is comparable to the age of the universe. Then it would leave no trace in BBN or the cosmic microwave background because of its longevity and low number density, but can be a viable DM candidate [11]. Luckily the feeble interaction required to achieve this implies that the light N I would effectively decouple from the system and can be ignored as far as leptogenesis and the seesaw mechanism are concerned. Therefore Eq. (3) is a useful working assumption. We fixed M exp min = 100 MeV in our numerical analysis, but effectively Eq. (10) is the lower bound. Since we found almost no valid points below this value, our constraints are self-consistent, but possibly too conservative for small masses. This leaves the possibility open that M exp max < M J < M max . The main theoretical prejudice in our analysis remains the choice of M max . We will in all cases assume that its value lies below the W-mass and analyse how our conclusions depend on the choice. This brings us in the position to define benchmark scenarios.
Scenario A: The νMSM scenario -This scenario corresponds to the situation described after Eq. (3), i.e. one sterile neutrino N 1 is light and long-lived while the remaining two N 2 and N 3 are in the mass range between 100 MeV and 80 GeV considered here. This scenario has been described in Ref. [6,13,14,77] and references therein. We will not repeat or refine this analysis here. Note that the phenomenology of N 2 and N 3 is the same as in the model with n = 2 as far as neutrino masses, baryogenesis and direct searches are concerned because N 1 effectively decouples. Two aspects of this specific scenario should be recalled in the present context. First, baryogenesis requires a mass degeneracy between M 2 and M 3 at level ∼ 10 −3 [12,55]. The small mass splitting introduces a long oscillation time scale, which can give rise to resonant phenomena that require special treatment in the early universe [78][79][80][81][82] and possibly in the laboratory [83,84]. There is no specific reason for a mass degeneracy in general n = 3 scenarios. In this case standard methods pertinent to short oscillation time-scales can be used to describe the propagation of the N I even in the early universe [62,[85][86][87]. Second, in this scenario there exists a strict lower limit on the mixing of N 2 and N 3 with active neutrinos, which is a function of their mass [14,17,48,56]. This is not only related to the generally smaller parameter space in an (effective) n = 2 model, but also to the previous requirement M 2 ≃ M 3 . We note that in combination these two points imply that N 2 and/or N 3 need to have sizable mixings with active neutrinos in order to explain the observed neutrino masses. In more general scenarios it could always be that the N I with M exp min < M I < M exp max have a tiny (or zero) mixing and neutrino masses are generated by heavier Scenario B: The optimistic SHiP scenario -This scenario is defined by the choice For this choice the masses of all N I lie in the region where the SHiP experiment will have its highest sensitivity in U 2 αI , estimated to be ∼ 10 −9 . This parameter region is interesting for leptogenesis because the BAU can be explained without significant tuning in the ratio U 2 eI /U 2 µI [42]. The sensitivity drops considerably for M I larger than the D-meson mass m D = 1, 869.61 ± 0.1 MeV [1] because the N I are mainly produced in decays of D-mesons. Other experiments (such as LHCb and BELLE) also have their main sensitivity in this mass range, but it is unlikely that they can significantly improve bounds from past searches [59,[67][68][69][70][71], see also [17,47].

Scenario C: The optimistic b-factory scenario -This scenario is defined by the
For this choice the N I can all be produced in B-meson decays (with the charged B-meson mass m B = 5279.26 ± 0.17 MeV). This makes a detection of the N I in B-meson decays at LHCb or BELLE II possible. While the current bounds from these experiments [39,44] in most of the mass range are still weaker than the old bounds from DELPHI [71], they will soon become the front line in the exploration of the cosmologically motivated parameter space [42].
Scenario D: The general GeV seesaw at b-factories and SHiP -It is therefore instructive to analyse the perspectives of an discovery at b-factories are if only part of the N I mass spectrum is accessible to direct searches. Amongst possible near-future searches, SHiP would be the most powerful one. We define a benchmark scenario by That is, we assume at least one M I is in the range accessible to SHiP, LHCb and BELLE II and investigate what are the known experimental constraints on this particle's properties if we allow its siblings to have masses anywhere between 100 MeV and m W . During our parameter scan this is practically realised by randomising M 1 in the interval 100 MeV does not introduce any bias because all other parameters are being randomised without preference for any value in the parametrisation (7).
Scenario E: The general GeV seesaw Finally, we explore the possibility that we can test the entire mass range up to m W experimentally, which in the near future is only possible at ATLAS and CMS. So far both collaborations have published constraints only for M I > m W [37,38]. For M I < m W the experimental perspectives are thought to be much better because the N I can be produced in gauge boson decays. Even further improvement could be made at a future circular collider [43]. We define the last scenario with With the present work, we initially intended to study the potential to discover the N I in meson decays, hence scenarios A-D are focused in the mass range M I < m B . We were made aware of the potential of CMS to probe mixing angles as small as U 2 I ∼ 10 −9 [88] only after large parts of our analysis had already been finished. Because of this we added the scenario E, but only completed the analysis for normal hierarchy and m 1 = 0 in order to avoid further delay of the publication (while for all other scenarios we investigated both choices of neutrino mass hierarchy and two choices of the absolute mass scale). Given the experimental perspectives, scenario E certainly deserved further investigation in the future.

Experimental and observational constraints 4.1 Neutrino oscillation data
Neutrino oscillation data constrains the differences between the neutrino masses m i as well as the mixing angles and phases in U ν . It also imposes constraints on the mixings U 2 αI . That is, once the m i and U ν are fixed, not all values for the U 2 αI can be realised. On one hand this is because at least some of the N I have to have a significant mixing with active neutrinos to yield eigenvalues for m ν m † ν that are large enough to explain the observed On the other hand the U 2 αI cannot be too large. Otherwise radiative corrections would lead to physical neutrino masses that are too large even if the tree level m i are small.

Minimal mixing
An obvious question to pose in any of the scenarios defined above would be What is the minimal mixing U 2 eI with ν e that N I with M exp min < M I < M exp max can have?
Unfortunately the answer is U 2 eI = 0 (and same for mixing with ν µ ). This statement in some sense is trivial: Since only two neutrino mass differences have been observed, we can explain all observational data with n = 2 (in which case the lightest neutrino is massless). This is e.g. effectively realised for n = 3 with U 2 e1 = U 2 µ1 = U 2 τ 1 = 0, hence this scenario must be allowed. Even if we require that all light neutrinos as massive, we can still choose the remaining parameters such that individual U 2 αI vanish. For illustrative purposes we consider the special case of vanishing phases δ = α 1 = α 2 = 0 and ω 23 = ω 13 = 0 and normal hierarchy. Then all mixings Θ α3 ∝ m 1 , i.e. N 3 entirely decouples if the lightest neutrino is massless (irrespectively of the choice of the M I -spectrum). For ω 12 = −θ 12 one can in addition suppress Θ e2 up to corrections of order m 1 /m 3 , i.e. the ratio between the lightest and heaviest neutrino. We have checked analytically that there are choices for the ω ij for which a given Θ eI = 0 or Θ µI = 0. There exist analytic solutions for ω ij that set the coupling of two of the N I to electrons exactly to zero, but the expressions are rather lengthy and not illuminating. Hence, we cannot impose a lower bound on the mixing of any individual N I . These findings are consistent with what has been described in [56]. Leptogenesis can also provide a lower bound on U 2 I [14,42], but this relies on the extra assumption that the N I should explain the BAU and depends on whether the seesaw model (1) is embedded into a bigger extension of the SM, hence we do not consider it here.
Of course, some of the N I need to have a significant mixing with active neutrinos to explain neutrino masses. If we demand that all masses M I lie within the experimentally accessible region (as we do in scenarios B, C and E), then we can pose the question Assuming N 1 is the sterile neutrino with the largest mixing with ν α of a given flavour α (i.e. U 2 1α > U 2 2α , U 2 3α ), what is the smallest value that U 2 1α can take as a function of M 1 ? 2 The answer to this question is some non-zero U 2 αmin , which is a function of M 1 . Its meaning is the following: Assume that we have an experiment that covers a mass range between M exp min and M exp max , and that all three N I have masses in this interval. Then we are guaranteed to find at least one sterile neutrino if we push the sensitivity in U 2 αI below U 2 αmin across the whole mass range. Otherwise we have excluded the model. Unfortunately also this approach is of limited use unless at least two of the M I are within reach of the experiment, as otherwise the N I with sizable mixing could simply be too heavy to be seen. Hence, in order to derive a non-zero lower bound this way, we would have to impose an upper bound on two of the M I by hand, and this upper bound would have to be relatively low (below M exp max ). There is no reason why this mass spectrum should be realised in nature, and a lower bound that is based on this condition is pretty much meaningless.
We can, however, put a lower bound on the quantity without extra assumptions. If all three m i are non-zero, then a lower bound on U 2 I is enforced by neutrino oscillation data: To generate three non-zero neutrino masses, all three N I have to mix with the active neutrinos. 3 This means that each N I must mix with at least one active flavour, meaning that U 2 I is bound from below for given M I . If the lightest neutrino is massless, then there is no such lower bound. Two non-zero m i can be generated by the seesaw mechanism with two sterile flavours; the third sterile flavour is allowed to have an arbitrarily small mixing (even zero). Hence, for given M I there is no lower bound on U 2 I . In this case U 2 I is bound from below only by cosmology, coming from the requirement to have a lifetime of less than 0.1s described in section 4.6. This bound is unavoidable in the framework of the minimal seesaw. In the n = 3 scenario BBN can only constrain U 2 I , while in the n = 2 scenario it also allows to derive bounds on the individual U 2 αI , as too large hierarchies amongst them are forbidden by neutrino oscillation data.

Maximal mixing
The active neutrino masses (5) are generated at tree level, but also receive radiative corrections from loops with virtual N I . If all entries of the neutrino mass matrix are of the same order, then these corrections are always negligibly small because the Yukawa couplings F αI are all of the order of the "naive seesaw expectation" F 0 ≡ m atm M M /v 2 , i.e. the wouldbe value of the coupling constant in a world without flavour, n = 1 and a neutrino mass of the order m atm . It is, however, possible that individual F αI are much larger than F 0 and the smallness of the m i is the result of a cancellation amongst their contributions to the neutrino mass matrix. Such a cancellation could either be accidental or due a symmetry, such as an approximately conserved lepton number. The existence of such cancellations at tree level does not guarantee that these persist after loop corrections, in particular there could be sizable radiative effects δm ν to m ν if individual F αI are much bigger than F 0 [89]. Therefore one can obtain an upper bound on the F αI (and hence U 2 αI ) from the requirement that the correction to the observed ∆m 2 sol and ∆m 2 atm from the one-loop contribution to m ν , given by [90] δ(m ν ) αβ = is smaller than the uncertainty quoted in Ref. [64]. Note that in this procedure, we first feed the observed neutrino mass differences into the parametrisation (7) of F and then demand that radiative corrections are small enough not to spoil the agreement with experiment. Using this method, we discard potentially viable points where the tree-level and one-loop contributions add to a mass matrix that is consistent with observation, even though the radiative contributions are large. Nonetheless, this approach is pragmatic since an efficient parametrisation of the see-saw model that fully captures the one-loop corrections is not available, and we should expect that this shortcoming only affects a small and fine-tuned region of parameter space.

Lepton flavour violation
The violation of lepton flavour in the neutrino sector also mediates flavour violation amongst the charged leptons. Contributions come from both, light and heavy neutrinos. This allows to constrain the N I parameter space by looking for rare processes. Amongst these, the experimentally most constrained is the decay µ → eγ, with a branching ratio B(µ → eγ) < 5.7 × 10 −13 [65]. The predicted branching ratio in the seesaw model (1) is given by [91,92] and the loop function We use the condition of consistency with the experimental upper bound as a criterion to constrain the U 2 αI . At this stage we have not implemented data from any other search for lepton flavour violation, which is usually thought to be less constraining [34]. The nonobservation of µ → eγ imposes rather strong constraints on U 2 αI in models with M I ∼ TeV [34]. For GeV masses, however, the bounds are known to be much weaker already for n = 2 [93,94]. Since n = 3 offers more freedom, it can be expected that the experimental limits are even less constraining. This is confirmed by our analysis, and it is easy to understand this qualitatively. Let us first consider the difference to the TeV seesaw studied in [34]. Our analysis shows that for large U 2 αI , Eq. (18) is dominated by the contribution from N I exchange. Each individual contribution from the light states ν i is much larger than those from the N I , but they cancel each other since the loop function is approximately constant, the light neutrino contribution is given by (V ν V † ν ) eµ . While individual |(V ν ) αI | 2 are of order one, the combination (V ν V † ν ) eµ is expected to be of order 10 −17 if one inserts the neutrino oscillation data reported in Ref. [64] and assumes GeV scale RHNs. While the loop function (19) only increases by a factor of order one if one decreases M I from TeV to GeV, the mixing angle θ ∼ M D /M M scales as θ 2 ∼ m i /M I even with the "naive seesaw" assumption F 0 , so we expect the contribution from N I exchange in Ref. (18) to be three orders of magnitude weaker. Moreover n = 3 leaves more room for cancellations.

Neutrinoless double β-decay
In the seesaw model (1) neutrinos are Majorana particles, and the Majorana mass term M M does not only violate lepton flavour, but also total lepton number. This makes neutrinoless double β-decays possible. The experimental constraints on neutrino properties are conventionally expressed in terms of the effective Majorana mass, which is defined as Here we follow Ref. [34] and use f A = (M A /M I ) 2 F A , M A ≃ 0.9 GeV and F A = 0.079 for 76 Ge, based on previous estimates in Refs. [19,95,96]. We can impose upper bounds on the U 2 eI by requiring that the rate for neutrinoless double β-decay is consistent with the experimental upper bound m ee < 0.2 eV [66]. Though the process is almost always dominated by the exchange of light neutrinos [97], this still imposes constraints on the seesaw model, especially in the lower mass region of our analysis.

Direct searches
In the present work we are interested in N I with M I < m W . In this case these particles can be produced in the decay of weak gauge bosons [98,99], as it has e.g. been looked for at the DELPHI experiment, which currently imposes the strongest bounds for m B < M I < m Z [71]. These could be improved at a future electron collider [43]. LHC so far has only published bounds for M I > m W [37,38], and being a hadron collider, it is not ideally suited to look for N I [12,14,34,58] due to backgrounds. For masses M I < m B the N I can also be produced in meson decays. Depending on the detector design, there exist different search strategies, either looking for missing 4-momentum or aiming to observe the N I -decay within the detector, see e.g. [6,13,36,40,47]. In the region m D < M I < m B the strongest bounds come from DELPHI [71], LHCb [44] and BELLE [39]. The latter will improve with the LHC's 14 TeV run and BELLE II, exploring parts of the leptogenesis parameter space [42]. For masses below m D the existing bounds from PS191 [67], NuTeV [68], NA3 [69] and CHARMII [70] are much stronger. In the lower mass range M I < m K new results have recently been published by the E949 collaboration [59]. In this mass region even small improvements are significant, especially for the n = 2 scenario or scenario A (which effectively is n = 2). The reason is that in this case BBN can impose a lower bound on U 2 µI (see section 4.6), which is not much below the upper bound from direct searches, such that the remaining window may be closed in the near future. Further improvement may also be possible with NA62 [60] and at T2K [40]. The most significant improvement for M I < m D (and possibly even for slightly larger masses) could be made with the proposed SHiP experiment [41]).
The interpretation of the experimental results is not always trivial because some collaborations made assumptions about the interactions of heavy neutrinos that are different from those of the Lagrangian (1), see Ref. [17] for a discussion. Moreover, in some cases it was assumed that the particles are Dirac fermions (while the N I are Majorana particles). For our analysis we can nevertheless directly use the bounds obtained for Dirac fermions: For Majorana particles the number of spin states (and hence the number of particles produced in meson decays) is half as large as for Dirac particles, but they can decay via charge conjugate processes (hence the decay is twice as fast), such that one gets same number of  Figure 1: Constraints on U 2 eI from the experiments DELPHI [71], L3 [100], PIENU [101], TINA [102], PS191 [67], CHARM [103], NA3 [69], and kaon decays [104]. For peak searches below the kaon mass we used the summary given in [47], for PS191 we used the re-interpretation given in [17]. events in the detector. Another issue is that different collaborations quote bounds with different confidence levels, which are partly summarised in [47]. In principle this should be modelled with a smooth likelihood function. We entirely ignored this issue to simplify the problems and treated all bounds as "hard", i.e. as if all forbidden points are forbidden with absolute certainty and all allowed points can be realised in nature with exactly the same probability. We display a summary of the existing direct search bounds in figures 1-3. They are also shown as red lines in figures 4-16.

Lepton Universality and electroweak precision data
If the N I are light enough to be produced in the decays of charged mesons, this can not only lead to directly observable signals discussed in section 4.4; the existence of these decay channels also modifies the decay rates into SM leptons [98,99]. While the individual decay rates are difficult to calculate precisely, most of the uncertainties cancel out if one forms their ratio. This is known as lepton universality. Currently the bounds imposed by this measurement [51,[72][73][74][75] are not competitive with direct searches as far as U 2 eI and U 2 µI are concerned, and we do not include these in our analysis. Universality constraints on U 2 τ I are stronger than the direct ones shown in Figure 3, see e.g. [47], where upper bounds ∼ 6 × 10 −3 are obtained. However, we do not find any points with such large U 2 τ I in our scan, which seems to justify neglecting lepton universality bounds. They should, however, be included in future studies because new measurements, for example at the NA62 , LHCb [44], BELLE [39], BEBC [105], FMMF [106], E949 [59], PIENU [101], TINA [102], PS191 [67], CHARMII [70], NuTeV [68], NA3 [69] and kaon decays in [104,107]. For the bounds from kaon decays we used the interpretation given in [47,108]. For NA3, BEBC and FMMF we used the estimates from [47]. For PS191 we compare the results according to the interpretation in Ref. [17] (solid line) as well as [59] for two different channels (dashed and dotted line).  Figure 3: Bounds on U 2 τ I based on the interpretation of CHARM data in [109], NOMAD [110] and DELPHI [71]. experiment, may change this situation [51]. The situation with respect to CKM unitarity and other electroweak precision tests is similar [47].

Big Bang nucleosynthesis
All constraints reviewed in the previous Sections come from laboratory experiments. In addition, the lifetime of the RH neutrinos is constrained by the requirement that they decay before BBN. If they decay during or shortly before BBN, then the decay products have energies of the order M I , much larger than the temperature of the primordial plasma at that time. This would affect the abundances of light elements that are produced, either by directly dissociating nuclei that have already formed or indirectly by causing a deviation from thermal equilibrium in the plasma that modifies the rates for the processes involved in nucleosynthesis. This can be used to impose a lower bound on the U αI .
Here we adopt a simple (but rather conservative) criterion and demand that all N Ilifetimes are shorter than 0.1s. To calculate the lifetime, we use the decay rates given in Refs. [36] and [14]. This imposes a lower bound on U 2 I that is unavoidable in the framework of the minimal seesaw. It can be softened in extensions of (1) if there are decay channels for N I other than via mixing with active neutrinos. For M I < 140 MeV, where this effect is most important because the N I -lifetime scales as M 5 I , this effect has been studied in detail in Refs. [18,111,112]. Their results show that the actual bounds for M I below the pion mass are weaker than our criterion, while for larger masses our approach seems justified as a first approximation. As already anticipated in (10), this implies that our results can only be applied to N I that are heavier than pions. Due to the weaker BBN bounds there may be allowed parameter regions in the mass range M I < m π , which seem excluded according to Figures 4-16. We postpone a more detailed study of this mass range to the future.
For n ≥ 3 BBN cannot impose bounds on the mixings U 2 αI of N I with individual active flavours. A given U 2 αI can be arbitrarily small (even zero) if another U 2 β =αI is sufficiently large to ensure a lifetime below 0.1s. The situation is different for n = 2 (and hence in our scenario A). In this case there are far less free parameters in F ; in particular, there is only one complex "Euler angle" ω in R that governs the mixing of N I with all active flavours, hence U 2 eI , U 2 µI and U 2 τ I cannot be too different in size. This allows to translate the BBN bounds into constraints on individual U 2 αI and to correlate the upper bounds on the individual U 2 αI to obtain a stronger upper bound on U 2 I . For example, in the mass range M I < m K and for n = 2 the lower bounds on U 2 eI and U 2 µI also strongly constrain the poorly measured U 2 τ I , allowing to impose an upper bound on U 2 I that is not too far from the lower BBN bound [51], and this window may soon be closed. For the n = 3 scenario under investigation here this is not possible because the three angles ω ij offer more freedom to change the size of the different U 2 αI . This is also the reason why leptogenesis with n = 3 does not rely on a mass degeneracy [62] and can be achieved for much larger U 2 αI [42].

Numerical results
In our analysis, we fix the mixing angles θ ij in U ν and the mass splittings m 2 i −m 2 j to the best fit values of the global fit presented in [64]. For each of the benchmark scenarios B, C, D and E we independently study the cases of "normal hierarchy" (m 2 1 m 2 2 < m 2 3 ) and "inverted hierarchy" (m 2 3 < m 2 1 m 2 2 ). Since the absolute scale of neutrino masses is not known, we focus on the two extreme cases that the lightest neutrino is massless and has a mass of 0.23 eV, the upper limit from cosmology [113], studying both situations independently. We randomise all remaining parameters in a numerical scan, where we used the parametrisation defined in Section 2 and restrict Imω ij to the interval −8 < Imω ij < 8. The dependence of the physical parameters on δ, 4 α 1 , α 2 and Reω ij is periodic. Each of the 13 scans consists of 5 × 10 8 randomly generated parameter choices. Amongst all these points in parameter space we keep only those that are consistent with the indirect constraints from neutrino oscillation data [64] 5 , non-observation of µ → eγ decays [65], neutrinoless double β-decay [66] and direct search constraints from the experiments summarised in Figures 1-3. In addition, we require that the lifetime of each of the N I is shorter than 0.1s to be consistent with constraints from BBN. For this we employ the N I decay rates calculated in [14,36]. Note that this is a conservative estimate, and a more detailed analysis [18] might relax the BBN bound a bit. Finally, we required that all diagonal entries (F † F ) II M I /(16π) are smaller than the smallest mass splitting amongst the heavy neutrinos, so that they are well defined individual particles. The highly mass-degenerate case in which this is not given is not excluded and can lead to the interesting resonant leptogenesis scenarios [53,[78][79][80][81][82], but the propagation in the laboratory may require a more sophisticated treatment [36,83,84]. Amongst the points that are consistent with all these constraints we identify those that maximise the mixing angles U 2 eI , U 2 µI , U 2 τ I and their sum U 2 I = U 2 eI + U 2 µI + U 2 τ I in order to determine upper bounds for these quantities. In addition, we identify the lower bound on U 2 I . As discussed in Section 4.1.1, it is not possible to impose a lower bound on the individual U 2 αI . The results of the scans are presented in  Present neutrino data already allows to identify a preferred region for δ [64]. However, since the significance of this preference is still low, we treated all values of δ as equally likely. 5 The lower bounds on U 2 αI from neutrino oscillation data are automatically imposed when using the Casas-Ibarra parametrisation, but the upper bound from radiative corrections (16) has to be imposed by hand.

Discussion
As expected, combining all known direct and indirect bounds self-consistently leads to stronger bounds than simply superimposing them in the M I − U 2 αI plane.

General remarks
Validity of the numerical approach -Since our approach is based on a stochastic scan, it is not possible to prove that the extreme points we obtain are indeed strict bounds.
• It could always be that there exist "funnels" or small isolated regions in parameter space that contain even more extreme allowed points. How narrow or extreme they appear of course depends on our parametrisation. This problem in particular arises for the upper bounds.
• As usual, we use the numerical convergence, along with the known parameter dependencies of the indirect bounds, as a criterion to judge whether there exists a strict upper bound near the maximal mixings we find, or whether the parameter space volume in some direction is simply shrinking so rapidly ("funnel") that the scan fails to find the "tuned" allowed parameter region leading to even larger mixing.
• For instance, in several scenarios there appears to be a small number of points near M I ∼ 100 MeV that suggest an upper limit on U 2 eI or U 2 µI that is much stronger than direct search bounds. We are almost certain that these are not valid upper bounds, i.e. for U 2 eI and U 2 µI there exist allowed points with mixings all the way up to the direct search bounds (red line) in this mass region.
• We use a scanning strategy in which the parameters are randomly chosen without correlation to each other or the previous choice. In order to explore narrow "funnels", an adaptive Monte Carlo approach would certainly be more promising. On the other hand, such method might be more likely to miss out on isolated regions. A comparison of different scanning techniques would be interesting in the future.
Parametric dependencies -There are a few simple rules to understand the parameter dependence of the mixings U 2 αI for given M I . For n = 2 these have e.g. been studied in detail in [17,34,48,97]. For n = 3 they are more complicated, but general tendencies can still be identified easily.
• In general, the dependence on the parameters δ, α 1 , α 2 and the Reω ij is relatively weak (though their precise value does matter if one aims to construct very special points, e.g. set one mixing U 2 αI exactly to zero).
• As already pointed out in [17,42,48,55], the magnitude of the active-sterile mixing is mostly controlled by the Imω ij because the elements of F αI depend exponentially upon them. Small |Imω ij | < 1 tend to give values for F αI near the "naive seesaw value" F 0 , typically leading to small U 2 αI and small branching ratios in experiments. This situation is considered "natural" by many authors, though this statement of course depends on the parametrisation and an arbitrary definition of "naturalness".
• Large |Imω ij | > 1 lead to comparably large individual U 2 αI . This makes it easier to find the N I experimentally. The smallness of the m i in this case is not (only) due to the smallness of individual elements on m ν , but requires cancellations amongst them that lead to small eigenvalues m 2 i of m ν m † ν . While arbitrarily large |Imω ij | are allowed by the seesaw relation (5) at tree level (which is only sensitive to a particular combination of the Θ αI ), different individual U 2 αI = |Θ αI | 2 are constrained by (16), (17) and (20). Strict bounds on the U 2 αI are of course obtained from direct searches, see section 4.4, but these are well-known, and it is the goal of this work so explore whether stronger bounds can be imposed by combining them with indirect constraints.

Summary
Neutrino mass spectrum • The choice of M max appears to have no significant effect on the upper bounds on U 2 αI . In the region where the mass range of our benchmark scenarios overlaps there is no difference between them. This can easily be understood qualitatively. The direct search bounds are independent of the choice of benchmark scenario. The indirect constraints from lepton number and flavour violation generally forbid large mixings unless there is a (fundamental or accidental) symmetry in Θ that leads to cancellations in (17) and (20), allowing to keep them small in spite of relatively large individual Θ αI . This is most easily realised if two of the N I have quasi-degenerate masses, then they can e.g. form a pseudo-Dirac spinor. The maximal mixings Θ αI that a given N I with fixed M I can have can be found in this parameter region. If the mass splittings are large, the entries of Θ tend to correspond to the "naive seesaw expectation". Therefore the maximal active-sterile mixings are rather independent of M max .
• The absolute neutrino mass scale affects the lower bound on U 2 I , as discussed in section 4.1.1.
• The combination of absolute neutrino mass scale and neutrino mass hierarchy also affects the upper bounds on the mixings. In the mass interval M I < m D , where direct search bounds dominate the combined constraints, the dependence on the absolute neutrino mass scale is most visible in the bounds on U 2 τ I because this element is much less constrained by direct searches, and the way how indirect bounds can be obtained from combining the known constraints depends on the m i . If the lightest neutrino is massless, then there exists a direction in sterile flavour space (a linear combination of the N I ) that does not couple to active flavours, i.e. neutrino masses are only generated by two linearly independent sterile flavour eigenstates. This means that one can express the U 2 τ I in terms of U 2 eI and U 2 µI , and the direct search bounds on these translate into a bound on U 2 τ I . It is dominated by the (weaker and hence more relevant) upper bounds on U 2 µI . This effect is stronger for inverted hierarchy. If the lightest neutrino is not massless, then there is more freedom. However, the remaining indirect bounds prohibit arbitrarily large Imω ij , therefore the bounds on U 2 µI still leave an effect on those on U 2 τ I . The bounds on U 2 I also depend on the absolute neutrino mass scale because weak bounds on U 2 τ I automatically imply weak bounds on U 2 τ I . For a massless lightest neutrino and normal hierarchy the combined bounds on U 2 eI seem stronger than the direct search bounds due to constraints from neutrinoless β-decay; for a high absolute mass scale and/or inverted hierarchy the bounds on U 2 eI for M I < m D coincide with the direct search bounds.
• In the interval m D < M I < m B the constraints from neutrinoless double β-decay have an even stronger effect on U 2 eI ; for the other flavours the combined bounds more or less coincide with the direct search bounds. If the lightest neutrino is massless, then the combined bounds are stronger than the direct search bounds, tough it is not entirely clear whether they impose a strict upper bound or simply reduce the viable parameter space volume. For normal hierarchy this effect is more prominent in the entire interval m D < M I < m B . For a massive lightest neutrino the combined bounds generally coincide with the direct search bounds, though it seems easier to get large mixings for inverted hierarchy.
Total mixing U 2 I • In the minimal seesaw model (1), the requirement for the N I to decay before BBN imposes a strict lower bound on U 2 I . This bound is very crucial for M I < m K , where it is not too far away from the upper direct search bound. This mass region could possibly be excluded in the near future. For m K < M I the BBN bound is far less important because it rapidly decreases with increasing M I due to the ∝ M −5 I scaling of the N I lifetime.
• If the lightest neutrino is massive, then neutrino oscillation data also imposes a lower bound on U 2 I , see section 4.1.1. If the absolute neutrino mass scale is near the cosmological limit 0.23 eV, then the lower bound in the region M I < m D is larger than 10 −10 , i.e. not far away from the anticipated SHiP sensitivity.
• If one adds the requirement to explain the observed BAU by leptogenesis, then this imposes an upper as well as a lower bound on U 2 I . These have been studied in detail in Scenario A [12,14,16] and partly been investigated in scenario C [42]. A detailed investigation of the lower bound would be particularly interesting for the case that the lightest neutrino is very light and no lower bound can be obtained from neutrino oscillation data. eI with electron flavour • Very strong constraints on U 2 eI have been derived from neutrinoless double β-decay shown in some past works [19,47] under the assumption that either there is only one sterile neutrino (n = 1) or all M I are degenerate. For n = 3 they are much weaker. 6 Though very large U 2 eI are ruled out once µ → eγ and/or direct search constraints are imposed, the resulting final bound is still much weaker than the one shown in [19,47].
• The upper bounds for M I < m D essentially coincide with the direct search bounds, except for the case of massless lightest neutrino and normal hierarchy discussed above.
Here SHiP would be the ideal experiment to search for N I .
• For m D < M I < m B the constraints from neutrinoless double β-decay and µ → eγ are stronger than those from direct searches if the lightest neutrino is massless, as discussed above.
• For m B < M I the constraints from neutrinoless double β-decay become less severe, but those from neutrino oscillation data are increasingly important, as the suppression of radiative corrections ∼ M 2 I /v 2 in (16) becomes inefficient for M I 10 GeV. This imposes an upper bounds that scales as ∝ M 2 I /v 2 , but is difficult to be resolved numerically because the phase space already shrinks exponentially in the vicinity of the bound.
Mixing U 2 µI with muon flavour • In almost the entire region below m B the upper bound can be obtained by simply superimposing the known direct search bounds. The only exception is a small mass interval near M I ∼ m D . In this region U 2 µI is indirectly constrained by the CHARM bounds on U 2 eI and their combination, which extend to slightly higher masses than the NuTeV bounds on U 2 µI itself.
• The lack of indirect constraints in the region m D < M I < m B (where they exist for U 2 eI ) is due to the fact that neutrinoless double β-decay is not directly sensitive to U 2 µI . The numerical convergence is excellent if the lightest neutrino is massless. If it is massive, the convergence is worse.
• For M I > m B the behaviour is similar to that of U 2 eI ; it is governed by the upper bound from the requirement to keep radiative corrections to m i small, which becomes increasingly efficient for M I 10 GeV.
Mixing U 2 τ I with tau flavour • For n = 2 the lower bounds on U 2 eI and U 2 µI strongly constrain the poorly measured U 2 τ I , allowing to impose an upper bound on U 2 I that is not too far from the lower BBN bound. This appears to already rule out the mass region below M I ≃ 240 MeV if one takes the BBN bounds literally [51]. The reason is that the parameter space for n = 2 is sufficiently small that neutrino oscillation data prohibits large hierarchies between the mixings of a heavy neutrino with different active flavours, see section 4.6. This is not the case for n = 3: The upper bound on U 2 I lies above the lower bound from BBN even for M I < 240 MeV because U 2 τ I can be comparably large. We find that these bounds meet each other only below M I ∼ 100 MeV; a region, where our treatment of BBN bounds is too conservative. Hence, we cannot deduce a lower bound on M I .
• In spite of this, combining all direct and indirect search bounds does allow to impose stronger constraints on U 2 τ I than any individual measurement. For M I < m D these U 2 τ I are orders of magnitude stronger than the known direct constraints. They come from the combination of neutrino oscillation data and direct search bounds on U 2 eI and U 2 µI . • For M I > m B the behaviour is very similar to that of U 2 eI and U 2 µI ; it is governed by the upper bound from the requirement to keep radiative corrections to m i small.

Conclusions
We have combined various direct and indirect constraints on the type-I seesaw model with three RH neutrinos and Majorana masses between the pion mass and W-boson mass. This extends previous studies, which were mostly performed under the assumption that there are only one or two RH neutrinos, and includes updated neutrino oscillation data as well new experimental data from the experiments GERDA, MEG, LHCb, BELLE and the E949. The main findings of our numerical analysis are the following.
• If the absolute neutrino mass scale is near its upper cosmological limit, then the lower bound on U 2 I from neutrino oscillation data for M I < m B is roughly 10 −10 , i.e. near the anticipated SHiP sensitivity. For M I ≫ m B it is roughly 10 −12 . If the lightest neutrino is massless, then the only lower bound on U 2 I comes from BBN and is much weaker.
• In comparison to the n = 2 scenario, for n = 3 the combination of neutrino oscillation data and direct search bounds on U 2 eI and U 2 µI leads to weaker constraints on U 2 τ I and consequently U 2 I . The combined bounds on U 2 τ I for M I < m D are nevertheless orders of magnitude stronger than the direct search bounds.
• As a result of these weaker constraints on U 2 τ I , the combined upper bound on U 2 I is above the lower bound from BBN for all masses we consider, and we cannot impose a lower bound on M I . This is in contrast to what is found for the n = 2 scenario, in which M I can be constrained from below if the same BBN criterion is applied as in the present work (lifetime < 0.1 s) [51]. It should, however, be kept in mind that the BBN constraints suffer from various uncertainties and may be too conservative.
• The combined bounds on U 2 eI in the mass range M I < m B coincide with the known direct search bounds if the absolute neutrino mass scale is relatively high. If the lightest neutrino is massless, then the combination of indirect bounds (in particular from neutrinoless double β-decay) and all direct bounds imposes a slightly stronger constraint on U 2 eI than direct searches alone. However, the constraints from neutrinoless double β-decay on U 2 eI are much weaker than in the n = 1 or mass degenerate scenarios studied in Refs. [43,114].
• The bounds on U 2 µI in the entire mass range M I < m B coincide with the known direct search bounds except in a small region near M I ∼ m D , where the combined bounds are much stronger than the direct search bounds.
• For masses M I ≫ m B neutrino oscillation data also disfavours large mixings U 2 αI , which would lead to large radiative corrections to the neutrino masses. For M I > 10 MeV this gives a stronger constraint stronger than direct searches in DELPHI and L3. This applies to all flavours. The reason is that cancellations in the tree level mass matrix (5) in general do not keep loop corrections (16) small unless there is an underlying fundamental symmetry in the Lagrangian, such as approximate lepton number conservation. It could, however, be that there exist specific regions in which such a symmetry is realised and allows for larger mixings, and that our scanning method failed to find these regions.
• Low scale leptogenesis is one key motivation to search for heavy neutrinos with masses below the electroweak scale. Previous studies that aimed to identify the parameter region where leptogenesis is possible in scenarios A [12,12,16] and C [42] included several indirect constraints, but did not incorporate all direct search bounds self-consistently. The latter were simply superimposed in the mass-mixing plane after the completion of the parameter scan. It would be interesting to perform a parameter scan that consistently incorporates all constraints we discussed here and the requirement to explain the observed BAU. This may reduce the viable leptogenesis parameter space in comparison to previous estimates. Such a study would be crucial to evaluate the perspectives of b-factories, the proposed SHiP experiment or future circular colliders to probe low scale leptogenesis scenarios.
While our analysis represents the state of the art in the given mass range, there are several improvement that can be made in the future. Given the expected experimental progress, a future analysis should incorporate more observables that are sensitive to lepton number and lepton universality violation and cover the experimentally interesting mass range M I < m π . A measurement of the absolute neutrino mass scale or the neutrino mass hierarchy would greatly reduce the uncertainty in the bounds on active-sterile mixings. To a lesser degree, also a measurement of the Dirac phase δ in the PMNS matrix would help to constrain the parameter space. On the cosmological side it would be desirable to better understand the precise bounds on sterile neutrino lifetimes from BBN in the mass range M I > m π . At a technical level, the use of more sophisticated Monte Carlo methods would help to determine the precise upper bounds on active-sterile mixings in parameter regions where highly tuned "funnels" exist. Eventually it would of course be desirable to understand these bounds analytically in the scenario with n = 3 RH neutrinos, which has partly been achieved for n = 2.