A Frequentist analysis of three right-handed neutrinos with GAMBIT

The extension of the Standard Model by right-handed neutrinos can not only explain the active neutrino masses via the seesaw mechanism, it is also able solve a number of long standing problems in cosmology. Especially, masses below the TeV scale are of particular interest as they can lead to a plethora of signatures in experimental searches. We present the first full frequentist analysis of the extension of the Standard Model by three right-handed neutrinos, with masses between 60 MeV and 500 GeV, using the Global and Modular BSM (beyond the Standard Model) Inference Tool GAMBIT. Our analysis is based on the Casas-Ibarra parametrisation and includes a large range of experimental constraints: active neutrino mixing, indirect constraints from, e.g., electroweak precision observables and lepton universality, and numerous direct searches for right-handed neutrinos. To study their overall effect, we derive combined profile likelihood results for the phenomenologically most relevant parameter projections. Furthermore, we discuss the role of (marginally) statistically preferred regions in the parameter space. Finally, we explore the flavour mixing pattern of the three right-handed neutrinos for different values of the lightest neutrino mass. Our results comprise the most comprehensive assessment of the model with three right-handed neutrinos model below the TeV scale so far, and provide a robust ground for exploring the impact of future constraints or detections.


Introduction
The observation of neutrino flavour oscillations is one of the strongest hints for the existence of particle physics beyond the Standard Model (SM). The oscillations imply that neutrinos have small masses, while the minimal SM predicts that they are massless. At the same time neutrinos are the only elementary fermions that are only known to exist with left handed chirality ν L . If right handed neutrinos ν R exist, one could immediately add a Dirac mass termν L M D ν R to the SM Lagrangian in analogy to all other known fermions. The fact that the ν R have not been seen yet could easily be explained because they are "sterile", i.e., not charged under any known gauge interactions. The same property also makes it possible for them to have a Majorana mass term ν R M M ν c R in addition to the Dirac mass. For eigenvalues of M M that are much larger than the observed light neutrino masses, the smallness of the neutrino masses can be explained via the seesaw mechanism [1][2][3]. Neutrino oscillation data is, however, not sufficient to pin down the value of M M , known as seesaw scale, because it is primarily sensitive to the combination M D M −1 M M T D . The range of allowed values spans from a few eV [4] up to the scale of Grand Unification [5]. For specific choices of their Majorana mass the ν R could in addition solve a number of long standing problems in cosmology. For instance, they could explain the baryon asymmetry of our Universe via leptogenesis during the decay [6] or production [7,8] of the heavy neutrinos or provide a viable dark matter candidate [9,10]. An overview of the cosmological implications of different choices of M M can e.g. be found in Ref. [11].
Experiments can directly search for heavy neutrinos if M M is below the TeV scale. Such searches have been performed in various different facilities, including high energy colliders and fixed target experiments. This is the mass range we consider in the present article. In addition, the ν R would indirectly affect precision observables or searches for rare processes. A summary of different existing constraints can be found in the reviews [11][12][13][14][15]. For the future a wide range of different searches have been proposed, an overview can be found in Refs. [16][17][18]. In order to decide about the best possible search strategy is it important to understand which parameter region is already ruled out by past experiments. This is in fact a non-trivial question because different observables are correlated in the seesaw model, and the requirement to simultaneously respect all known experimental results imposes stronger constrains on the model parameter space than superimposing individual bounds. Such global constraints can only be derived within a given model. An important quantity in this context is the unknown number n of right handed neutrino flavours. The minimal number that is required to explain the light neutrino oscillation data is n = 2, which would necessarily require the lightest SM neutrino to be massless. The minimal number that is required to generate masses for all three SM neutrinos is n = 3. This choice is also somewhat appealing in view of the fact that there are three fermion generations in the SM, and it is mandatory for anomaly freedom in many gauge extensions of the SM. The goal of the present work is to impose global constraints on the parameter space of the model with n = 3, based on the combination of direct, indirect and cosmological constraints summarised in section 3.
Several authors have previously imposed global constraints on the properties of right handed neutrinos. Here we exclusively focus on models in which the right handed neutrinos can explain the light neutrino oscillation data. 1 This e.g. excludes most sterile neutrino Dark Matter models because the feeble coupling of such particles that is required to ensure their longevity implies that its contribution to the light neutrino mass generation can be neglected [20]. 2 One of the most complete studies of indirect constraints on the parameter space for n = 2 in the last few years was presented in Ref. [23], where multiple electroweak precision observables and flavour-violating decays were included, along with tests of lepton universality and the unitarity of the CKM matrix. Loop corrections to some of these relations were considered in Ref. [24]. The authors of [25] included direct search constraints and those from big bang nucleosynthesis. The model with n = 3 is much less studied. Recent analyses of indirect constraints include Refs. [26,27], direct search constraints and BBN have been added to this in Ref. [28].
In this paper, we present the first full frequentist analysis of the n = 3 RHN extension of the SM, for a wide range of RHN masses from about 60 MeV to 500 GeV. We opted for a frequentist analysis rather than a Bayesian analysis since this is best suited to fully explore the valid parameter space while avoiding prior 1 The authors of Ref. [19] considered a single heavy neutrino, but made the conservative assumption that this particle may predominantly decay into a dark sector via new interactions. 2 We refer the reader to Refs. [21,22] for recent reviews on sterile neutrino Dark Matter. dependence and volume effects of the parameter space. We improve on different aspects of earlier analyses by combining all the strongest limits exerted by experiments as well as indirect signatures in a statistically consistent manner. Previous studies that examined the parameter space for n = 3 either used a subset of the constraints included here [26,27] or used less rigorous statistical methods [28] and focused on specific regions of the parameter space [29].
-While most previous studies fixed the mixing angles and mass differences in the active neutrino sector to the best fit values as presented in [30], we take into account likelihoods for the active neutrino observables. -Electroweak observables require precise calculations for its comparison with the extremely accurate measurements. We therefore use the calculation of the SM prediction for sin θ ef f w up to two-loop order [31]. -Most studies of lepton flavour violation in neutrino models focus exclusively on the most constraining processes, such as µ → eγ and µ → eee [23,28].
In this work we include all lepton flavour violating processes, in particular all leptonic τ decays, for which we use the most recent average of experimental results provided by HFLAV [32], as well as µ − e conversion in two nuclei (Pb and Ti). -For neutrinoless double-beta decay, in comparison with [28], we opt to carry out our analysis conservatively; in addition, the upper limit on the effective Majorana mass and hence the mixing is encoded in the form of a (one-sided) Gaussian likelihood, not as a strict cut. -Lepton universality tests are often centered on leptonic decays of mesons, K and π, τ -leptons and W -bosons [23]. We supplement these tests of universality with the recently observed semileptonic decays of B-mesons [33][34][35]. -We improve the treatment of CKM unitarity with respect to the discussion in Ref. [28]. -Concerning direct searches, previous studies have used only a subset of the experiments considered here [19,36], or chose to place a hard cut at the upper limits presented in the individual papers [25,28]. We implement the strongest constraints over the mass range as likelihoods. The statistical combination of these likelihoods also leads to more accurate profile likelihood contours in comparison to simply overlaying individual limits. -We study in detail the flavour mixing pattern of the three RHN, for different values of the lightest neutrino mass. We discuss the limit where the lightest neutrino is massless and the connection to the n = 2 case.
We use here the open-source software package GAM-BIT [37]. It includes an interface to Diver [38], a differential evolution-based scanner that provides efficient sampling performance for frequentist scans. This paper is organised as follows. In section 2, the model, parametrisation used and essential quantities are defined. All the observables and experiments that are considered are subsequently discussed in detail in section 3. Our scanning strategy, parameter ranges and applied priors are mentioned in section 4. The results are presented in section 5 and we discuss the implications of the combined constraints for future searches in section 6. In Appendix A we comment on the details of the implementation in GAMBIT, in Appendix B we explicitly give the expressions for the different observables, in Appendix C we provide details on how we interpret our results in view of the criterion of technical naturalness, and in Appendix D we show the different partial likelihoods.

Basic definitions
The addition of three RHNs to the particle content of the Standard Model introduces in total 18 new parameters. In this section we summarise basic relations in the seesaw model and define our notation, following Ref. [28].
The most general renormalisable Lagrangian that can be constructed from SM fields and the ν R has the following form: Hereby, L = (ν L , e L ) T indicate the left-handed leptons 3 of the SM and Φ is the Higgs doublet withΦ = Φ * and being the Levi-Civita tensor. M M is the Majorana mass matrix for ν R and F is the Yukawa coupling matrix. We work in a flavour basis where M M = diag(M 1 , M 2 , M 3 ). After electroweak symmetry breaking (EWSB), the complete neutrino mass term reads with where M D = F v, v being the Higgs vacuum expectation value (v = 174 GeV in the ground state). We include the one loop corrections δm 1loop ν and δM 1loop N as we aim for performing an analysis to be consistent at second order in the Yukawa couplings F . The mass matrix (3) can be diagonalised by a matrix of the form [24] U = cos(θ) sin(θ) − sin(θ † ) cos(θ † ) with Hereby, θ indicates the matrix that mediates the mixing between the active neutrinos ν L and the sterile neutrinos ν R . We can generally write with The additional complex conjugation of U N ensures that the relation among mass and flavour eigenstates will be analogous for LHNs and RHNs within the notation.
In the second relation in eq. (8) we have neglected the difference between the eigenvalues of M M and M N , which is of second order in θ. This is justified for the present purpose because of the experimental constraints on the magnitude of the elements θ αI , which we discuss further below.

The seesaw limit
The limit of small θ αI is usually referred to as the seesaw limit, it corresponds to M D M M (in terms of eigenvalues). It allows the approximation and leading to with m tree The loop correction to the light neutrino mixing matrix is given by [39]: where l(M I ) is a loop function given by The light and heavy neutrino mass eigenstates are described by the flavour vectors and respectively. We can define the matrices V ν and V N that represent the mixing between mass and interaction eigenstates in the respective sectors as while mixing between the two sectors is encoded in the matrix This quantity is of primary interest because it controls the interactions of the heavy neutrinos with the physical Higgs field h and the gauge bosons W and Z, Here g is the weak gauge coupling constant and θ W the Weinberg angle. For convenience, we introduce the notation From the relations (3) and (7) it is straightforward to derive the relation  [40]. This means that the production cross sections for heavy neutrinos are not affected. However, the branching ratio between lepton number violating (LNV) and lepton number conserving heavy neutrino decays is affected by U N [41]. This has no big effect on our scan because constraints from searches for LNV are sub-dominant in almost the entire mass range that we consider, but it may have important implications for future searches.
U N can have a big impact on the individual mixings U 2 αI of each heavy neutrino if all three Majorana masses are degenerate. This can be accommodated in technically natural scenarios discussed in the following section 2.5, cf. in particular footnote 6. The practical impact on experimental searches is, however, limited because most experiments are not able to kinematically resolve small mass splittings and therefore only probe U 2 α in this regime (rather than the couplings U 2 αI of individual heavy neutrino flavours). Also in this case observables that are sensitive to LNV are the only ones that are likely to be affected.
Finally, if the degeneracy between the heavy neutrino masses is accidental, then the proof in Ref. [40] does not apply, and U N can have a significant effect on the U 2 αI even if only two heavy neutrinos have degenerate masses. Our results contain a significant number of points of this kind because we performed several scans 4 Note that the approximation U N = I also allows to neglect δM 1loop N because it only amounts to a change in the matrix U N [40].
with "agnostic" parameter ranges that do not suppress fine-tuned points, cf. table 5. However, the fact that experiments are unlikely to resolve the individual resonances in this case implies that they are only sensitive to the quantities U 2 a , where the summation is to be taken over the mass degenerate heavy neutrino flavours only. As in the previous two cases, the effect of U N on the total production rate is minor because the matrix mainly re-distributes coupling between the mass degenerate states. The main affect would again be on LNV observables.
In summary, if any heavy neutrinos are discovered in the future, a comparison between the branching ratios of lepton number violating and lepton number conserving decays will give important insight into the mechanism of neutrino mass generation and will be crucial to identify any underlying symmetries.

Casas-Ibarra parametrisation
In the current work, we use the Casas-Ibarra parametrisation [42], generalised to include the 1-loop correction to the left-handed neutrino mass matrix [43]. This provides a simple way to impose constraints from light neutrino oscillation data in our scan. This parametrisation is based on the observation that m ν in eq. (12) can be expressed as with Since the loop function is smooth we can neglect the difference in the eigenvalues of M M and M N , In this scheme the sterile neutrino mixing matrix, i.e. the matrix encoding the mixing among LHNs and RHNs (20) can be written as where U ν is the PMNS matrix introduced above, m diag ν is the diagonalised, one-loop-corrected LHN mass matrix andM diag is the analogous RHN mass matrix, given by (28). Furthermore, R is a complex, orthogonal matrix that is parametrised by complex angles ω ij where R ij has the non-zero elements Since we work in the flavour basis in which the Yukawa couplings of the charged leptons are diagonal, U ν can be parametrised as where U ±δ = diag(e ∓iδ/2 , 1, e ±iδ/2 ) and V ij , parametrised by the LHN mixing angles θ ij , has non-zero elements analogous to R. Furthermore, α 1 , α 2 and δ are CP-violating phases.
The C-I parametrisation scheme generates by construction Yukawa couplings and mixing angles Θ that are consistent with light neutrino oscillation data up to second order in θ. This has two disadvantages. First, one may find it unsatisfactory that we treat light neutrino oscillation data differently from other constraints. Second, the C-I is a "bottom up" parametrisation. There is usually no simple relation between the C-I parameters and parameters that may be well-motivated from a model building viewpoint, and any theory-motivated prior on the RHNs' mixings and masses would acquire a rather convoluted form in the C-I parametrisation. In particular, there is no simple way to distinguish "natural" from "fine tuned" parameter choices. Hence, we refrain from performing Bayesian scans in the current work, and instead concentrate on a likelihood-based frequentist treatment. In view of the high dimensionality of the parameter space and the complicated functional form of the different constraints, the disadvantages of the C-I parametrisation are, however, compensated for by the numerical advantage that one gains.

The symmetry protected scenario
The smallness of the light neutrino masses m i can be explained in different ways by the seesaw relation (26). One possibility is that the N I are very heavy, i.e., M I v, in which case the smallness of m i is due to the smallness of the ratio v/M I . This choice for the mass scale(s) M I is well-motivated by Grand Unified Theories, 5 but raises the question of radiative corrections to the Higgs potential from the Yukawa couplings of the RHNs [46].
This "hierarchy problem" can be avoided in low scale seesaw scenarios. Low values of M I are natural because in the limit M I → 0 the B − L symmetry in the SM is restored. In this case, however, the smallness of m i can no longer be explained efficiently by the suppression of v/M I , as it typically requires couplings that are very small, in particular for seesaw scales as low as 100 MeV.
Such small values for fundamental parameters are considered 'unnatural' by many theorists [47], though some possible explanations have been proposed [48]. However, this estimate relies on the underlying assumption that there are no cancellations (accidental or otherwise) in the seesaw relation (26), which would allow for much larger U 2 αI = |Θ αI | 2 than the naive estimate (35) suggests while keeping the eigenvalues Hence, a technically natural [49] way to obtain small neutrino masses m i can be realised if the Lagrangian (1) approximately respects a B −L symmetry [50,51], whereL is a generalised lepton number under which combinations of the ν Ri are charged. Such B −L symmetry is exact if the Yukawa coupling and mass matrix take the form [52] (36) in which case the light neutrinos are exactly massless m i = 0. In order to generate non-zero light neutrino masses this symmetry has to be slightly broken, i.e., where the entries of the matrices µ and are small symmetry breaking parameters. If the symmetry is not exact M M can have offdiagonal elements, see for example Ref. [53]. Throughout this work we use a basis in which M M is diagonal. The diagonalisation affects the form of the Yukawa matrix F , but as long as the off diagonal elements of µ are small, this only leads to a small modification of the flavour structure. For the following discussion we will therefore adapt the simpler form [54] 6 6 An important exception is the case µ 1,M M . In that situation even small off-diagonal elements µ ij can lead to a comparably large misalignment between the basis in which F has the form (38) and the heavy neutrino mass basis, which means which that all heavy neutrinos have unsuppressed Yukawa couplings ∼ Fa in spite of the fact that a 1, cf. ref. [41] for a discussion. However, in this case all three mass eigenstate N I have approximately the same massM and cannot be distinguished kinematically. In this case the experimentally relevant mixing is U 2 a , the magnitude of which is controlled by the large entries Fa. Heavy neutrino oscillations in the detector [55][56][57][58][59][60][61][62][63][64] could provide an indirect way to access the small mass splitting and phenomenologically study this specific case.
with α , α , µ, 1 being small symmetry breaking parameters and F α being of the order of one. This means that one heavy neutrino practically decouples while the other two approximately form a Dirac spinor with mass M .
In this symmetry protected scenario there is no upper limit on U 2 αI from neutrino oscillation data. In the mass range considered here the upper limit comes from the experimental constraints, while for larger masses there is a theoretical bound U 2 αI < 4π(n − 1)(v/M ) 2 from the requirement that the Yukawa couplings remain perturbative [5]. This provides a theoretical motivation for a low scale seesaw with experimentally accessible mixings U 2 αI . Specific examples that motivate this limit include "inverse seesaw" [65][66][67][68], "linear seesaw" [69,70], scale invariant [71] and some technicolour-type models [72,73] and also the νMSM [8,50].

Connection to the model with n = 2
The parametrisation (38) suggests that the B −L symmetric limit for the model with n = 3 should contain the model with n = 2, as the third heavy neutrino decouples for a → 0. This is, for example, observed in the νMSM. However, some care is required when taking this limit if one wants to be consistent with neutrino oscillation data.
First, it is clear that not all seven symmetry breaking parameters a , a , µ can be set to zero because this would give exactly massless light neutrinos. Which of these parameters are non-zero and how small they are with respect to each other depends on the way how the symmetry is broken and thus on the particle physics model in which the Lagrangian (1) is embedded. It is not possible to make a model independent statement about the relative size of the a in relation to other model parameters.
Second, the parametrisation (38) is not the most general one: If we allow for small off diagonal elements in the general form (37), then all three heavy neutrinos can have unsuppressed interactions ifM M , cf. footnote 6. Hence, ifM is degenerate withM , one cannot expect to recover the n = 2 model even if α 1. Finally, as discussed in more detail in Appendix C, there are Casas-Ibarra parameter choices that yield small values of m ν0 , but correspond to highly fine-tuned scenarios where the this smallness is due to accidental cancellations. These solutions can imitate the symmetry protected scenario and can also circumvent the seesaw upper limit and thus reach high values of U 2 αI .

Observables, experiments and likelihoods
Models with heavy right-handed neutrinos, as described above, will alter the SM predictions for different observables that are already significantly constrained by experimental results. In this analysis, we implemented all relevant constraints such as active neutrino likelihoods ( In this section, we will focus on the physics and statistics aspects of our likelihood functions. The corresponding implementation of GAMBIT capabilities and module functions associated with the various observables are discussed in detail in Appendix A.

Electroweak precision observables
The leptonic charge currents are modified by the RHNs, and hence the value of G µ that is measured via the muon decay will differ from the actual Fermi constant G F which is defined in terms of the fine structure constant and mass of the Z boson. The correction can be written as [28] and is caused by the non-unitarity of the flavour mixing matrix V ν , see Eq. (18), which leads to a slight suppression of the muon decay. Both the weak mixing angle θ w and the mass of the W boson m W depend on G µ at one loop, which means they also get a correction from the active-sterile mixing matrix Θ, which is given by [23] where s 2 w = sin 2 θ w . Since experiments typically measure the effective Weinberg angle s 2 ef f , and assuming the QCD corrections factorize from the leptonic corrections [95], we use for the SM prediction the highly accurate calculation, including corrections up to twoloops, from [31] [s 2 Other electroweak precision observables affected by the presence of the heavy neutrinos are the decays of the Z and W bosons, in particular the invisible decay width of the Z boson, Γ inv , and the leptonic decays of W . Under the assumption that the radiative corrections factorize from the heavy neutrino contribution, at least up to order θ 2 [24,95], one can write the invisible decay width of the Z as [96] where we have neglected the contribution from Z → N i N j due to being of order θ 4 , and for the SM decay Z → ν i ν j we use the 2-loop calculation from [97]. The contribution of heavy neutrinos to the W decay widths to leptons can be written as [23] where we defined x α ≡ m 2 lα /m 2 W . We construct Gaussian likelihoods for these observables using the experimental measurements and uncertainties displayed in Table 1. All these observables depend on G µ (eq.(39)) either directly or through another observable (s w or m W ). Since the experimental measurements of these quantities are independent of each other, we assume them to be uncorrelated.

Lepton flavour violation
Flavour changing neutral processes, such as lepton flavour violation (LFV), are strongly suppressed in the Standard Model at one loop due to the GIM mechanism [99]. Hence, any non-trivial contribution to these processes from physics beyond the Standard Model would dominate over the SM contribution, which in turn makes the experimental determination of these observables a smoking gun of new physics. Several experiments have attempted to measure LFV processes with outstanding precision and they have imposed a set of upper limits on their branching fractions. In Table 2 we list the most significant of these observables, along  with the experiments that provided that bound. When more than one experiment is cited, the HFLAV average is used [32]. All upper bounds are given at the 90% C.L..
with the experimental upper bound on their branching ratios and the experiment that provided it. The experimental upper bounds for LFV µ and τ decays in Table 2 are given as branching fractions with respect to the total decay width of the respective lepton [98,110], In the model with three heavy neutrinos the leading contributions to these observables arise from dipole and box diagrams with mixing between the active and sterile neutrinos, given by the active-sterile mixing matrix Θ. The relevant LFV processes containing these diagrams are of the form The associated decay widths can be found in Appendix B.1.
Lastly, LFV processes can result in a neutrinoless µ − e conversion inside a nucleus. Muons captured by a nucleus typically decay in orbit providing a continuous spectrum of energy for the electron in the final state. In coherent flavour violating conversion, µ − N → e − N , final state electrons have a discrete energy spectrum, corresponding to the mass of the decaying muon. Consequently experiments measure the rate at which this conversion happens, with respect to the rate of capture by the nucleus, The corresponding expressions for the conversion ratio, as well as the nuclear parameters for the two nuclei studied, Ti 48 22 and Pb 208 82 , can be found in Appendix B.1.
The likelihoods for these LFV observables are all Gaussian upper limit likelihoods. They are computed as using the experimental data from Table 2. More specifically, we assume a measured value of x 0 for all observables 7 , and set σ = v/1.64 for full Gaussians and σ = v/1.28 for one-sided Gaussians, where v is the quoted upper 90% C.L. limit.

Lepton universality
Recent measurements of meson decays [33][34][35] have put into question the flavour-independence of leptonic charged currents, as predicted by the SM. Previous tests of lepton universality performed by LEP and SLC, using lifetime measurements of the tau and muon as well as the partial decay widths of the Z boson, showed no such deviation. This has lead to the formulation of many BSM theories attempting to explain the deviation shown in meson decays [96,111]. The presence of right-handed neutrinos modifies the leptonic currents and thus triggers a contribution to processes testing lepton universality such as in the fully leptonic decays of charged mesons, X + → l + ν, or the semileptonic decays of B mesons B 0/± → X 0/± l + l − .
In order to cancel the considerable hadronic uncertainties present in the decays of pseudoscalar mesons, lepton universality tests are best formulated using ratios between lepton species. For fully leptonic and semileptonic decays of mesons, these ratios are expressed as respectively. In case of fully leptonic decays, one can express the test of lepton universality in terms of deviations from the SM prediction as where the sterile neutrino contribution can be calculated from the active-sterile mixing matrix Θ as [28,112] where we used with ϑ being the Heaviside step function, r α ≡ m 2 lα /m 2 X and r I ≡ M 2 I /m 2 X . The SM predictions used in eq. (49) for the tests of lepton universality for pions and kaons are R π eµ,SM = 1.2354×10 −4 and R K eµ,SM = 2.477×10 −5 , respectively [113].
The contribution from heavy right-handed neutrinos to the semileptonic decays of B mesons is much less significant than to the leptonic decays. As argued in Ref. [96], the effect on B decays to charmed mesons, B ± → Dlν, is completely negligible. Semileptonic decays to K mesons are more affected, particularly the decays B + → K + l + l − and B 0 → K * 0 l + l − . Assuming that m l m K ( * ) and that the Wilson coefficient C 7 C 9 , C 10 , one can approximate the ratios R K and R K * as [114] and the BSM contributions to the Wilson coefficients ∆C α 9 and ∆C α 10 can be expressed as [115] ∆C α NNL calculations for the Standard Model contribution to the Wilson coefficients C 9 and C 10 used in Eq. (52) gives C SM 9 = 4.211 and C SM 10 = −4.103 [116,117]. In addition to meson decays, other common tests of lepton universality include the decays of the W boson to leptons as well as τ decays. The ratio of decay widths of W to charged leptons l α and l β can be written as [23] Deviations from the SM for the lepton universality test in τ decays follow the same form as in Eq. (50) and the SM prediction is R τ µe,SM = 0.973 [118].  These tests of lepton universality are implemented as Gaussian likelihoods centered on the experimentally measured value. The SM contributions to these observables, when relevant, as well as the experimental measurements, with their corresponding uncertainties 8 , are shown in Table 3, where two experimental measurements are shown for R K * corresponding to two regions of the dilepton invariant mass 0.045 < q 2 < 1.1(GeV 2 /c 4 ) for (1) and 1.1 < q 2 < 6.0(GeV 2 /c 4 ) for (2).

CKM unitarity
The determination of the CKM matrix elements (V exp CKM ) i ab is usually done under the implicit assumption of a zero active-sterile mixing matrix, Θ = 0. The measurements of the (V exp CKM ) i ab therefore need to be adjusted to take into account effects of RHNs.
Firstly, the smallest element of the CKM matrix, (V CKM ) ub , can be neglected in our study as its absolute value |(V CKM ) ub | 2 ∼ 10 −5 is much smaller than our sensitivity to the Θ parameter. Hence, under the assumption of the unitary of the CKM matrix, one can derive the following relation: 8 The experimental uncertainties for R B K ( * ) are obtained as the sum in quadrature of the statistical and systematic uncertainties provided by [33] and [35].
Thus, we use the various experimental measurements of (V exp CKM ) us [123][124][125] and (V exp CKM ) ud [126] to simultaneously constrain the true value of |(V CKM ) us | and active-sterile mixing matrix Θ.
Following Refs. [23,28], the experimental measurements and true value of CKM matrix element (V CKM ) us,ud are related via where we defined the functions f i to encode the contribution of RHNs to the process considered in each experiment. The decay processes considered to extract the value of |(V exp CKM ) us |, and the f (Θ) functions, are given by [23] The situation is simpler in the determination of the |(V exp CKM ) ud | element as the uncertainty is dominated by the superallowed 0 + → 0 + nuclear beta transitions measurements, which need to be modified accordingly to: The experimentally measured values of |(V exp CKM ) i us | in each of the decay processes above are listed in Tab. 4, and the value of |(V exp CKM ) ud | = 0.97417 ± 0.00021 is taken from the world average [126].
We thus construct the likelihood for this constraint from a chi-squared function, 2 ln L = −χ 2 , where Parameter Process Value Ref. Table 4: Experimental values of (V CKM )us and the average value of (V CKM ) ud used in the calculation of the CKM likelihood. The factor f+(0) = 0.959 ± 0.005 is taken from [127].
the discriminant measures the deviation of the true value (V CKM ) us,ud and the experimental measurements (V exp CKM ) i us,ud , and is given by Due to the unitarity relation in Eq. 56, the value (V CKM ) ud is obtained from (V CKM ) us for every parameter point, and thus the only free floating parameters are the value of (V CKM ) us and the active-sterile mixing matrix, Θ. For simplicity, and since this is the only constraint to depend strongly on the value of (V CKM ) us , we optimise on its value for each Θ, which removes the necessity of making (V CKM ) us part of the scanning model. This approach is similar to the discussion in [28], but we improve upon it by optimising on the true value (V CKM ) us , including the Θ corrections, for each parameter point, rather than the value measured experimentally.

Neutrinoless double-beta decay
Double-beta decay refers to the decay of two neutrons into two protons while emitting two electrons and two anti-neutrinos. In case of neutrinos having a Majorana nature, lepton number would be violated and neutrinoless double-beta decay (0νββ) induced. Besides the exchange the light neutrinos, the exchange of RHNs is similarly possible and would alter the expected effective neutrino mass m ββ . The effective mass is constrained by half life measurements of 0νββ decay. The most stringent limits are currently set by the GERDA experiment (Germanium) [128] with m ββ < 0.15−0.33 eV (90% CL), and KamLAND-Zen (Xenon) [129], m ββ < 0.061 − 0.165 eV (90% CL). The effective mass m ββ , can be theoretically evaluated in term of the mixings and masses of the light and right handed neutrinos [130] Hereby, the first term denotes the contribution from LHNs, the second the one from RHNs. With a typical momentum exchange of around 100 MeV in 0νββ decay, RHNs with a mass above this threshold participate in the process only virtually. This suppression is taken into account by the factor [130] The typical momentum exchange p 2 depends not only on the specific isotope in consideration but is also subject to the theoretical model in which the constraints are derived and the value of the nucleon axial-vector constant. An overview is given in [131]: For our analysis, we use the "Argonne" model and the lower of the two values for p 2 (quenched), which yields the most conservative constraints: p 2 = 178 MeV for xenon, and p 2 = 159 MeV for germanium. A more dedicated analysis of the impact of different limits due to nuclear uncertainties is beyond the scope of this work. Since we are focusing on profile likelihood for our results, this approach is largely equivalent to profiling over systematic uncertainties assuming a flat prior that spans the entire range of values p 2 in ref. [131]. For our analysis we use the experimental values, as stated above, as one-sided Gaussian likelihoods, choosing the higher of the two values in order to remain conservative.

Big Bang Nucleosynthesis
If RHNs decay shortly before or during BBN, the typical energy of decay products, here ∼ M I ≥ 50 MeV, is significantly higher than the plasma temperature at that time, ∼ 100 keV. Therefore, either by dissociating formed nuclei, or by causing deviations from thermal equilibrium, they will affect the abundances of primordial elements, which are however observationally well constrained. The requirement that the RHN decay happens sufficiently early enough before BBN implies an upper limit on the lifetime (τ I ) of RHNs, or equivalently, a lower bound on the mixing U 2 I [132]. However, in the presence of multiple RHN species, BBN cannot constrain individual mixing angles U 2 αI (22) but only the total mixing U 2 I (23). We consider decay channels of RHNs until the mass range of B mesons; beyond this mass, the constraint of BBN is considerably weaker. The expressions for the decays are taken from Ref. [133] and Ref. [134], as well as the values of the necessary decay constants. The decay width for each topology is listed in Appendix B.2.
In the current study, we require the lifetime of each RHN to be less than 0.1s [135], which is implemented as a hard cutoff in the parameter space. In principle, this limit can be weakened if the lightest active neutrino has a mass < O(10 −3 ) eV, since the RHNs do not necessarily thermalize in this case [136]. We leave, however, the implementation of refined BBN constraints in GAMBIT for future work 9 .

Direct RHN searches
Different experiments search with various approaches directly for RHNs. One can distinguish between three types: peak searches (PIENU), searches at beam dump experiments (PS-191, CHARM, E949, NuTeV), and searches at e + e − or pp colliders (DELPHI, ATLAS, CMS).
One possibility to look for RHN, is to search for peaks in the lepton energy spectrum of a meson decay. If, for example, a meson of mass m X decays into an RHN of mass M I and an electron/muon with mass m lα , this peak will be approximately at Even in situations where backgrounds are sizeable, a peak search can hence be used to impose constraints on the mixing. In beam dump experiments, the large background signal that is usually present near the target hinders the detection of charged particles that are produced along with the RHNs. On the other hand, RHNs with mass below the D meson scale can be long-lived enough to travel macroscopic distances. Looking for their charged decay products some distance away from the target leads to (almost) background-free experimental situations.
In collision experiments (e + e − or pp), vector bosons or mesons get produced that subsequently can decay leptonically. The bounds on these processes are then able to constrain the corresponding active-sterile mixing angles in a certain mass range.
To implement the direct detection constraints as likelihoods, we follow two different approaches, depending on the information that is provided in each study. Firstly, some of the experiments found no signal events and had no background counts after cuts (DELPHI, CHARM, PS191 and NuTeV). In this case, since the processes in the experiments are essentially Poissonian, we construct the likelihood (to observe n events) as a Poisson distribution. The number of expected counts, µ, is a function of the RHN masses and mixings, i.e. µ = µ(M I , U 4 αI ) (assuming the experiment does so as well, the fourth power takes both production and decay of RHNs into account). For expected µ events and background b, the likelihood is: With no reported detections (n = 0) and background cuts reducing b to approximately zero, To connect µ with our model parameters, we use the fact that the expected signal counts are proportional to the LHN-RHN mixing, µ ∝ U 4 αI . The factor of proportionality is set to reproduce the results from the experimental papers (assuming that these limits are based on the common Feldman-Cousins procedure [137], where e.g. a 95% CL upper limit would correspond to an expected number of signal counts of µ = 3.09).
On the other hand, for the experiments which either quote non-zero signal events and/or backgrounds, or if this information is ambiguous (CHARM (ν τ reinterpretation), PIENU, ATLAS and E949), we model the constraint likelihood as Gaussian upper limits, i.e. we model them as half-Gaussians with zero mean and error set according to the confidence level at which the results are presented. For example, in the case of an experiment that presents limits at 90% CL, for a half Gaussian, this lies within 1.28σ of the mean.
It is worth noting that collider experiments often use simplified model assumptions to compute the confidence level intervals presented in their results. Since we use these to construct our likelihoods, we are incorporating these assumptions as well, in spite of the fact that our confidence intervals are computed by profiling over the multidimensional parameter space. Given that a full collider simulation is beyond the scope of this study, we employ the provided simplified model limits as given. We acknowledge, however, that the true limits may be slightly weaker due to, e.g a reduction of the production cross-section, and we defer the exploration of the differences between the collider predictions of simplified and full models to future work.

PIENU
The PIENU experiment [138] sought to detect RHNs in the mass range of 68 − 129 MeV by searching for peaks in the energy spectrum of the decay process π + → e + ν. It was, hence, sensitive to the mixing |Θ eI | 2 ≡ U 2 eI and µ in eq. (72) is also taken to scale as U 2 eI in our analysis. Although no peaks were found, exact information on the number of background events is unavailable. Further, production processes in peak searches are, in general, unaffected by the Majorana/Dirac nature of the RHNs; hence, no correction is necessary here.
The constraints on U 2 eI are at 90% CL, so it is implemented in GAMBIT as a half-Gaussian with zero mean and error set at 1.28σ.
After our analysis was complete we became aware of the slightly stronger updated constraints presented in Ref. [139], which are not included in our scan.

PS-191
This experiment [140] was designed for the purpose of detecting neutrino decays. RHNs would be produced via either of the following mechanisms: π + /K + → e + ν e , or π + /K + → µ + ν µ , and would then decay via Thus, PS-191 could constrain the quantities U 4 eI and U 4 µI for RHNs with a mass between 20 − 450 MeV.
Having found no signal or background events, it placed constraints on these quantities at 90% CL. We deviate from the original analysis in two ways. The first is necessitated by the fact that in the original analysis, the constraints were derived under the assumption that the RHNs interact only through the charged current. In [36], these limits were re-interpreted with the inclusion of neutral current interactions. Thus, instead of the signal count being proportional to the fourth power of the relevant flavour mixing, it is proportional to U 2 e/µI × α c α U 2 αI , with the coefficients given by We use these revised bounds here. The limits are encoded in likelihood form as in eqn. (72), with the aforementioned proportionality factor being 2.44.

CHARM
RHNs were searched for in CHARM [141] using two strategies, one with a neutrino beam from dumping protons on copper (BD) and another using a wide-band neutrino beam (WBB) from primary protons. In BD, the production of RHNs was assumed to occur through the decay of D mesons. They would then decay via µ + e − ν µ (and the anti-particle counterparts) and the decay products were looked for.
In WBB, RHN production was assumed to occur via neutrino-nucleus neutral current scattering ν µ N → ν R X. The subsequent decay ν R → µR, R representing hadrons, was then searched for. The limits from the WBB analysis are, however, weaker than those exerted by other experiments in the same mass range, and are not considered here.
The BD analysis yielded no candidate events or background and hence placed limits on U eI and U µI at 90% CL. Further, the original analysis assumed the possibility of RHNs interacting solely via the charged current; we use the results re-interpreted after the inclusion of neutral current interactions [36] as discussed in section 3.4.2, i.e. the signal count is proportional to U 2 e/µI × α c α U 2 αI and once again use eqn. (72) to represent the likelihood, with the proportionality factor being 2.44.
In [142], the data from the CHARM experiment was re-analyzed assuming that RHNs mix solely with tauflavoured leptons, and was able to place limits at 90% CL on U τ I , which we implement as a half-Gaussian with zero mean and error set at 1.28σ.
Dirac RHNs were assumed in both the original and tau-specific analyses, so the limits presented are also re-scaled by dividing them by √ 2.

E949
In this experiment [143,144], RHNs were searched for in the decay of kaons produced in a beam dump: K + → µ + ν R . Constraints on U µI were placed at 90% CL in the mass range 175 − 300 MeV; we also divide the limits by a factor of √ 2 to account for the Majorana nature of RHNs in our model.
The likelihood is modeled as a half-Gaussian with zero mean, error set at 1.28σ and µ ∝ U 2 µI .

NuTeV
The NuTeV experiment [145] searched for RHNs through their decay into the following final states: µeν, µµν, µπ and µρ. They were assumed to be produced in the decay of mesons. 90% CL limits on U µI were placed for RHNs with a mass between 0.25 − 2 GeV. Information about the assumed Dirac or Majorana nature of the RHNs is not present, so we take the conser-vative route and presume Majorana RHNs were considered in the analysis. No candidate events or background were detected, so the likelihood is modeled as in eq. (72), with a proportionality factor of 2.44 and µ scaling as U 4 αI .

DELPHI
At DELPHI [146], e + e − → Z 0 → ν Rν was the dominant RHN production mechanism; the process Z 0 → ν RνR would be suppressed due to the additional U 2 factor. The products of the RHN decaying via the weak and neutral current were then searched for, according to: DELPHI could constrain Θ eI , Θ µI and Θ τ I for RHNs having a mass between 0.5 − 80 GeV.
Since the RHNs could have existed long enough to travel macroscopic distances of upto 100 cm, different signatures had to be considered and the analysis was split to tackle the short-and long-lived cases separately.
In the short-lived RHN case, depending on the particle mass, two signatures were looked for. For masses less than about 30 GeV, due to the large boost received by the RHNs, the signature would be a monojet. Background coming from leptonic Z boson decays or γγ processes were accounted for. Higher masses open the decay channel into qq (and a lepton, depending on the channel), and the signature in this case would be two acollinear jets which are also acoplanar with respect to the beam axis. Most of the background in this scenario came from hadronic Z decays with missing energy; a neural network was used to remove all of them from the final data.
Longer-lived RHNs were looked for using displaced vertices and calorimeter clusters. The former was useful in tracking RHNs with an intermediate lifetime; however, a cluster finding algorithm along with vertex reconstruction did not find any signals. Calorimeter clusters were used to detect the longest-lived RHNs, whose decay products would interact with the outermost layers/components of the experimental setup: the signature would be a cluster of hits in a small angular region coincident with the beam collision, which could be traced back to the initial interaction point.
The analysis was carried out assuming Majorana RHNs and yielded one candidate event and no background events. In our analysis, this means the proportionality factor is 3.09 and µ scales as U 4 αI . A caveat must be mentioned here: the DELPHI analysis presented bounds on the mixing in a flavourindependent manner: the limit on U 2 , as presented in the paper, applies equally to U 2 e , U 2 µ and U 2 τ , as they mention. In the mass range under consideration, the mass of the tauon will, of course, influence the strength of the limit and, as they quote, the presented bounds become weaker for masses below ∼ 4 GeV. However, the extent of the kinematic suppression due to the tauon mass is not quantitatively discussed; we use the limits as is, noting that it is highly likely that NA62 will subsume these bounds in the near future [147].

ATLAS
The process relevant for RHN production in AT-LAS [148] is pp → (W ± ) * → l ± ν R . The RHNs were taken to be heavier than the W boson, allowing it to decay to a lepton a W boson: ν R → l ± W ∓ ; the W boson would then decay predominantly into a quarkantiquark pair, and the signature of this decay chain was searched for, with either two electrons or muons in the final state. 10 Hence, in our analysis, µ ∝ U 4 αI , α = e, µ. The original analysis was carried out under the assumption of Majorana RHNs, so no additional correction is necessary.
The analysis placed 95% CL limits on the two mixing angles in the mass range of 100−500 GeV. Details on the number of observed/expected events and background is available and could be cast into a likelihood function combining Poissonian and Gaussian errors; however, we find that implementing the limits in GAMBIT as a half-Gaussian with zero mean and error set at 1.64σ reproduces the experimental limits well enough for the purpose of a global fit.

CMS
With the LHC having run with a center-of-mass energy of 13 TeV, the CMS detector searched for different event signatures of the same process as ATLAS. 95% CL limits were calculated for U eI and U µI for RHNs with mass between 1 GeV and 1.2 TeV [149].
As before, Majorana RHNs were assumed in the analysis, and our implementation of the limits mirrors that of ATLAS.
Note that updated bounds from ATLAS [150] and CMS [151,152] have been released, but are not included, 10 There is an ongoing dispute in the literature on whether the rate of LNV processes at collider experiments are always suppressed by the small parameters i and µ in Eq. (38) and therefore unobservably small (roughly of the order of the "naive seesaw estimate") [51,52] or whether coherent flavour oscillations can lead to LNV signatures in spite of the smallness of these parameters [57,60,61]. In the range of M I below the electroweak scale under consideration here, the strongest direct search constraints do not come from experimental signatures that rely on LNV, and our results are therefore only mildly affected by the outcome of this discussion.
since these papers came out after our scans were completed. However, the new bounds from ATLAS are comparable to those from DELPHI, and the newer dilepton search from CMS only produces stronger bounds for RHN masses above ∼ 500 GeV, which is beyond our range of study.

LHCb
LHCb has performed direct searches for heavy neutrinos. The most recent results [153] were derived with an inconsistent model and have been corrected in Ref. [154]. They are subdominant in the mass range considered here. In Ref. [155] the results of a generic long lived particle search [156] has been re-interpreted in the context of heavy neutrinos. We do not include these results here because the conservative interpretation does not yield stronger bounds than the ones we include.

Scanning strategy and parameter ranges
In this work, we focus on the exploration of the RHN parameter space using frequentist statistics. Our main goal is to establish the ranges of RHN parameters that are not yet explored by experiments, and a frequentist approach delivers a suitable and prior-independent method. We are dealing with a high dimensional parameter space, which we have to project into two-dimensional plots. To this end, the central quantity of interest is the profile likelihood, which is, for fixed parameters of interest θ 1 and θ 2 , the maximum value of the (log-)likelihood function that can be obtained when maximizing over the remaining parameters η.
Our scanning strategy is designed in order to explore the complex parameter space of the RHN model such that we obtain reliable results for the projections shown in this work. To this end, we perform a large set of scans with different settings which we then merge into a single dataset. We study the normal and inverted hierarchy  independently, in order to avoid artificially favouring one over the other due to the different normalisation of the active neutrino likelihoods (c.f. Section 3.1). Hence, we make independent scans for each of the neutrino mass hierarchies, normal and inverted, for the full set of scans described below.

Parameters and priors
The parameter ranges and priors for the original scans can be seen in Table 5. We emphasize that 'priors' do here not correspond to priors in the Bayesian sense, but rather determine the efficiency with which different regions of the parameter space are explored. For convergent scans, the results are prior-independent.
We have chosen to split the complex angles ω ij into their real and imaginary parts. The active-sterile mixings depend strongly on the imaginary parts of ω ij and large values of Imω produce mixings that are too large to pass any constraints, so we take a conservative range Imω ∈ [− 15,15], and also pre-emptively disallow choices that lead to |Θ| 2 ij > 1. As discussed in 2.5, a condition for an approximate B −L symmetry to be realized is for two RHNs to have almost degenerate masses, which extends the range of the mixings so that they can be probed by experiments. This provides motivation for using a logarithmic prior on the RHN masses, also allowing the scanner to sample better the region close to the limits of the most constraining experiments/observables. The C-I parametrisation, as defined in Section 2.4, together with the particular parametrisation choice of R in eq. 30, was found to not fully cover the entire parameter space. To circumvent this and ensure that all possible couplings are covered by the scans, we introduce an additional parameter to the scan R order with discrete values [1,6] corresponding to each of the possible permutations of the definition of R in terms of R ij . This allows full coverage of the coupling space and, since the likelihood is conceptually independent of the order in R (and confirmed by the data), it ensures an uniform distribution of values in the parameter R order .
Out of the active neutrino parameters, only α 1 and α 2 are unconstrained by oscillation data, hence they are allowed to vary freely from 0 to 2π with flat priors. The ranges for the other neutrino phases and angles are taken from the 3σ ranges from the NuFit collaboration [74], also with flat priors. The mass of the lightest active neutrino, m ν0 , has a definite impact on the lower bound of U 2 I (23) [28], so we choose a logarithmic prior, which enables us to examine this impact in greater detail than a flat prior would allow and keeps the BBN limits relevant [136]. The upper limit on m ν0 is chosen as the broad cosmological bound given by Planck [94], m ν < 0.23 eV 11 . In order to better fit the active neutrino data, the mass splittings ∆m 2 21 and ∆m 2 3l are chosen as scan parameters, where l = 1 and ∆m 2 3l > 0 for normal hierarchy and l = 2 and ∆m 2 3l < 0 for inverted hierarchy.
Since the construction of the mixing matrix in the C-I parametrisation depends on m H (1-loop correction), as seen in 2.4, we take m H as a nuisance parameter with a Gaussian distribution around its averaged measured value [98] and a flat prior. Other SM parameters are fixed to their PDG values [98].

Targeted scans
We encountered a number of challenges while sampling the full RHN parameter space. One reason is connected to the behaviour of the likelihood function over the whole parameter range. The adopted scanning algorithm (Diver, see below for details) is designed to find regions of maximum likelihood across the parameter space. However, as we will discuss later when we study the effect of each individual observable, most constraints have flat contributions to the likelihood in a large portion of the parameter space. Hence, the scanner often does not fully explore large regions with equal or worse likelihood. This happens especially near the experimental bounds. Furthermore, although high couplings are possible between active and sterile neutrino sector, they often lie in the symmetry protected regime, as described in Section 2.5 and/or require severe fine-tuning of the parameters. Again, exploring these regions turned out to be challenging. Therefore, we designed and performed a large set of targeted scans to fully saturate the experimental bounds, the list of which can be found in Table 6. The design strategies we adopted for these targeted scans can be summarised as follows.
First, all targeted scans were performed using a differential RHN model, where the parameter M 2 is replaced by ∆M 21 , with a logarithmic prior. This allows the exploration of the symmetry protected region, with near degenerate masses for two right-handed neutrinos.
Most of the experimental bounds occur at high couplings, thus in order to encourage the scanner to explore the high coupling regions, we added an artificial likelihood to the scan to drive the scan to the unexplored boundaries. To saturate the experimental bounds for each coupling U 2 αI , α = e, µ, τ , different targeted scans were performed using this coupling slide likelihood on each of the couplings, of the form s log U 2 αI + m log M I . Table 6 shows the parameter that is optimised in each scan, α, and the coefficients, (s, m). This contribution was later removed from the data in the postprocessing stage.
The targeted scans were further split along the M I axis following the limits of the various experimental constraints (mostly from direct searches). This ensures that each coupling (with the selection above) saturates the most relevant experimental upper bound in each mass range. Additionally, some scans used different values of ∆M 21 and/or m ν0 to further force the scan into fine-tuned regions of parameter space. The ranges used for M I , ∆M 21 and m ν0 for each scan are specified in Table 6.
We found that some of the experimental likelihoods provide positive contributions to the total likelihood in specific regions of the parameter space. This forced the scan towards those regions, leaving others unexplored. Although this is a rather interesting feature, and will be discussed in detail later, it prevented a thorough exploration of the full parameter space. We thus chose to remove the likelihood contribution of R K eµ from the total likelihood that drives the scan, adding it later in postprocessing. Other likelihoods with positive contributions, Γ inv , CKM and R τ eµ , tended to force the scan towards large U 2 τ I couplings. Although desirable to saturate the limits, this also left regions with low τ coupling undersampled. Thus, a cut on the coupling U 2 τ I was enforced in some scans to fully sample all regions. The adopted strategy for scanning was driven by the need to fully sample the parameter space. The results from all the diverse scans were combined into a single dataset after some postprocessing (see below). This does not pose a problem for the statistical interpretation, since we are interested in the profile likelihood, which only becomes more accurately estimated when adding additional chains.

Scanning framework
To perform the detailed scans, we make use of the GAM-BIT framework, as described in Appendix A, and the differential evolution scanner Diver, version 1.0.4 [174], which is a self-adaptive sampler, capable of sampling the profile likelihood more efficiently than other scanners. We choose a population size of NP = 19200 and a convergence threshold of convthresh = 10 −10 . After some tests, we have concluded that the aggressive λjDE setting in Diver provides an improvement on the sampling of the parameter space, since it is more suited for sampling fine-tuned regions.
These scanner settings, including the very low convergence threshold, together with the scanning strategy described above, ensure a thorough exploration of the parameter space, albeit at the price of CPU time. Despite the fact that none of the observables used required heavy computation or simulations, most scans took between 2 and 10 hours of running time on a large number of supercomputer cores varying between 250 and 780. All tests and scans were carried out across several supercomputer facilities, including the MareNostrum supercluster in Barcelona, Marconi in Bologna, LISA/Surfsara through the University of Amsterdam and Prometheus in Krakow.

Data postprocessing
Upon completion of the scans, a number of postprocessing tasks were performed on the data to prepare it for plotting. As previously mentioned, the first of these tasks was to remove the artificial coupling slide likelihood used to drive the scans to high couplings.
Due to the large amount of scans performed and the low convergence threshold used, the size of the samples surpassed 1TB for each hierarchy, rendering them unmanageable for most plotting routines. We hence performed a few operations on the scan results prior     Table 5.
to combining them. With the target of showing profile likelihood plots in the M I vs U 2 αI planes, we hence extracted a subset of the data points optimised in these planes, with a resolution of 10 −5 . Since most scans were targeted to saturate the limits for a particular coupling (see Table 6) we perform this reduction of the data in the respective mass vs coupling two-dimensional planes. The combined set will hence be optimised for all couplings. Additionally, and independent reduction of the data is performed on the planes m ν0 vs U 2 αI , since we intent to study the effect of m ν0 cuts on the coupling limits.
The C-I parametrisation is inherently biased towards the different right-handed neutrinos, in particular through their masses M I . The relevant physics processes should be invariant to permutations of the three heavy neutrinos. Hence, to remove this bias, after combining the reduced datasets for all the scans, we conduct a symmetrization procedure over the combined datasets. We therefore symmetrize over M I as well as U αI , which will increase the size of the datasets six fold.
Lastly, in order to compare with the n = 2 case, two further datasets were obtained, for normal and inverted ordering, where the data points are required to lie in the symmetry protected region.
Out of the incalculable amount of data points we collected through our scanning procedures, a total of 37.7 million valid data samples were used for plotting. Of which 9.4M correspond to normal hierarchy and 8.7M for inverted hierarchy, optimised on M I vs U 2 αI planes, and 9.9M for normal and 9.7M for inverted hierarchy, optimised on m ν0 vs U 2 αI planes. The datasets with points in the symmetry protected region have over 71k and 20k valid data samples for normal and inverted hierarchy, respectively. These samples can be found in Zenodo [175].

Capped likelihood
The figures in this article show the so-called capped profile likelihood (unless stated otherwise), which is defined in each of the scanned point to an equal or worse fit than the SM: L = min[L SM , L RHN ]. It can thus be interpreted as exclusion-only likelihood. Capped likelihoods have been used in previous studies, particularly in the context of collider searches [176,177]. The rationale behind the use of this capped likelihood is the presence of positive (above SM) contributions to the log likelihood from various observables. Importantly, these 'excesses' would not show up as localized features in the total profile likelihood, as there is enough of freedom to add points in the M I − U 2 αI plane to find M J , J = I with values that would saturate the excess likelihood. Thus a very large fraction of the parameter points would have the maximum allowed likelihood from the combination of all excesses. This effect forces to separate the exclusion studies from the possible signal observation. Thus, in most of the paper, we use the capped likelihood to present parameter constraints. The excess likelihoods will be discussed separately in Sec. 5.4.

General constraints on the RHN mass and mixing
The constraints are shown for the couplings U 2 αI to the active neutrino flavours α = (e, µ, τ ), as well as their combination U 2 I = α U 2 αI , as functions of the heavy neutrino masses M I . Here, the second index can refer to any of the heavy neutrino flavours I = (1, 2, 3), because their labelling is not physical. The allowed profile likelihood regions are flat for most of the parameter space, in particular for small couplings U 2 αI , and drop smoothly at high couplings following the relevant upper limits. The white lines around the experimental limits mark the 1σ and 2σ contours.
The largest values of mixings U 2 αI for all flavours are allowed for M I above the masses of the weak gauge bosons. In this regime the direct searches at colliders are sub-dominant, and the heavy neutrino properties are primarily constrained from above due to electroweak precision observables, lepton flavour violation and CKM constraints. The upper limits on the couplings U 2 αI within 2σ of the highest likelihood for each hiearchy and flavour in the high mass region can be found in Table 7. It can be readily noticed that the upper limits for the τ couplings is much larger than for the other two flavours, which can be understood because the limits from EWPO and LFV are stronger for e and µ (see also Sec. 5.4).
For M I between the masses of the D mesons and the W boson the limits from direct searches dominate because the heavy neutrinos can be produced efficiently via the s-channel exchange of on-shell W bosons. In the range between the D meson masses and the W boson mass, the limits from the DELPHI [146] and CMS [149] experiments compete to impose the strongest bound.  Below the D meson mass the constraints on U 2 eI and U 2 µI are dominated by direct search constraints from fixed target experiments, in particular CHARM [141] and NuTeV [145] above the kaon mass, PS-191 [140] and E949 [144] between the pion and kaon mass and pion decay experiments at even lower masses. In this regime the global constraints on U 2 eI and U 2 µI are in good approximation given by the direct search constraints, as discussed in Sec. 5.2 and Figs. 5-7. This is in contrast to the model with n = 2, where the global fits rule out a significant mass range below the kaon mass that appears to be allowed if one simply superimposes the direct constraints in the mass-mixing planes [25]. For U 2 τ I , the direct search constraints are much weaker, the limit from long-lived particle searches by DELPHI remains the most significant one in our scans. Figure 3 shows that direct searches become subdominant for the τ coupling and the EWPO limit is saturated for a considerable range of masses below the kaon mass.
For masses below roughly 0.3 GeV the global constraints are stronger than the sum of their ingredients due to an interplay of the lower bound from BBN on the mixings, the upper bounds on U 2 eI and U 2 µI from direct searches and the constraints on the heavy neutrino flavour mixing pattern from neutrino oscillation data (discussed further below in Sec. 5.3). The latter disfavours large hierarchies amongst the couplings to individual SM flavours, though these constraints are weaker than in the model with n = 2 [25,178]. This implies that upper bounds on combinations of U 2 eI and U 2 µI indirectly constrain U 2 τ I . The BBN constraint on the lifetime does not impose a constraint on any individual coupling U 2 αI , but requires at least some of them to be sizeable and practically translates into a lower bound on U 2 I that is visible in Figure 4. Both, the BBN constraint and the constraint on the flavour mixing pattern (that will be discussed in more detail in Sec. 5.3 and is visible in Fig. 9) leads to the lower and upper bounds on U 2 τ I that are visible in Figure 3.
The upper bound on the total mixing U 2 I from the global constraints can roughly be identified with the bound on U τ I across the entire mass range as it is  constrained the weakest. The lower bound is again given by the lifetime constraint from BBN. In addition, there is a lower bound from the requirement to explain the light neutrino oscillation data that depends on m ν0 and is therefore only visible if one imposes a cut on this unknown quantity. Our results agree with the analytic estimates made in Ref. [40] and are illustrated in Figure  8. , with the coefficients c α from eq. (73). As we profile over the other two couplings, the strongest limit for the α flavour for these experiments would correspond to (U exp I ) 2 /c α , with (U exp I ) 2 being the reported limit by the experiment. Hence the former ratio is what is shown in the figures as the PS-191 and CHARM limits. As observed, U 2 eI is constrained by PS-191 and CHARM in the lowest and next-to-lowest mass regions, whereas they are superseded by the limits from E949 and NuTeV for U 2 µI . In the lowest mass region for the µ coupling it would appear that the E949 bound is in fact not saturated as the experimental limit falls below the data. This is however just an artifact of binning and interpolation in that region and the fact that the E949 limit is quite jagged. To illustrate this, we show in Figure 7 a zoom into the lowest mass region from Figure 6, where it can  be seen clearly that the profile likelihood follows the limits of E949.

Discussion of individual bounds
As mentioned above, EWPO (including Γ inv , m W , W -decays and s w ), LFV and CKM constraints become relevant for very large couplings and are thus the dominant limit in the high mass region, as well as a small region at small masses for the τ coupling ( Figure 3). Besides providing constraints, in particular Γ inv and CKM observables are also responsible for the slight excesses in the total likelihood, which we will discuss in Sec. 5.4. Other constraints included in the analysis have little to no effect on the profile likelihoods as shown above.
Among the leptonic decays, only R K eµ has some impact on the likelihood, with a negative contribution at masses below 0.45 GeV. Both R K eµ and R τ eµ show minor excesses in total likelihood, which again will be discussed later. Other lepton universality constraints have only little effect on the likelihood.

Neutrinoless double beta decay is often a relevant constraint and sets strong upper bounds on U 2
eI . This rules out large volumina in the 18 dimensional model parameter space. However, in the limit where lepton number is approximately preserved the expected signal from 0νββ is suppressed. Since many of our parameter points are in this symmetry protected scenario, particularly at high couplings, the impact of this constraint on the likelihoods in the projection on the mass-mixing plane is minimal.
The effect of BBN can be seen in the lower limits of Figure 3 and 4. The lower limit on U 2 I is a direct consequence of BBN, as lower couplings would mean that right-handed neutrinos would not decay before BBN  and thus affect the abundance of primordial elements. Although no individual limits are imposed by the BBN constraint on the couplings, the strong upper limit on the e and µ flavours at low masses has the side effect of setting a lower limits on U 2 τ I , as seen in Figure 3. For a more detailed explanation of the effect of each partial likelihood, and associated figures, we refer to Appendix D.

Lightest neutrino mass and flavour mixing
Oscillation data strongly constrains most of the active neutrino parameters, in particular the mass splittings ∆m 2 21 and ∆m 2 3l , the mixing angles θ ij and CP phase δ CP , whereas the lightest neutrino mass m ν0 remains unknown. There are upper bounds from cosmology on the sum i m i that depend on the active neutrino mass hierarchy, the data set used and the underlying cosmological model. The value quoted by the Planck collaboration for a standard cosmological model is i m i < 0.12 eV [179], a discussion of how this value changes with differ-ent assumptions can e.g. be found in the Particle Data Group Report [180]. In fact, using the best fit values for the mass splittings from the NuFit data [74] and the conservative value i m i > 0.23 eV, we can infer the upper limits of m ν0 < 7.12 × 10 −2 eV for normal and m ν0 < 6.55 × 10 −2 eV for inverted hierarchy.
The value of m ν0 strongly impacts on the lower limit on U 2 I . One can obtain a reliable estimate of the lower bounds on U 2 I by setting R = 1 [40]. This makes the PMNS matrix unitary, and the lower limit one the smallest mixing can be estimated as U 2 I m ν0 /M I . Using this approximation, we show in Figure 8 the lower limits on U 2 I that we obtain in our scans for different values of m ν0 = (0.05, 10 −2 , 10 −3 , 10 −4 ) eV. In the case of m ν0 = 0 there is no absolute lower limit on the coupling from the seesaw mechanism, and the residual lower limit on U 2 I is due to the BBN constraint.
The lightest neutrino mass has an important effect on the pattern of flavour mixing. In the limit of large m ν0 , there is almost no constraint on the allowed flavour ratios U 2 αI /U 2 I . This is shown in Figure 9 by the black   solid (dashed) contours, which indicate the allowed region within 1σ (2σ) where the lightest neutrino mass is m ν0 < 10 meV (close to the cosmological bound stated above). In this case, there is no visible upper limit on U 2 µI /U 2 I or U 2 τ I /U 2 I for normal hierarchy, whereas U 2 eI /U 2 I is constrained 0.95. Conversely, for inverted hierarchy there is an upper limit for the µ and τ flavours, but none for the e flavour. However, for smaller values of m ν0 , the allowed range for the flavour mixing pattern becomes significantly constrained. 12 This is shown 12 When constraining mν 0 to very small values, we almost decouple one right handed neutrinos. The contribution of this feebly coupled state to the generation of light neutrino masses is negligible, which in return implies that its properties are almost unconstrained by neutrino oscillation data, and such is its flavour mixing pattern. Thus extreme ratios U 2 αI /U 2 I can in principle occur for this particular heavy neutrino, although the absolute values of U 2 I remains negligible, and it has no effect on any near future experiment. Since our focus is primarily on heavy neutrinos that make a measurable contribution to by the lines for m ν0 < 1 meV (blue), m ν0 < 0.1 meV (green) and m ν0 < 0.01 meV (red). For masses lower than 0.01 meV the constraints saturate and the size of the ellipse remains almost constant. This can be also seen in Figure 10, where the largest coupling ratio is plotted for each flavour as function of neutrino mass.
It is instructive to compare our results to the constraints on the flavour mixing pattern in the scenario with n = 2 that were found in Refs. [25,178,182]. For this purpose it is not sufficient to simply insert very the generation of light neutrino masses and/or may be discovered in experiments, we applied a cut on M I U 2 I > 10 −10 GeV in Figure 9 and  small values for m ν0 in the parameterisation (29) because such values can also be achieved due to accidental cancellations in the light neutrino mass matrix (without decoupling of any of the heavy neutrinos), cf. Section 2.6. To remove such fine tuned points we impose the following cuts Here is an arbitrarily small number, which we choose as = 0.01 for convenience. In addition, we work in the limits as defined by |Imω 23 | 1 and Reω 13 ∼ π/2 for normal hierarchy, and |Imω 12 | 1 for inverted hierarchy (c.f. Appendix C). Note that we randomised the order of the matrices R ij , and hence for normal hierarchy we can only reproduce the true symmetry protected regime for the permutation R = R 23 R 13 R 12 . The inverted hierarchy limit is independent of permutations as two of the ω ij are zero. In Figure 11, we show the triangle plots for NH and IH in the symmetry protected region after applying the aforementioned cuts to remove fine-tuned points. The results are consistent with what was found in Ref. [178] for n = 2 RHNs.

Discussion of excesses likelihoods
In the previous subsections, we have made use of an exclusion-only 'capped' profile likelihood to study the constraining effect of the various observables on the parameter space (for a justification see Sec. 4.5). The total likelihood, however shows a pattern of excesses in some small regions of the parameter space. As discussed in Sec. 4.5, experimental results with a preference for specific heavy neutrino masses and mixings would in general not show up as localized excesses in the total profile likelihood. This is due to the fact that for each value of M I it would be in general possible to find a value of M J (with J = I) and associated couplings that would maximize the excess likelihood, irrespective of the values of M I . In order to extract the specific masses and couplings preferred by an excess likelihood, we adopt throught this subsection the following strategy. We only allow one of the three RHNs (which we take to be I = 1) to acquire the required masses and couplings, while disallowing the other two RHNs to enter the preferred region. This is emphasized in the plots by specifying M 1 and |U α1 | 2 instead of M I and U 2 αI . Mind that these results would be identical for M 2 and M 3 .
The invisible width of the Z boson is modified by the presence of the right-handed neutrinos through their mixing, as described in Section 3.2. For very high τ couplings, U 2 τ I > 10 −3 , the prediction from the RHN model is actually a better fit to the experimental measurement than the SM, and thus there is a slight (< 2σ) excess. A similar effect occurs for the CKM and R τ eµ constraints, where the modified contribution on the neutrino mixing in the decay products of K-mesons and τ , enhances the prediction with respect to that of the SM. Figure 12 shows the excesses on the total profile likelihood in the M 1 vs U 2 τ I plane, zoomed in at high couplings (as discussed above, we excluded M 2 and M 3 from entering the excess regions). Since there are no constraints from direct searches at masses above M 1 > 80 GeV or in the range 0.3 < M 1 < 0.5 GeV, there is a combined excess shown of about 2σ.
In order to study the impact of the different partial likelihoods on the total likelihood excess, we show in Figures 13 and 14 the partial one-dimensional likelihoods for Γ inv (blue), CKM (green) and R τ eµ (pink) with respect to the total likelihood (red) for M 1 < 1 GeV and M 1 > 60 GeV, respectively. All likelihoods are normalised so that they show up as a bump over the combination of all other likelihoods L 0 (grey). These plots show that the combination of excesses from all three sources amounts to a deviation of around (high mass) or above (low mass) 2σ with respect to the background. As observed in the figures, the effect of R τ eµ is rather negligible compared to the other two relevant likelihoods. Even larger couplings are severely penalised by the steep drop in the Γ inv likelihood. Figures 12-14 in |U τ 1 | 2 , for both low and high masses, are the most significant excesses arising in our three RHN scenario, but not the only ones. At masses around the K-meson resonance, there is an even dimmer excess in |U e1 | 2 , arising from the constraint on fully leptonic decays of K-mesons, R K eµ . As seen in Figure 15, for both normal and inverted hierarchy, there is a ∼ 1σ excess at M 1 ∼ 0.45 GeV. As before, we show in Figure 16 the one-dimensional likelihoods for R K eµ (purple) with respect to the total likelihood (red), over the background of the combination of the rest of constraints (grey). Although the R K eµ likelihood keeps increasing for larger values of |U e1 | 2 , the total likelihood drops at the limit shown in the figures  due to the constraints from the CHARM experiment (orange).

The excesses shown in
Although the identified excesses provide interesting hints towards specific regions of the RHN parameter space, they should not be over-interpreted, since their significance remains rather small and probably consistent with statistical fluctuations. The presence of such excesses was already observed before, identified in EW-POs [23] (cf. also [183]) and CKM constraints, and particularly in τ → s transitions [28].

Conclusions & Outlook
We presented here the first frequentist global analysis of the extension of the Standard Model by three heavy right-handed Majorana neutrinos for a large range of their masses, from 60 MeV to 500 GeV, and for normal and inverted hierarchy of the active neutrino masses. As detailed in Section 1, our analysis improves on previous studies in numerous ways. Most notable is the inclusion of a larger number of experimental constraints than in previous studies, such as EWPOs, all LFV decay channels, active neutrino mixing and masses, as well as many direct searches. Furthermore, we have performed a proper statistical combination of all constraints using a composite likelihood approach, and studied the overall constraints on the parameter space using robust profile likelihood methods. To this end, we have used the advanced BSM inference tool GAMBIT [37], which  we appropriately extended with the relevant model specifications and experimental constraints.
The results shown in Section 5 cover the full studied mass range for all couplings down to U 2 I ∼ 10 −10 . The profile likelihood contours are consistent with the results found in previous studies. The upper limits on the heavy neutrino mixing with electron and muon flavour mostly follow the confidence levels provided by direct search experiments, while the mixing with the third generation is constrained by a combination of different direct and indirect searches. We explicitly studied the limit of vanishing lightest neutrino masses, where we have shown that the flavour mixing pattern becomes significantly constrained for small values of m ν0 . For m ν0 < 0.01 meV these constraints become independent of the precise value of m ν0 in both mass hierarchies, which suggests that one heavy neutrino has effectively decouples. In this regime we demonstrated that one can recover the results that have previously been found in the model with only two RHN in earlier works.
Furthermore, we identified a few excesses in the profile likelihood, which are due to the invisible decay width of the Z-boson, the CKM unitarity constraint and R K eµ . Our best fit has a significance (w.r.t. SM) slightly above 2σ. Although these excesses are not significant enough to favour the n = 3 right-handed neutrino model in favour of the SM at the moment, an improvement on the measurements of the relevant observables will increase/decrease their significance in the future. Future e + e − colliders, such as the ILC, FCC-ee or CEPC, might measure EW observables, including the Z decay width, with higher precision [184] than the current value from LEP [185]. The NA62 experiment, which targets kaon decays, might be able to improve the measurments of the CKM matrix elements V us and V ud , as well as the lepton universality ratio R K eµ through more precise measurement of the fully leptonic decays of kaons [186].
Since the strongest constraints on the absolute value of the couplings come from direct searches, it is expected that the results obtained in this analysis will change significantly with the next generation of direct search experiments. An overview of projected sensitivities can e.g. be found in Refs. [16][17][18]. Many of these searches can be performed at existing facilities, including the LHC, NA62, T2K or the DUNE near detector. The sensitivity of the LHC will soon be upgraded with the recently approved FASER experiment [187] and other proposed dedicated detectors [188,189,[189][190][191][191][192][193][194]. In the more distant future the SHiP experiment [195,196] can search for heavy neutrinos in the GeV mass range [197], while future folliders such as FCC [198] or CEPC [199] can explore larger masses. These experimental perspectives make the study of right handed neutrinos an exciting topic for the years to come. Additional motivation for such searches comes from cosmology because the baryon asymmetry of the universe can be explained by low scale leptogenesis for all experimentally allowed values of the mixing angles in the model considered here if the heavy neutrino masses lie below the electroweak scale [54]. If any heavy neutral leptons are found in experiments then our results for their properties, such as the flavour mixing pattern as a function of light neutrino parameters, provide a powerful test to assess whether these particles are responsible for the generation of light neutrino masses and/or the baryon asymmetry of the universe [200], and to distinguish the model with three heavy neutrinos considered here from the model with two heavy neutrinos or other extensions of the SM.

Appendix A: GAMBIT Implementation
GAMBIT 13 (the Global and Modular BSM Inference Tool) [37] is a global fitting software framework that allows for extensive calculations of observables and likelihoods in particle and astroparticle physics. It provides, out-of-the-box, a suite of statistical methods and parameter scanning algorithms, together with a hierarchical model database, a strong interface to external tools and a host of other utilities that make it one of the most powerful global fitting tools on the market.
In a nutshell, the fundamental building blocks of GAMBIT are its module functions, which calculate all physical and mathematical quantities revelant to an analysis. Each module function provides a capability which, together with the return type of the function, unequivocally specifies the quantity calculated.
These module functions are sorted in the different physics modules, according to their purpose, e.g. functions calculating dark matter relic density lie in Dark-Bit [201]. Since most observables and likelihoods computed for this analysis do not belong naturally to any of the existing GAMBIT modules, we introduce a new GAMBIT physics module, NeutrinoBit, which contains all calculations relating to (active and sterile) neutrino physics. All of the module functions and capabilites described below are implemented in the new module NeutrinoBit, unless otherwise stated.

A.1: Neutrino models
In GAMBIT, models are defined by a set of parameters and relations to other models [37]. All the SM and active neutrino parameters are defined in a model called StandardModel_mNudiff, a daughter model of the GAMBIT model StandardModel_SLHA2, which includes SM parameters written in the SLHA2 convention [202]. The former contains the parameters mNu_light, dmNu21 and dmNu3l, which give the lightest neutrino mass and mass splittings, respectively. Other parameters in this model that are relevant for this study and are scanned over, as described in Section 4, are alpha1, alpha2, delta13, theta12, theta23 and theta13.
The right-handed neutrino sector is defined in another model, RightHandedNeutrinos, which contains the RHN masses M I and the real and imaginary parts of the ω ij parameters in the C-I parametrisation (c.f. Section 2.4). These are M_1, M_2, M_3, ReOm12, ReOm13, ReOm23, ImOm12, ImOm13, ImOm23. To better explore the symmetry preserved region (Section 2.5), a differential model, inheriting from RightHandedNeutri- 13 gambit.hepforge.org. nos, has been defined, RightHandedNeutrinos_diff. This model swaps the parameters M_2 for delta_M21, and defines a translation function from the parameters of the daughter model to the parent model as Lastly, the parameter Rorder encodes the ordering of the matrices R ij in Eq. (30), which allows us to fully cover all the parameter space with the C-I parametrisation.
There is a number of useful quantities and observables that can be constructed from the neutrino parameters, and these are all implemented in NeutrinoBit.cpp. In the active neutrino sector we calculate the neutrino mass matrix, m_nu, and their mixing matrix, UPMNS, as well as useful quantities such as the type of hierarchy, ordering, the squared mass splittings, md21, md31 and md32, and the minimal neutrino mass min_mass. It is worth noting that it is possible to fix the hierarchy of a scan by providing the option ordering to the capability m_nu in the configuration file. For example, in order to scan the normal hierarchy we would define in the YAML file The right-handed neutrino sector also contains a couple of useful capabilities, SeesawI_Vnu, which is the active neutrino mixing matrix in type-I seesaw, effectively UPMNS corrected by the presence of the right-handed neutrinos, and SeesawI_Theta, the active-sterile mixing matrix in type-I seesaw, currently implemented using the C-I parametrisation.
Another useful capability defined here is Unitarity, which is fulfilled by two module functions according to whether the model scanned is the SM or a RHN model, and checks whether the full mixing matrix is unitary. All these capabilities relating to neutrino masses and mixings and the module functions that fulfill them, along with their dependencies and options can be seen in Table  8.
Lastly, NeutrinoBit.cpp also contains likelihoods for the active neutrinos, which are implemented following the results from the NuFit collaboration (c.f. Section 3.1). The capabilities associated with these are md21_lnL, md3l_lnL for the mass splittings, deltaCP_lnL, theta12_lnL, theta23_lnL and theta13_lnL for the phases and mixing angles, and sum_mnu_lnL for the cosmological limit on the sum of neutrino masses. All these capabilities, with their module functions and dependencies are listed in Table 9. Calculates the active-sterile mixing matrix in seesaw type-I using the C-I parametrisation.

Unitarity_SeesawI(bool):
Checks for unitarity in the full neutrino mixing matrix in seesaw type-I.  A.2: Right-handed neutrino likelihood functions Every observable and likelihood described in Section 3 has an assigned capability within GAMBIT. Most of these have been implemented in the new GAMBIT module NeutrinoBit, since they concern mostly neutrino physics. Their module functions are coded in the file RightHandedNeutrinos.cpp, to keep them separated from the likelihoods and observables concerning only active neutrinos in NeutrinoBit.cpp. The exception to this is the LFV observables and semileptonic lepton universality tests, which can be found in FlavBit [203], implemented in FlavBit.cpp and the electroweak precision observables, which are coded up in PrecisionBit.cpp in PrecisionBit [204]. The implementation details for each specific observable are as follows:

Electroweak precision observables
Mainly, the EWPO capabilities lie in the physics module PrecisionBit and the associated module functions are implemented in PrecisionBit.cpp. These capabilites, prec_sinW2_eff and mW, can be seen in Table 10 along with their module functions and dependencies. The log-likelihoods are provided by the capabilities lnL_sinW2_eff and lnL_W_mass can also be seen in the same Table. Additionally, the module DecayBit contains the capabilities for the invisible width of Z, which are Z_gamma_nu and lnL_Z_inv, and leptonic decays of the W boson, W_to_l_decays and lnL_W_decays, all of which can be seen as well in Table 10.

Lepton flavour violation
The capabilities related to lepton flavour violation can be found in FlavBit and are muegamma, tauegamma, taumugamma, mueee, taueee, taumumumu, taumuee, taueemu, tauemumu, taumumue, mueTi and muePb. Table 11 shows these capabilities, the module functions that provide them, implemented in FlavBit.cpp, and their dependencies. The likelihoods, shown in Table 12, are collated into three capabilites, according to the type of process, l2lgamma_lnL for l → lγ, l2lll_lnL for l − → l − l − l + and mu2e_lnL for µ−e conversion in nucleii.

Lepton universality
The observables and likelihoods associated with lepton universality constraints are spread between the NeutrinoBit and FlavBit modules. Those involving fully leptonic decays are implemented in RightHandedNeutrinos.cpp and those for semileptonic decays of B mesons are in FlavBit.cpp. The capabilities for leptonic decays are R_pi, R_K, R_tau and R_W, and for semi-leptonic RK, RKstar_0045_11 and RKstar_11_60. They can be seen in Table 13 together with the module functions that provide them and their dependencies. The capability LUV_M collates all semileptonic universality observables into the FlavBit-defined class FlavBit::predictions_measurements_covariances 14 . The capabilites that compute the likelihoods for lepton universality tests are lnL_R_pi, lnL_R_K, lnL_R_tau and lnL_R_W for leptonic decays, and LUV_LL for semi-leptonic, and they, the module functions and dependencies, can be seen in Table 14.

CKM unitarity
The NeutrinoBit capability calc_Vus, implemented in RightHandedNeutrinos.cpp, calculates the value of V us that maximizes the likelihood for a given 14 For more details about FlavBit types, see [203]. mixing matrix Θ. The capabilities lnLckm_Vus and lnLckm_Vusmin compute the log-likelihood using V us as a scan parameter and as calculated by the profiling of calc_Vus, respectively. The capabilities, module functions and dependencies defined in GAMBIT for the calculation of the observable and the likelihood connected to CKM unitarity are listed in Tab. 15.

Neutrinoless double beta decay
In NeutrinoBit there are two computations of the likelihood for neutrinoless double beta decay, one based on the half-life and one based on the invariant mass of the two electrons m ββ . The capabilities Thalf_0nubb_Xe and Thalf_0nubb_Ge calculate the half-life of the 0νββ process as studied with Xe and Ge detectors. Equivalently, the capabilities mbb_0nubb_Xe and mbb_0nubb_Ge compute m ββ for the process in Xe and Ge detectors. The log-likelihoods are computed according to the experiments: lnL_0nubb_KamLAND_Zen and lnL_mbb_0nubb_KamLAND_Zen calculate the log-likelihood for the KamLAND-Zen experiment based on the half-life and m ββ , respectively; and lnL_0nubb_GERDA and lnL_mbb_0nubb_GERDA for the GERDA experiment. Lastly, the total log-likelihood is given by the capabilities lnL_0nubb, constructed from the half-life, and lnL_mbb_0nubb, from m ββ . Tab. 16 shows the defined capabilities, associated module functions and dependencies related to neutrinoless double beta decay that are responsible for the calculation of observables and likelihoods.

Big Bang nucleosynthesis
There are a number of processes that contribute to the decay width of the right-handed neutrinos, relevant for Big Bang nucleosynthesis, and each of them is computed by a capability. These are Gamma_RHN2pi0nu, Gamma_RHN2nuddbar and Gamma_RHN2ludbar. The total decay width of each of the right-handed neutrinos is given by Gamma_BBN and the log-likehood for BBN by lnL_bbn.
Tab. 17 shows the capabilities, as defined in GAMBIT, that pertain to Big Bang nucleosynthesis, and the module functions that satisfy them, along with dependencies that other module functions fulfill. The decay process considered in each function is mentioned below its name.   and Ut3_phase. These can be seen in Table 18 where

RK RHN_RK(double):
Calculates the test of lepton universality R K .

Other capabilities
The theoretical constraint for perturbativity of the Yukawa couplings has been implemented in NeutrinoBit as well. The capability perturbativity_lnL calculates a step function likelihood for this constraint. Tab. 20 shows the module function that provides this capability and its dependencies. Lastly, the artificial coupling slide likelihood that was introduced to drive the scan towards high couplings, as described in Section 4, was also implemented in Neutrino-Bit with capability coupling_slide. The module function and dependencies of this capability can also be seen in Table 20.   Calculates the log-likelihood.   The decay widths of LFV processes, as described in Section 3.3.1, are given by [205,206]

Gamma_BBN
The couplings e, g 1 , g 2 correspond to the electromagnetic, hypercharge and weak couplings of the SM.
The form factors K X 1 , K X 2 , A S XY and A V XY are taken in the flavour basis where the charged lepton mass matrix is diagonal.
The dipole form factors K X 1 and K X 2 are given as [206] K L 1 = 0 (B.6) (B.8) The four lepton form factors A V XY and A S XY corresponding to the process l − α → l − β l − γ l + δ , with a vector or scalar mediator respectively, are [206] where sum over a and c is assumed, x a = m 2 νa /m 2 W , s w = sin θ w and g ± = g 1 sin θ w ± g 2 cos θ w .
The µ − e conversion ratio described in section 3.3.1, for a general nucleus N Z+N Z , can be written as where B K XY ↔ C K XY for up-type quarks and the numerical factors G K are given in [207]. The nuclear form factor F p , the effective atomic number Z eff and the capture rate Γ capt of the nucleus [208] can be seen, for the cases we are studying, Ti 48 22 and Pb 208 82 , in Table 21.  The form factors K X 1 , K X 2 are defined above and B K XY and C K XY are given by Lastly, the loop functions used in (B.7)-(B.24) are defined as [206]

Appendix C: Distinguishing symmetry protected from tuned parameter choices
One goal of the present work is to fully understand the experimentally allowed range of parameters for heavy neutrinos, with minimal theoretical bias. To achieve this we employ the Casas-Ibarra parametrisation (29) and adapt agnostic priors for the parameters in table 5. On the other hand it is also instrictive to understand what fraction of the parameter space can only be realised at the cost of fine tuning in the parameters. This requires to distinguish fine-tuned parameter choices from symmetry protected ones. The Casas-Ibarra parametrisation (29) is inherently motivated from "bottom up", and it is not easy to see directly from the values of its fundamental parameters whether they exhibit a symmetry protection. A full analytic exploration of all the possible solutions and their classification between symmetric and fine-tuned would be a useful exercise, but lies outside the scope of this work. In our numerical scan we take a more pragmatic approach. We first generate a huge amount of parameter choices by randomising the parameter values and the order of the matrices R ij to ensure a maximal coverage of the parameter space. We then use the cuts (75) to distinguish the symmetry protected points a posteriori. This cut practically enforces the structure (38) on the masses and couplings.
Using the cut (75) requires some care for two reasons. First, the form (38) does not capture all symmetry protected points, cf. footnote 6. We may therefore misidentify some symmetry protected points as tuned. We find, however, that the number of such points is small. Second, when using the Casas-Ibarra parametrisation, it is possible to generate points that mimic the form (38) and hence pass the cut (75), but in fact exhibit a significant amount of tuning.
To illustrate the second point we work at tree level and approximateM diag M M , which yields For the inverted hierarchy, it is straightforward to show that one qualitatively gets a pseudo-Dirac pair of heavy neutrinos for M 1 = M 2 =M , (ω 12 , ω 13 , ω 23 ) = (ω, 0, 0) (C. 61) with |Imω| 1. In addition we have to set m ν0 = 0 to find the symmetry protected region, as a non-zero lightest neutrino mass is not consistent with α → 0. In this case the upper left block of the matrix R √ M M in (29) is large and the third row that multiplies m 3 = m ν0 is small, thereby mimicking the structure in (38). One can therefore interpret the decoupling F α3 = 0 and the vanishing mass of the lightest neutrino physically as results of the symmetry.
If we choose normal ordering, then choosing (C.61) still yields a structure that passes the cut (75), but it is in fact a tuned solution that just mimics this structure. In that case m ν0 = m 1 multiplies the large components of R √ M M −1 in (29). Hence, the approximate symmetry makes the wrong light neutrino mass small (m 3 instead of m ν0 = m 1 ). Of course one can set m 1 = 0 by hand in (29), but this choice cannot be justified by the symmetry. Though the limit (C.61) leads to a pseudo-Dirac structure amongst the ν Ri as predicted by the B −L symmetry, the vanishing mass of the lightest neutrino is not a result of that symmetry. The problem is that the Casas-Ibarra parametrisation allows on to set m ν0 = 0 by hand and gives no warning if this leads to accidental cancellations.
Realising the symmetry requires to choose the eigenvalues of M M and the non-zero ω ij consistently in a way that the large block in the matrix R √ M M −1 in (29) multiplies the two non-zero light neutrino masses. For normal ordering this is achieved with again with m ν0 = 0. In particular, one cannot choose ν R3 to be the particle that decouples. This is clearly no fundamental problem because the labels of the ν RI have no physical meaning, but it means that the labelling and the order of the matrices R ij have to be taken into consideration when applying a cut to identify symmetry protected points in the numerical data. The situation is yet more tricky if one considers small perturbations around the choices (C.62) or (C.61). In our numerical scan we randomise the order of the three matrices R ij in (30) to generate more points. If one exactly takes the choice (C.62) or (C.61) for normal or inverted neutrino mass ordering, respectively, then the approximate B −L conserving limit is reproduced irrespectively of the ordering of the R ij . However, the effect that small perturbations around this limit have strongly depends on this ordering. The effect of perturbing R is the smallest if the matrices R ij are ordered in a way that the one with large entries (controlled by ω) directly multiplies M −1 M in (C.60). For normal ordering this is the case with R = R 23 R 13 R 12 , and for inverted ordering with R = R 12 R 13 R 23 . This procedure was crucial to reproduce the constraints on the heavy neutrino flavour mixing pattern in the n = 2 model found in ref. [178], cf. fig. 11.
It is worth noting that (C.61) and (C.62) are not the only combinations of parameters that yield the symmetry protected scenario, but rather the simplest. The non-trivial structure of the complex rotation matrix R yields many solutions to the required block layout described before. In fact, we will take advantage of this fact further below to recover the n = 2 case from the n = 3 Lagrangian for normal ordering by taking M 1 = M 2 =M , (ω 12 , ω 13 , ω 23 ) = (0, π/2, ω), (C. 63) since in this work we will focus mainly on the case where M 1 and M 2 are almost degenerate.

Appendix D: Partial likelihoods
The final result of a global fit shows the combined effect of all likelihoods on the parameter space of the model. It is, however, often useful to understand the effect on the individual partial likelihoods. Therefore, we show here a comprehensive set of scatter plots that show the contribution of each relevant partial likelihood in the M 1 vs |U α1 | 2 . As we have seen before, away from the massless neutrino limit there is little difference between NH and IH, and thus we will only show the partial likelihoods for inverted ordering.  As seen in Figure 21, direct searches from CMS compete in a small mass range with DELPHI prompt searches, the statistical combination of the two setting stronger limits than each of them individually. Recent  and future results from CMS and ATLAS not included in this study are expected to dominate in this range. Figure 22 shows that the larger mass range is unconstrained by direct searches, hence electroweak precision observables, in particular sin θ W , are responsible for the upper limits in this range.
Similar to the case above, the coupling |U µ1 | 2 is constrained from above by several direct and precision  the results from CHARM; the long-lived particle search from DELPHI remains unchallenged for M I ∼ (2, 4) GeV whereas, as before, the DELPHI prompt search competes in the range M I ∼ (4, 80) GeV, with searches at CMS.
Larger masses are not constrained by direct searches, but rather by a combination of precision limits. Contrary to |U e1 | 2 , where only sin θ W dominated at large masses,  upper values of |U µ1 | 2 are also mildly constrained by lepton flavour violating decays, particularly µ → eγ. Hence, for M I 80 GeV, the combination of EWPO, sin θ W , and LFV decays, are the most constraining, as seen in Figures 30 and 31.
The couplings of heavy neutrinos to the τ flavour, |U τ 1 | 2 , are not as strongly constrained from above by direct searches. In Figures 32-33, one can see that for low masses, M I 0.3 GeV, only the direct searches from CHARM in the τ channel set an upper limit on the couplings. In the mass range M I ∼ (0.5, 80) GeV, longlived and prompt searches by DELPHI dominate. At low masses, the |U τ 1 | 2 coupling is constrained from below, as seen in Figure 35. This lower bound is a consequence of BBN, which sets a lower limit on the sum of couplings  |U 1 | 2 , and PS191, which forces the e and µ couplings to be small at low masses.
In the mass range M 1 ∼ (0.3, 0.5) GeV, as well as for large masses M 1 80 GeV, direct searches do not constrain |U τ 1 | 2 . Hence in these ranges, the strongest constraints come from the invisible decay of the Z boson, as seen in Figure 36. This figure uses the "capped" likelihood defined previously, so the excesses in Γ inv discussed in Section 5 will not be visible.