Global constraints on heavy neutrino mixing

We derive general constraints on the mixing of heavy Seesaw neutrinos with the SM fields from a global fit to present flavour and electroweak precision data. We explore and compare both a completely general scenario, where the heavy neutrinos are integrated out without any further assumption, and the more constrained case were only 3 additional heavy states are considered. The latter assumption implies non-trivial correlations in order to reproduce the correct neutrino masses and mixings as observed by oscillation data and thus some qualitative differences can be found with the more general scenario. The relevant processes analyzed in the global fit include searches for Lepton Flavour Violating (LFV) decays, probes of the universality of weak interactions, CKM unitarity bounds and electroweak precision data. In particular, a comparative and detailed study of the present and future sensitivity of the different LFV experiments is performed. We find a mild $1-2\sigma$ preference for non-zero heavy neutrino mixing of order 0.03-0.04 in the electron and tau sectors. At the $2\sigma$ level we derive bounds on all mixings ranging from 0.1 to 0.01 with the notable exception of the $e-\mu$ sector with a more stringent bound of 0.005 from the $\mu \to e \gamma$ process.


I. INTRODUCTION
The present evidence for neutrino masses and mixings demands an extension of the Standard Model (SM) to account for them and represents one of our best windows to new physics. The simplest and one of the most appealing possibilities, given its symmetry with the quark sector, is the addition of fermion singlets, right-handed neutrinos, to the SM particle content. However, even this simplest extension points towards the existence of a new type of term in the SM Lagrangian: a Majorana mass term for the right-handed neutrinos, allowed by the SM gauge symmetry. Such a term would imply direct breaking of the otherwise accidental Lepton number symmetry and the introduction of a mass scale not directly related to the Higgs mechanism of electroweak symmetry breaking. Depending on the actual scale of this mass term, interesting phenomenological consequences follow, such as the possible explanation of the baryon asymmetry of the Universe via the Leptogenesis mechanism [1] or of the mysterious dark matter component via these sterile neutrinos [2][3][4][5].
A popular assumption for this mass scale is that it is above that of electroweak symmetry breaking. This choice indeed leads to the well-known Seesaw mechanism [6][7][8][9] that nicely accommodates the strikingly tiny neutrino masses, as compared with the rest of the SM fermion content. In particular, for neutrino Yukawa couplings ranging in value from the electron to the top quark, a Majorana mass scale between the electroweak (EW) and grand unification scales can correctly reproduce our present constraints. Unfortunately, this huge hierarchy of scales also suppresses any other observable consequence of the model, beyond the leading order Weinberg d = 5 operator [10] that explains neutrino masses, rendering its experimental verification extremely challenging.
An interesting alternative is that of explaining the smallness of neutrino masses, not through a large hierarchy of scales, but rather via an approximate symmetry [11][12][13][14][15][16]. In particular, there are choices for the Majorana mass and Yukawa matrices such as the inverse [11,12] or linear [17] Seesaw mechanisms that, for a given assignment of the charges among the extra states, approximately conserve B − L. Therefore, the Majorana masses obtained via the Weinberg operator by the light neutrinos are necessarily suppressed by the small B − L-violating parameters. Interestingly, higher order operators that would, a priori, be more strongly suppressed than neutrino masses, are not necessarily protected by this symmetry and can thus lead to sizable signals. In particular, apart from the Majorana nature of neutrino masses, the most characteristic signals of the Seesaw mechanism with right-handed neutrinos are deviations from unitarity of the lepton PMNS mixing matrix generated by the only d = 6 operator present at tree level. This in turn would lead to signals in lepton flavour violating processes (LFV), non-universality of weak interactions and/or affect electroweak precision observables .
In this work we will combine results from all these probes to derive updated constraints on the presently allowed mixing among the extra massive neutrinos and the SM flavour eigenstates. We will present our results both, for a completely model-independent parametrization without any further assumption about the extra massive states and for the more restricted assumption of only three massive neutrinos, in analogy to the three generations for all other fermions. In the latter case, since the three extra neutrinos must also reproduce the correct pattern of masses and mixings as observed in neutrino oscillations, correlations among the potentially observable effects are predicted and constraints are qualitatively different from the general case. This paper is organized as follows. In Section II we introduce the parametrizations adopted for our studies for the general and three-heavy-neutrino cases. In Section III we describe the set of observables used to probe for the heavy extra neutrinos. In Section IV we present and describe our results and finally we conclude in Section V.

II. PARAMETRIZATIONS
Starting from the usual type-I Seesaw Lagrangian: where φ denotes the SM Higgs field, M N the Majorana mass allowed for the right-handed neutrinos N i R and Y N the Yukawa couplings between the neutrinos and the Higgs field. The vev of the Higgs v EW will, in addition, induce Dirac masses m D = v EW Y N / √ 2. In the usual Seesaw limit, for M N m D , the three light and mostly-active neutrinos observed in the neutrino oscillation phenomenon will be clearly separated from the heavy and mostly-sterile new states. Upon integrating out these heavy states, their low energy phenomenology will be encoded in a series of effective operators. The first such operator is the well-known d = 5 Weinberg operator [10] that, upon electroweak symmetry breaking, induces the Majorana masses for the light neutrinos: where U PMNS = U 23 (θ 23 )U 13 (θ 13 , δ)U 12 (θ 12 )diag(e −iα 1 /2 , e −iα 2 /2 , 1) is the Unitary mixing matrix that diagonalizes the symmetric mass matrixm generated from the Weinberg operator. At tree level, the only d = 6 operator obtained upon integrating out the heavy neutrinos induces non-canonical neutrino kinetic terms for the three SM active neutrinos when the Higgs develops its vev [44]. After diagonalizing and normalizing the kinetic terms, the mixing matrix appearing in charged current interactions will thus contain, not only the two Unitary rotations to diagonalize the d = 5 and d = 6 operators respectively, but also the necessary rescaling to bring the neutrino kinetic term to its canonical form. Thus, in all generality, the matrix describing the mixing between the light neutrino mass eigenstates and the SM charged leptons via W interactions will not be Unitary and to stress this feature we will dub it N . Since any general matrix can be parametrized as the product of an Hermitian and a Unitary matrix, these deviations from unitarity have been often parametrized as [45]: where the small Hermitian matrix η (also called in other works) encodes the deviations from unitarity in neutrino mixing. This parametrization is very convenient from a phenomenological point of view. Indeed, since the particular neutrino mass eigenstate is never identified in physical observables, its index is always summed upon, while the flavour index labeling the charged leptons participating in the process is normally fixed. Thus, most observables depend on the combination: and can thus be expressed only though the parameters contained in the Hermitian matrix η. Moreover, the physical interpretation of η is also very transparent in terms of the mixing between the extra heavy neutrinos and the SM flavours. Indeed, if the full mass matrix is diagonalized as: where m and M are diagonal matrices containing respectively the masses of the 3 light ν i and heavy N i mass eigenstates. The diagonalizing matrix U can be written as [46]: where and Θ ∼ m † D M −1 N is the general matrix that describes the mixing between the heavy mass eigenstates and the active neutrino flavours. Thus, the non-unitary correction I − η can be identified with the first term of the cosine expansion 1 − ΘΘ † /2 such that: Furthermore, η is also (1/2 of) the coefficient of the d = 6 operator obtained upon integrating out the heavy neutrino fields: In all generality the d = 6 operator η is completely independent from the d = 5m and thus from the measured neutrino masses and mixings in oscillation experiments [47,48]. However, bothm and η are ultimately built from m D and M N and thus, in particular cases, may not be fully independent. Apart from the completely general parametrization through η, here we will also investigate one such case. Namely, we will focus on the particular scenario in which: • The SM is only extended through 3 right-handed neutrinos.
• The three extra neutrino mass eigenstates are heavier than the EW scale.
• Large, potentially observable, η is allowed despite the smallness of neutrino masses.
• The small neutrino masses are radiatively stable.
The only way to simultaneously satisfy these requirements is through an underlying L symmetry [49,50] (see also Ref. [40,51]) which leads to: with all i and µ j small lepton number violating parameters (see also Ref. [52] for a particular scenario where these small parameters arise naturally). By setting all i = 0 and µ j = 0, lepton number symmetry is indeed recovered with the following L assignments L e = L µ = L τ = L 1 = −L 2 = 1 and L 3 = 0. Alsom = 0 (3 massless neutrinos in the L-conserving limit), M 1 = M 2 = Λ (a heavy Dirac pair) and M 3 = Λ (a heavy decoupled Majorana singlet), but: So that large η is possible even in the limit of massless neutrinos when L is conserved. Upon switching on the L-violating parameters in Eq. (10), neutrino masses and mixingsm that can reproduce the observed neutrino oscillations are generated. However, these are not completely independent from η and the following relationship between the θ α in Eq. (11) andm follows [53]: Thus, this extra constraint will lead to correlations among the heavy-active mixing parameters θ α and therefore also η αβ through Eq. (11), not present in the completely general scenario with more than 3 heavy neutrinos. From now on we will refer to the unrestricted scenario as G-SS (general Seesaw ) and to the particular case with 3 extra heavy neutrinos as 3N-SS. The parameters characterizing the heavy neutrino mixing and the correlations between them in each case are summarized in Table I. In particular, the constraints on η for the G-SS come from the fact that η is positive definite (see Eq. (8)).Regarding θ τ in the 3N-SS case, its value is fixed by θ e and θ µ through Eq. (12) once the SM neutrino masses and mixings encoded in the d = 5 operatorm are specified. In our analysis we will thus scan the allowed parameter space of the 3N-SS by leaving θ e and θ µ free in the fit, together with the remaining unknown values characterizingm: the Dirac phase δ, the Majorana phases α 1 and α 2 , the absolute neutrino mass and the mass hierarchy (normal or inverted). Regarding the absolute neutrino mass scale we will add the constraint from Planck on the sum of the light neutrino masses m i < 0.23 at a 95% CL [54]. The rest of the oscillation parameters free free fixed by Eq. (12) fixed by θ e , θ µ fixed by θ e , θ τ fixed by θ µ , θ τ  The rest of the oscillation parameters are fixed to their best fits from Ref. [55].
are fixed to their best fits from Ref. [55] since they are well-constrained by present neutrino oscillation data. When presenting the results of the global fit in Section IV we will derive constraints on the mixing of the heavy neutrinos with the SM active flavours θ α in Eq. (11) for the 3N-SS. Regarding the G-SS, we do not specify the number of heavy neutrinos with which the SM is extended since all the observable effects are simply encoded in the matrix η. Thus, each heavy neutrino can have a different mixing Θ αi and, to ease the comparison with the results from the 3N-SS, we will use the combination √ 2η αα which represents the total mixing from all the additional heavy neutrinos with the flavour α and an upper bound on the individual mixings Θ αi :

III. OBSERVABLES
Global constraints on the mixing between the heavy and active neutrinos will be derived through a fit to the following 28 observables: • Ratios of weak decays constraining EW universality: R π µe , R π τ µ , R W µe , R W τ µ , R K µe , R K τ µ , R l µe and R l τ µ • 9 weak decays constraining the CKM unitarity The dependence of each observable on the non-unitarity mixing matrix N αi and the parameters η αβ will be presented and discussed in this section. In Ref. [53] it was recently shown that loop level corrections involving the new degrees of freedom can be safely neglected. However, many SM-mediated loop corrections are relevant for these precision observables and will therefore be accounted for [56]. Notice that, in principle, these SM loop corrections also contain an indirect dependence on the non-unitarity parameters, notably through their dependence on G F as determined in muon decay. This subleading dependence of the observables will be neglected and only the corrections from non-unitarity affecting the tree level relations will be discussed in the following expressions. The numerical analysis, however, contains all relevant SM loop corrections when comparing with the corresponding observables. The loop-corrected SM expectation, together with the leading non-unitarity correction and the experimental measurements that will be the inputs of our global fit are all summarized in Tab.(II).
A. Constraints from µ decay: G F , M Z , M W and θ W As usual, all SM predictions will be made in terms of the very accurate measurements of α, M Z and G F as measured in µ decay, G µ [56]: However, a non-unitary mixing matrix N αi would modify the expected decay rate of µ → eνν. Indeed, since the final state neutrinos are not determined, their index must be summed upon obtaining: Thus, G F as determined through muon decay (G µ ) acquires a non-unitary correction that will propagate to most observables: In particular, the relation between G µ and M W allows to constrain η ee and η µµ through kinematic measurements of M W : .
Similarly, the weak mixing angle s 2 W will be modified and independent determinations of s 2 W will be used to further constrain η ee and η µµ :  Regarding different measurements of s 2 W it is important to note that in some low energy determinations, such as from the weak charge of the proton or Møller scattering, the dependence on this parameter appears through the following combination −1/2 + 2s 2 W . Since the value of s 2 W is close to 1/4, there is a partial cancellation in this observables that, in the SM, allows for a very accurate determination of s 2 W , since small changes in its value significantly affect the degree of the cancellation and hence the size of the observable. For the same reason, we find that these observables are also very sensitive to corrections of the order of SM loop corrections times the non-unitary parameters η. Indeed, including some of these corrections we find that the corresponding coefficients in front of the η parameters in Tab.(II) would vary up to a factor 2, indicating that our approximation of neglecting these terms is not good enough for these precision observables. Since the inclusion of these corrections is beyond the scope of this work, we choose not to include these particular determinations of s 2 W in the list of observables for our global fit.
B. Constraints from Z decays

Z decays into charged fermions
The Z decays into charged fermions are not directly modified in presence of heavy neutrinos or a non-unitary lepton mixing matrix at tree level. However, these measurements depend on G F and s W and, as such, an indirect dependence on the non-unitarity parameters appears through its determination via muon decay, as described above. In particular: where the vector and axial-vector form factors are given by: with N C the color factor, N C = 3 (1) for quarks (leptons) and where Q f and T f are the electric charge and third component of the weak isospin of the fermion f . Notice that an additional dependence on η ee and η µµ will be present in g V through s 2 W and Eq. (18). The usual combinations of decay rates will be used as observables for the global fit: where Γ had ≡ q =t Γ q .

Invisible Z width
In presence of a non-unitary lepton mixing matrix N αi , the Z coupling to neutrinos is directly affected and becomes non diagonal since (N † N ) ij = δ ij . Thus, apart from its indirect dependence through G F , the invisible width of the Z, from which the number of active neutrinos can be determined, is directly sensitive to the mixing of heavy neutrinos: Notice that, since η αβ is positive definite from Eq. (8), the number of active neutrinos as measured through the invisible Z width will be smaller than 3 in presence of mixing with heavy neutrinos, to be compared with the present determination of N ν = 2.990 ± 0.007 from LEP [57].

C. Constraints from weak interaction universality tests
The lepton flavour universality of weak interactions is strongly constrained through ratios of lepton and meson decays differing in the charged lepton generation involved, such as π → µν i vs π → eν i . Since the final state neutrino cannot be determined, these processes are proportional to i |N αi | 2 ≈ 1 − 2η αα , where α is the flavour of the charged lepton. Thus, a flavour dependence is induced in presence on non-unitary mixing and the weak interaction universality constraints become powerful probes of heavy neutrino mixing: where the ratio of the SM expectations for the decay widths Γ SM α will be given by a function of the charged lepton masses involved containing the corresponding phase space and chirality flip factors as well as the different loop corrections. Thus, at tree level and for the particular case of π decays: Constraints on the values of the ratios of weak coupling constants R αβ as defined in Eq. (23) have been derived through ratios of different decays [58] and are summarized in Tab. (II).

D. Unitarity of the CKM matrix
The presence of extra heavy neutrinos leads to unitarity violations of the lepton PMNS mixing matrix leaving the CKM quark mixing unaffected. However, the processes through which the elements of the CKM matrix V are determined are affected both directly (for processes involving leptons) and indirectly (through the determination of G F in muon decays). In particular, the unitarity relation among the elements of the first row of the CKM matrix is very strongly constrained and reads: For the present accuracy on V us , the value of V ub = (4.13 ± 0.49) × 10 −3 [56] can be safely neglected in Eq. (25). This relation, together with the measurements from the different processes used to constrain V ud and V us will thus also present indirect sensitivities to η αβ .
In particular we will rewrite through Eq. (25): and use the following experimental constraints to fit for V us and the η αβ parameters on which they depend. In our final constraints on η αβ the dependence on V us has been treated as a nuisance parameter and the χ 2 has been minimized with respect to it.

Superallowed β decay
Superallowed β decays provide the best determination of |V ud |. However, in presence of a non-unitary PMNS matrix it will receive a direct correction with (1 − 2η ee ) from the electron and neutrino coupling, as well as the indirect correction from G F in Eq. (16). All in all the value of V ud extracted from this process corresponds to: The most recent update on V β ud based on 20 different superallowed β transitions [59] is listed in Tab. (II) and will be an input for our fit.
2. |V us | |V us | can be determined through τ decays and semileptonic or leptonic K decays. The values of f + (0) and f K /f π involved in these observables have been taken from [60].
• K decays Kaon decays offer a direct way to determine |V us |. Apart from their sensitivity to this parameter, decays with µ (e) final states also have a direct dependence on η µµ (η ee ) which cancels against the indirect dependence through G µ leading to: The present determinations of V K→πeνe us and V

K→πµνµ us
are listed in Tab.(II) and have been obtained from [61,62] together with f + (0) from [60], the correlation matrix among observables from [61] has also been taken into account.
An alternative determination of |V us | stems from the ratio of the branching fractions B (K → µν) /B (π → µν). Notice that in this ratio any direct or indirect dependence on leptonic non-unitarity cancels allowing to constrain the ratio |V us | / |V ud | as in the SM. Since this measurement is latter combined with V β ud from Eq. (27) to obtain V K,π→µν us the same (1 + η µµ ) correction as for V β ud is finally present: • τ decays An alternative constraint on |V us | can be obtained from the τ → Kν τ decay rate.
In presence of non-unitary leptonic mixing, a direct correction by (1 − 2η τ τ ) will be present from the τ coupling as well as the indirect correction from G F leading to the following dependence: The value of V τ →Kντ us is given in Tab. (II) [63].
Another possibility is to constrain |V us | from the ratio B (τ → Kν τ ) /B (τ → πν τ ). In complete analogy to Eq. (30), the sensitivity to the non-unitarity parameters takes the form: All these observables with the values listed in Tab. (II) will be used to fit for η ee , η µµ and η τ τ . Regarding |V us |, its value will be free to vary in the fit and will be treated as a nuisance parameter, choosing the value of |V us | that minimizes the χ 2 for each value of η ee , η µµ and η τ τ .

E. LFV observables
Flavour transitions α → β in presence of non-unitary mixing such that (N † N ) αβ = −2η αβ = 0 are no longer protected by the GIM [64] mechanism. Thus, the stringent constraints that exist on lepton flavour violating (LFV) processes translate into strong probes of the PMNS unitarity, in particular on the off-diagonal elements η αβ . Notice that from Eq. (8) η is a positive-definite matrix and its off diagonal elements subject to the Schwarz inequality: as summarized in Table I. Thus, the direct constraints on the diagonal elements of η stemming from the processes discussed above also constrain indirectly the size of the off-diagonal entries. Moreover, for the 3N-SS, Eq. (11) implies that the Schwarz inequality is saturated to an equality. Therefore, in the G-SS a global fit to constrain the diagonal elements of η with the list of observables described above will be performed. Then, constraints on the off-diagonal entries will be derived indirectly through the Schwarz inequality and compared with the direct bounds from LFV processes. For the 3N-SS, the LFV observables will be added directly to the global fit since they also constrain the diagonal elements through the saturation of the inequality. Below we list and describe the set of LFV transitions that would take place through non-unitary leptonic mixing. The present experimental bounds and future sensitivities are summarized in Table. III. A comparison summarizing the present relative importance of these observables constraining the off-diagonal elements of η (solid lines) is presented in Fig. 1. Since the LFV observables typically depend on the value of the heavy masses, we  Table. III. The red-shadowed region represents the non-perturbative region with |Y N | 2 > 6π. In the bottom panel, given the preference for non-zero h → τ µ [65,66] we show the preferred value in blue and the the 1σ region in yellow.
have performed the comparison for the 3N-SS, since there is only a common scale that simplifies the comparison. As can be seen, radiative decays l α → l β γ presently dominate the existing bounds and will thus be added to the global fit in the 3N-SS. However, regarding future expectations (dotted lines), the constraints on |η eµ | will be dominated by µ → eee or µ − e transitions in nuclei rather than by µ → eγ. On the other hand, the present and future sensitivity to |η eτ | and |η µτ | is completely dominated by the radiative decays l α → l β γ. In particular, the constraints on |η αβ | from the LFV decays of the Z and Higgs bosons, Z → l α l β and h → l α l β , are at least one or three orders of magnitude weaker than the bounds from radiative decays respectively. Unfortunately this precludes the explanation of the present mild preference for non-zero h → µτ [65,66] through heavy neutrino mixing (see yellow band in the lower panel of Fig. 1). Indeed, the values of the Yukawas required to explain these events are, not only excluded by the other observables depicted in the third panel of Fig. 1, but also fall into the non-perturbative region, shaded red in the figure. µ → e (Al) − < 10 −17 [78] µ → e (Ti) < 4.3 · 10 −12 [79] < 10 −18 [80] TABLE III: Summary of the present constraints and expected future sensitivities for the different LFV observables considered.
As shown in Fig. 1, at present l α → l β γ is able to set bounds much stronger than through this process.

LFV h decays
In the case of the LFV Higgs decay the expression at O η 2 αβ for the branching ratio is much more involved than in the Z → l ∓ α l ± β case. In Fig. 1 we have used the complete computation presented in [82][83][84]. Nevertheless, we instead present here an approximate expression which can be useful in order to understand the dependence on the parameters in the 3N-SS. where and . This approximate result is reasonably accurate for scales above few TeV and works very well for Λ 10 TeV. However, since here we are neglecting O (M 2 W /Λ 2 ) contributions, it fails for Λ 1 TeV. In any case, the full calculation shows that the constraints on |η αβ | are still very far from the present radiative bounds, falling indeed in the non perturbative region. 3. l α → l β l β l β decay Another LFV observable that would be induced by heavy neutrino mixing is the l α → l β l β l β . Its branching ratio, for the 3N-SS, is given by [85] Notice that, while additional non-unitarity corrections from G µ and s 2 W (also through Γ α when α = µ) would be present, these are higher order in η and therefore subleading since the whole process is already proportional to |η αβ | 2 . Fig. 1 shows that the present µ → eee decay bound on |η eµ | is quite competitive with the one coming from µ → eγ. The constraint is presently dominated by µ → eγ, but it is expected to be overcome by µ → eee in the future. On the other hand, the present and future sensitivity to |η eτ | and |η µτ | is dominated by the radiative decays.

µ → e conversion
In the 3N-SS, the ratio between µ → e conversion rate over the capture rate Γ capt in light nuclei is given by [34] where A corresponds to the mass number, Z (Z eff ) stands for the (effective) atomic number, F p is a nuclear form factor and The bounds shown in Fig. 1 have been obtained from µ → e conversion transitions in 27 13 Al and 48 22 Ti. The input values for the nuclear parameters F p , Z eff and Γ capt have been extracted from [86,87] and are summarized in Table 1 of [34].
According to the forecasted performances the future sensitivity to |η eµ | will be dominated by this observable. Remarkably, future µ → e searches [80] could improve the present bound by three orders of magnitude making it a very promising channel to probe for new physics signal in LFV decays.

Radiative decays
In the G-SS, the branching ratio for the radiative decays l α → l β γ is given by: For M k M W the limit can be simplified to: This expression shows how the non-unitarity induced in the PMNS by the heavy neutrinos and the separation of the two scales prevents the GIM cancellation. Indeed, the cancellation is recovered in the limit x k 1. These radiative decays are the observables dominating the present constraints on η αβ as shown in Fig. 1 and will thus be the ones introduced in the fit through Eq. (42) for the 3N-SS. In the G-SS, these constraints will be compared with the bounds stemming from the Schwarz inequality Eq. (33) from the outcome of the global fit to the diagonal entries.

IV. RESULTS
With the list of observables described in the previous section and under a Gaussian approximation we construct a χ 2 function to scan the parameter spaces of the G-SS and the 3N-SS. For the G-SS the free parameters of the fit are directly η ee , η µµ and η τ τ without further constraints and all the observables listed in Section III except for the LFV transitions will be used to constrain them. The LFV radiative decays rather constrain the off-diagonal elements of the matrix η. Therefore, to obtain the global constraints on the off-diagonal elements, the LFV radiative decays will be combined and compared with the indirect bounds implied by the Schwarz inequality Eq. (33) from the lepton flavour conserving observables.
Regarding the 3N-SS, the free parameters for the fit are θ e and θ µ (modulus and phase) while θ τ is given by Eq. (12) once the light neutrino masses and mixings are specified through the d = 5 operatorm. Thus, we also take as free parameters of the fit the values of the unknown phases of the PMNS matrix Dirac (δ) and Majorana (α 1 and α 2 ) as well as the mass of the lightest neutrino mass eigenstate for both a normal and an inverted neutrino mass ordering. The rest of the oscillation parameters are fixed to their best fits from Ref. [55] since they are well-constrained by present neutrino oscillation data. Notice that, a priori, the number of free parameters we fit for in the 3N-SS case is larger than in the G-SS. However, this larger number of parameters is only included to take into account the constraints affecting θ τ (and therefore η τ τ ) via Eq. (12) that are absent in the G-SS. Indeed, as we will see from the results of the fit, these constraints imply extra correlations between the parameters of the 3N-SS and there is in fact less freedom in the relevant parameters η ee , η µµ and η τ τ to fit for the observables. Since for the 3N-SS the Schwarz inequality Eq. (33) is saturated |η αβ | = √ η αα η ββ , the LVF radiative decays also imply non-trivial constraints on the values of θ α and the diagonal elements η αα and will hence be included in the list of observables of the global fit. Notice that, under the approximation of Eq. (42), the LFV radiative decays do not depend on the Majorana mass scale. Therefore, since none of the observables for the G-SS or 3N-SS cases depend directly on the Majorana masses, the bounds on the mixing derived apply for any choice of the heavy neutrino masses above the electroweak scale.
In Fig. 2 we present our results from the global fit, performed by scanning the relevant parameter spaces through a Markov chain Monte Carlo algorithm. The results presented here correspond to the frequentist confidence intervals for 1σ, 90% and 2σ significance. We present the results directly in the heavy-active neutrino mixing θ α for the 3N-SS under the assumption of a normal neutrino ordering (middle panels) and inverted neutrino ordering (lower panels). To ease the comparison of the constraints, we present the results for the G-SS (upper panels) in the variable √ 2η αα , which can be identified with the total effective mixing of the different heavy mass eigenstates with the flavour α, see Eq. (13), and an upper bound on the individual mixing Θ αi of any additional heavy neutrino N i . As can be seen, while the bounds on the individual parameters are comparable in strength for the two scenarios, the constraints imposed by Eqs. (11) and (12) for the 3N-SS reflect in strong correlations for their allowed regions. In particular, µ → eγ imposes a very stringent constraint in the product θ e θ µ leading to the hyperbolic constraints in the middle-left and bottom-left panels of the figure and absent in the upper for the G-SS. On the other hand, in the middle and bottom-right panels of the figure non-trivial correlations between θ e and θ τ , absent in the upper-right panel for the G-SS, can be observed. This stems from the fact that θ τ is not free to take any value preferred by the observables, but constrained by θ e , θ µ and the neutrino masses and mixings through Eq. (12).
To summarize the results of the global fit we present in Fig. 3 the profiles of the ∆χ 2 obtained as a function of the individual θ α and minimized over all the other parameters. The 1 and 2σ regions are colored in red and blue respectively. As can be seen, the observables considered (notably the invisible width of the Z and M W ) overall show a mild (between 1 and 2σ) preference for some degree of non-unitarity θ ∼ 0.03 − 0.04. The constraints on the universality of the weak interactions, particularly from ratios of pion and lepton decays, prefer these unitarity deviations with non-vanishing mixing with the heavy neutrinos to take place in the electron and tau sectors. This preference is clear in the upper panels of Fig. (3), which show the constraints for the unbounded G-SS. But, even in the more constraint case of a 3N-SS (middle panels for normal hierarchy and lower panels for inverted), there is enough freedom to accommodate this general preference shown by the datasets considered. The more characteristic feature that distinguishes the 3N-SS from the G-SS in Fig. 3 is the constraint in θ µ which, for the 3N-SS shows a very non-Gaussian behaviour with a very stringent 1σ limit and a much milder 2σ bound comparable to the one found for the G-SS. The reason for the comparatively much stronger 1σ constraint stems from the very stringent constraint from µ → eγ, which for the 3N-SS imply either a very small θ e or θ µ . Together with the 1σ preference for non-vanishing θ e , this implies a very strong 1σ upper bound for θ µ . On the other hand, at the 2σ level θ e can be arbitrarily small and thus the bound on θ µ from µ → eγ is evaded. Regarding the G-SS, µ → eγ only constrains the element η eµ and not η ee or η µµ since, contrary to the 3N-SS, the Schwarz inequality Eq. (33) is not saturated. Regarding θ e and θ τ , the limits for the 3N-SS and the G-SS are much more similar between them. Indeed, despite the constraint from Eq. (12) on θ τ , the preferred value for this parameter in the 3N-SS does not show significant deviations with respect to the G-SS. However, non-trivial correlations among the Majorana phases α 1 and α 2 as well as among the phases of θ e and θ τ : α e and α τ when a normal neutrino mass ordering is assumed are required to satisfy Eq. (12). These phase correlations are shown in Fig. 4. Two interesting features can be observed: (i) The values of the PMNS Majorana phases such that α 1 − α 2 ∼ 2nπ are favoured (left plot); (ii) The data prefers values for the phases of θ τ and θ e which satisfy α τ − α e ∼ (2n + 1) π (right plot). In the IH case, we have not found any significant correlation among the phases.
Regarding the off-diagonal elements |η αβ |, we present in Fig. 5 the limits obtained from the combination of all observables as a function of 2|η αβ | and marginalized over all the FIG. 5: Bounds on the off-diagonal entries of η αβ (|θ α θ β | for the 3N-SS). The upper panels are for the G-SS, and the middle and lower panels for the 3N-SS for a normal and inverted hierarchy respectively. For the G-SS the strongest limit between the direct bound from radiative LFV decays and the indirect limit from the diagonal entries through the Schwarz inequality is shown for each element.
other parameters for the G-SS (upper panels) and the 3N-SS for NH (middle panels) and IH (lower panels). As in Fig. 3, the 1 and 2σ regions are colored in red and blue respectively. For the G-SS the strongest limit between the direct bound from radiative LFV decays and the indirect limit from the diagonal entries through the Schwarz inequality is shown. For |η eµ | the constraint from µ → eγ gives the most stringent bound while for |η eτ | and |η µτ | the indirect constraints from the lepton flavour conserving (LFC) processes included in the global fit together with the Schwarz inequality Eq. (33) rather dominate. Moreover, the bound on the product |θ e θ τ | for the 3N-SS shows a 1σ preference for a non-zero value. This mild hint can be translated into a prediction for LFV τ − e transitions, in particular, to a branching ratio of τ → eγ of ∼ 2.5 · 10 −10 for |η eτ | ∼ 6 · 10 −4 . This is rather challenging to probe but not very far from the future sensitivities expected at Super-B factories.

V. DISCUSSION AND CONCLUSIONS
A global fit to lepton flavour and electroweak precision data has been performed to constrain the size presently allowed for the mixing of the extra heavy Seesaw neutrinos with the SM leptons. The analysis has been performed both in a completely general Seesaw (G-SS) with the effects of the extra neutrinos encoded in effective operators with no assumed correlations and for the particular case where only three heavy neutrinos are considered (3N-SS). The results of the fit are summarized in Table IV < 0.026 < 4.9 · 10 −3 < 4.9 · 10 −3 < 4.9 · 10 −3 √ 2η eτ , |θ e θ τ |  For the G-SS with an arbitrary number of extra heavy neutrinos the bounds are expressed in the quantity 2|η αβ | = i Θ αi Θ * βi (see Eq. (13)). Thus, the diagonal elements √ 2η αα correspond to the sum (in quadrature) of all mixings Θ αi of the individual extra heavy neutrinos N i to a given SM flavour α and represent an upper bound on each individual mixing. The off-diagonal entries, on the other hand, are the combinations that can mediate LFV transitions and even provide extra sources of CP-violation. Notice that, from this definition, η is a positive definite matrix and its off-diagonal elements subject to the Schwarz inequality |η αβ | ≤ √ η αα η ββ . In the case of the 3N-SS, only one mixing parameter θ α per SM flavour α can be large enough to saturate the bounds derived here, so as to comply with our present constraints on light neutrino masses and mixings from neutrino oscillation data (see discussion in Section II). Thus, the Schwarz inequality is saturated to an equality for the 3N-SS. Furthermore, some non-trivial correlations between the parameters θ α are also present (see Eq. (12)).
As shown in Table IV the data show a mild, between 1 and 2σ preference, for non-zero heavy-active mixing of order ∼ 0.03 − 0.04 in the e and τ sectors. At the 2σ level, upper bounds in all mixing parameters are found. The most stringent one ∼ 0.02 is found for the mixing with muons, followed by ∼ 0.05 for electrons and ∼ 0.07 for taus. Regarding the off diagonal entries, for the G-SS the indirect bounds from LFC processes can be compared with the direct constraints from LFV observables. Interestingly, the constraint from µ → eγ strongly dominates over all others leading to a bound one order of magnitude better ∼ 0.005 in the e − µ entry, while the e − τ and µ − τ values are rather dominated by the indirect constraints from the Schwarz inequality (comparison between the LFC and LFV columns). Regarding the 3N-SS, even though the necessity of correctly reproducing the observed neutrino mass and mixing pattern introduces non-trivial correlations among the θ α and the neutrino masses and mixings (dependence on normal or inverted hierarchy assumptions shown in the comparison of the third and fourth columns), there is still enough freedom to obtain very similar bounds to those found for the G-SS. This however implies some non-trivial correlations preferred at 1σ notably among the PMNS matrix Majorana phases as well as among the phases of θ e and θ τ as shown in Fig. 4.
The bounds derived here represent the most updated set of constraints and compare well with previous studies. Notably, it is interesting to compare with another recent global fit presented in Ref. [39] were bounds to the G-SS were also studied. We find that the agreement between the two sets of constraints is generally good. The same preference for non-zero mixing in the electron and tau sectors was found but in their case the preferred value is slightly (∼ 20 − 30%) larger. Similarly the upper bound on muon mixing is weaker in Ref. [39]. Conversely the limits on the off-diagonal elements are slightly (∼ 20 − 40%) stronger in Ref. [39] for the e − τ and µ − τ sectors. The only very noticeable difference is in the e − µ sector where the limit from µ → eγ is almost a factor 3 stronger than the one presented here (despite not being yet updated to the final MEG result). This difference can be attributed to not considering the propagation of the heavy neutrinos in the loop for the process which tends to restore the GIM cancellation (given the Unitarity of the full mixing matrix) and to therefore slightly weaken by the corresponding factor the bound stemming from the process. This extra contribution was not taken into account in Ref. [39] since a more agnostic source of the non-unitarity of the PMNS matrix was adopted while here we concentrate in constraining heavy neutrino mixings. The rest of the discrepancies can stem from small differences in our analyses. For example our observables for weak lepton universality and CKM unitarity are more updated and our bounds correspond to frequentist confidence regions while Ref. [39] rather presented Bayesian credible intervals. Regarding the 3N-SS, the closest study of a similar setup in the literature is that of Ref. [40]. This work is rather complementary to our results focusing instead in the region between 10 to 250 GeV, where more stringent constraints are derived since the extra neutrinos would be kinematically accessible.
It is also interesting to translate the bounds derived here to other common parametrizations, useful in particular for the analysis of neutrino non-standard interactions (see e.g. Ref. [88]). Indeed, the non-unitary PMNS matrix induced by the mixing with the extra heavy neutrinos modifies the neutrino production and detection processes, which can be encoded in production/detection NSI [32,45]. In particular: Furthermore, neutrino interactions with matter are also affected and these effects can also be described by matter NSI [32]: where n e and n n are the electron and neutron densities of the matter traversed by the neutrinos.
Finally, an alternative parametrization of the non-unitarity of the PMNS matrix of the form N = T U with T a lower triangular matrix [89][90][91]: is also considered appropriate to study the effects of non-unitary PMNS mixing in neutrino oscillation searches [45,[92][93][94][95]. Comparing Eqs. (3) and (45) it is easy to see that α ββ ≈ 1−η ββ , while |α βγ | ≈ 2|η βγ | = | βγ | so that the bounds derived here can be trivially translated to this parametrization too. All in all this level of non-unitarity (or equivalently NSI as in Eq. (43)) is extremely tough to probe at present or near-future neutrino oscillations facilities and its effects would be negligible. However, prospective very precise neutrino oscillation facilities such as the Neutrino factory [96,97] could probe beyond this very stringent present limits for some elements [45,92].
Notice that the bounds derived here apply for any heavy neutrino mass above the electroweak scale. For lighter heavy neutrino masses, the LFV radiative decays start to be suppressed by the restoration of the GIM mechanism (see Eq. (42)) and therefore the constraints shown in the LFV column of Table IV are not valid. The rest of the bounds summarized in the LFC column of Table IV do apply down to O(500 MeV) with the only exception of the invisible width of the Z, since for masses below ∼ M Z /2 the heavy neutrinos can be kinematically produced and unitarity is restored. Therefore, in the region between the Kaon mass and the EW scale we do not expect any significant change in the G-SS bounds shown in the LFC column of Table IV. Nevertheless, at these lower energies were the extra neutrinos can be directly produced, more stringent constraints than the ones derived here, from direct searches [93,[98][99][100][101][102] and cosmology [103][104][105][106][107][108][109][110][111][112][113][114][115] apply.
In summary, we have combined present probes on weak lepton universality, searches for LFV processes and precision electroweak observables to derive updated and global constraints on the allowed mixing of heavy Seesaw neutrinos with the SM fermions. These bounds apply for any value of the Majorana scale larger than the electroweak scale and have been computed both for a completely general scenario as well as for the case in which only 3 extra heavy neutrinos are considered. At the 1σ level a mild preference for non-zero mixing in the electron and tau sectors around 0.03 − 0.04 was found, which could be probed for by improving the LFC searches that currently lead to that preference, as well as through τ − e LFV transitions. At the 2σ level, upper bounds between 10 −1 and 10 −2 for all elements were derived with a most stringent constraint on the mixing in the e − µ sector an order of magnitude better from the µ → eγ process. While this is by far the present dominant bound, it will be superseded in the future by µ → eee and/or µ − e conversion in nuclei searches. Apart from this and other improvements in the datasets considered, this level of mixing is challenging but still plausible to probe at future collider [40,70,116,117] and dedicated neutrino oscillation searches [45,92].