IFT-UAM/CSIC-23-60

FTUV-23-0531.7594

IFIC/23-19

Bounds on lepton non-unitarity and heavy neutrino mixing

Mattias<EMAIL_ADDRESS>Enrique Fernández-
<EMAIL_ADDRESS>Josu Hernández-
<EMAIL_ADDRESS>
Jacobo<EMAIL_ADDRESS>Xabier
<EMAIL_ADDRESS>and Daniel Naredo-
<EMAIL_ADDRESS>
1Department of Physics, School of Engineering Sciences, KTH Royal Institute of
Technology,
AlbaNova University Center, Roslagstullsbacken 21, SE-106 91 Stockholm, Sweden
2Departamento de Física Teórica and Instituto de Física Teórica UAM/CSIC,
Universidad Autónoma de Madrid, Cantoblanco, 28049 Madrid, Spain
3Institute for Theoretical Physics, ELTE Eötvös Loránd University,
Pázmány Péter sétány 1/A, H-1117 Budapest, Hungary
4 Instituto de Física Corpuscular, Universidad de Valencia and CSIC
Edificio Institutos Investigación, Catedrático José Beltrán 2, 46980 Spain

We present an updated and improved global fit analysis of current flavor and
electroweak precision observables to derive bounds on unitarity deviations of
the leptonic mixing matrix and on the mixing of heavy neutrinos with the
active flavours. This new analysis is motivated by new and updated
experimental results on key observables such as $V_{ud}$, the invisible decay
width of the $Z$ boson and the $W$ boson mass. It also improves upon previous
studies by considering the full correlations among the different observables
and explicitly calibrating the test statistic, which may present significant
deviations from a $\chi^{2}$ distribution. The results are provided for three
different Type-I seesaw scenarios: the minimal scenario with only two
additional right-handed neutrinos, the next to minimal one with three extra
neutrinos, and the most general one with an arbitrary number of heavy
neutrinos that we parametrize via a generic deviation from a unitary leptonic
mixing matrix. Additionally, we also analyze the case of generic deviations
from unitarity of the leptonic mixing matrix, not necessarily induced by the
presence of additional neutrinos. This last case relaxes some correlations
among the parameters and is able to provide a better fit to the data.
Nevertheless, inducing only leptonic unitarity deviations avoiding both the
correlations implied by the right-handed neutrino extension as well as more
strongly constrained operators is challenging and would imply significantly
more complex UV completions.

###### Contents

  1. 1 Introduction
  2. 2 Observables
    1. 2.1 Constraints from $M_{W}$ and $s^{2}_{\text{eff}}$
    2. 2.2 Constraints from Z-pole observables
    3. 2.3 Constraints from lepton flavor universality (LFU)
    4. 2.4 Constraints from CKM unitarity
    5. 2.5 About the CDF-II $M_{W}$ measurement
    6. 2.6 Charged Lepton Flavor Violation
  3. 3 Global fit bounds for the 2 neutrino case (2N-SS)
  4. 4 Global fit bounds for the 3 neutrino case (3N-SS)
  5. 5 Global fit bounds for the general neutrino case (G-SS)
  6. 6 Global fit bounds for generic unitarity violation (GUV)
  7. 7 Discussion and Conclusions
  8. A Details of the analysis
  9. B About cancellations in the cLFV rates

##  1 Introduction

Adding right-handed (RH) neutrinos to the Standard Model (SM) particle
spectrum is one of the best motivated extensions to address several of the
open problems in particle physics. They constitute the simplest possibility to
accommodate the evidence for neutrino masses and mixing, and they can also
play an important role as dark matter candidates [1, 2, 3, 4] or portals to
dark sectors [5, 6, 7, 8] as well as in the generation of the matter-
antimatter asymmetry of the Universe [9, 10, 11, 12].

Due to the singlet nature of the RH neutrinos, also known as heavy neutral
leptons (HNLs), there is a new energy scale associated to their allowed
Majorana mass, which is unknown from the theoretical point of view. If this
scale is low enough, they could be produced at different laboratories,
motivating the strong experimental effort (see for instance [13]) that led to
a plethora of constraints (see e.g. [14, 15]). On the other hand, if the scale
is above the experimental energy, the RH neutrinos cannot be produced, but
their existence induces deviations from unitarity of the leptonic mixing
matrix, so they can still be probed for at the intensity frontier via
electroweak precision observables, universality ratios or charged lepton
flavor violating (cLFV) processes [16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26,
27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45,
46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56].

When the RH neutrino masses are above the experimental energy, it is useful to
consider an effective field theory (EFT) parametrization of their impact at
accessible energies. After integrating them out, the lowest order dim-5
Weinberg operator [57] generates neutrino masses and mixings after electroweak
symmetry breaking, while the only dim-6 induced at tree level leads to non-
unitarity effects in lepton mixing [58]. In a high-scale seesaw mechanism [59,
60, 61, 62], the naive expectation is that both dim-5 and 6 operators are
suppressed, explaining the smallness of neutrino masses at the price of
rendering the model virtually untestable. On the other hand, symmetry
protected [63, 64, 65, 66] low-scale seesaw mechanisms, such as the linear
[67, 68] or inverse [69, 70] seesaw mechanisms, suppress the dim-5 operator,
but not the dim-6 one, allowing for large deviations from unitarity.

When extending the SM with $n$ new RH neutrinos $N_{R}^{i}$, both a neutrino
Yukawa and a Majorana mass term for the right-handed neutrinos should be added
to the Lagrangian:

$\mathcal{L}=-\frac{1}{2}\overline{N_{R}^{i\,c}}M_{M}^{ij}N_{R}^{j}-Y_{\nu}^{\alpha
i}\overline{L_{\alpha}}\tilde{\phi}N^{i}_{R}+h.c.$ (1)

where $\tilde{\phi}=i\sigma_{2}\phi$ stands for the SM Higgs doublet,
$Y_{\nu}$ is a complex $3\times n$ Yukawa coupling matrix and $M_{M}$ is a
$n\times n$ symmetric Majorana matrix for the $N_{R}$ fields. When the RH
neutrinos are heavy, upon integrating them out, they induce light neutrino
masses via the dim-5 Weinberg operator [57],

$m_{\rm\nu}=-\Theta M_{M}\Theta^{T}\,,$ (2)

where $\Theta$ is the mixing between the active SM neutrinos and the RH ones:

$\Theta\equiv\dfrac{v}{\sqrt{2}}Y_{\nu}M_{M}^{-1}=m_{D}M_{M}^{-1}\,,$ (3)

and $v=246$ GeV is the Higgs vev. The symmetric mass matrix for the light
neutrinos $m_{\nu}$ will be diagonalized by a unitary matrix $U$.
Additionally, as already mentioned, the only dimension 6 operator generated at
tree level induces the non-unitarity effects we are interested in [58]. In
particular, the leptonic mixing matrix coupling the light neutrino mass
eigenstates to the SM charged leptons through CC interactions will be given
by:

$N=(\mathbb{I}-\eta)\,U\,.$ (4)

We have dubbed the leptonic mixing matrix $N$ to emphasize its _non-unitary_
character. Indeed, while $U$ is the unitary matrix that diagonalizes the
Weinberg $d=5$ operator in Eq. (2), $\eta$ parametrizes the unitarity
deviations of $N$ and corresponds to (half of) the coefficient of the $d=6$
operator obtained upon integrating out the heavy fields [58]:

$\eta=\dfrac{1}{2}\Theta\Theta^{\dagger}\,.$ (5)

Notice that, any general matrix can be parametrized through the product of an
Hermitian and a unitary matrix. Hence, Eq. (4) is a completely general and
convenient way to encode unitarity deviations through the small Hermitian
matrix $\eta$ [71]. As we will show in more detail later, most observables
depend on the combination $NN^{\dagger}$, since the flavor indices
corresponding to the charged leptons are fixed by the process and the physical
neutrino indices are summed over, since they are not measured. Being $\eta$ a
small parameter, we have

$\sum_{i=1}^{3}N^{\phantom{\dagger}}_{\alpha
i}N^{\dagger}_{i\beta}=\delta_{\alpha\beta}-2\eta_{\alpha\beta}+\mathcal{O}\big{(}\eta_{\alpha\beta}^{2}\big{)}\,.$
(6)

Thus, the observables do not depend on $U$ and the $\eta$ matrix encodes any
possible effect arising from non-unitarity, regardless of the UV completion
that originates the deviations and how many fields it contains.

In all generality, the dim-5 and 6 operators are independent, meaning that in
the most general case we can decouple $\eta$ from light neutrino masses and
mixing. Nevertheless, when induced by RH neutrinos, $\eta$ is a positive-
definite matrix, as shown by Eq. (5), and hence its parameters are not
entirely independent. In particular, the Schwarz inequality will be satisfied:

$\left\lvert\eta_{\alpha\beta}\right\rvert\leq\sqrt{\eta_{\alpha\alpha}\eta_{\beta\beta}}\,.$
(7)

Furthermore, since both operators are generated from different combinations of
the same Yukawa and mass matrices (see Eqs. (2) and (5)), it is possible to
find correlations between them depending on how many new heavy neutrinos are
considered. Thus, with the aim of covering from the minimal to the most
general case, we will consider the following 4 scenarios:

  * •

The minimal scenario that accommodates oscillation data adding only 2 RH
neutrinos. Here the dim-5 and 6 operators are fully correlated, since the
latter (and thus $\eta$) can be fully reconstructed from the former up to a
global scale [72].

  * •

The next to minimal scenario with 3 RH neutrinos, where the dim-5 operator
still imposes strong correlations on $\eta$ [73].

  * •

The most general RH neutrino scenario where we assume a general $\eta$ matrix,
independent of neutrino oscillation data, which encodes the low-energy effects
of an arbitrary number of heavy RH neutrinos.

  * •

Generic deviations from unitarity not requiring $\eta$ to be positive-definite
or subject to the Schwarz inequality. This would imply more elaborate
additions to the SM particle content beyond only RH neutrinos (see Ref. [55]
for a dedicated discussion).

In order to obtain a complete picture of the experimental situation, we will
perform a global fit analysis of an extended set of observables that are
affected by a non-unitary leptonic mixing matrix, and derive current limits
for all these scenarios. Compared to previous analyses [44, 74, 54] we improve
and extend the study in several ways:

  * •

We update and complete the list of observables taking into account all
correlations among them consistently. We also comment on the role of the CDF-
II anomalous measurement of the $W$ mass [75], given the importance of this
observable in the non-unitarity analysis.

  * •

We extend the number of scenarios considered. In particular and as mentioned
above, we specifically consider the minimal case with two heavy neutrinos and
obtain its corresponding bounds, which is an analysis that was previously
missing in the literature.

  * •

We improve the statistical treatment with an explicit calibration of the test
statistic, since deviations from Wilks’ theorem are expected for most of the
scenarios under consideration.

  * •

We avoid using the values provided in the PDG review for the SM prediction of
the different observables. Indeed, while this procedure is common and
convenient, these values are obtained from a global fit which includes also
the observables of interest in this study. This approach is therefore not
consistent as it double counts these observables and, by using them in the
prediction, it may artificially reduce any preference for new physics. We
instead choose an independent and accurate set of input observables distinct
from the others to make our predictions consistently.

This paper is organized as follows. In section 2 we describe the observables
considered in our global fits specifying their dependence on the unitarity
deviations of the leptonic mixing matrix and how we include them in our
statistical analysis. In sections 3, 4 and 5 we present the results for the
scenarios with 2, 3 and a general number of additional RH neutrinos
respectively, while in section 6 we show the impact of dropping the assumption
of $\eta$ being positive definite as it might be the case if not induced
(only) by mixing with RH neutrinos. Finally, we summarize our conclusions in
section 7.

##  2 Observables

In this work, we will derive global bounds on deviations from unitarity, and
consequently on heavy neutrino mixing, via a fit to the following observables:

  * $\bullet$

Four determinations of the W-boson mass: $M_{W}^{\text{LEP}}$,
$M_{W}^{\text{Tev}}$, $M_{W}^{\text{LHCb}}$, $M_{W}^{\text{ATLAS}}$.

  * $\bullet$

Two determinations of the effective weak angle: $s^{2\hskip
4.26773pt\text{LHC}}_{\text{eff}}$ and $s^{2\hskip
4.26773pt\text{Tev}}_{\text{eff}}$.

  * $\bullet$

Five LEP observables measured at the Z-pole, plus a determination of the $Z$
invisible width from CMS: $\Gamma_{Z}$, $\sigma^{0}_{\text{had}}$, $R_{e}$,
$R_{\mu}$, $R_{\tau}$, $\Gamma_{\text{inv}}^{\text{LHC}}$.

  * $\bullet$

Five weak decay ratios constraining lepton flavor universality: $R_{\mu
e}^{\pi}$, $R_{\tau\mu}^{\pi}$, $R_{\mu e}^{K}$, $R_{\mu e}^{\tau}$,
$R_{\tau\mu}^{\tau}$.

  * $\bullet$

Ten weak decays constraining CKM unitarity.

  * $\bullet$

cLFV observables.

With the obvious exception of cLFV processes, all these observables are lepton
flavor conserving (LFC) and thus constrain the diagonal entries of the $\eta$
matrix. Nevertheless, through the Schwarz inequality of Eq. (7), LFC
observables will also impose constraints on the off-diagonal entries. As we
will see, these limits will be often stronger than those imposed by direct
contributions to cLFV processes.

### SM input parameters and predictions

For the SM input parameters the very precise determinations of the Z-boson
mass $M_{Z}$, the fine structure constant $\alpha$ and the Fermi constant
($G_{F}$) extracted from $\mu$ decay $G_{\mu}$ will be used. The corresponding
values are taken from the PDG [76]:

$\displaystyle M_{Z}$ $\displaystyle=91.1876(21)\text{ GeV}\,,$ (8)
$\displaystyle\alpha^{-1}$ $\displaystyle=137.035999180(10)\,,$ $\displaystyle
G_{\mu}$ $\displaystyle=1.1663788(6)\cdot 10^{-5}\text{ GeV${}^{-2}$}\,.$

However, since the extraction of $G_{F}$ comes from a charged current weak
decay and the $W$ boson couplings to leptons are modified due to the non-
unitarity of the leptonic mixing matrix, $G_{\mu}$ will no longer be equal to
$G_{F}$. In fact, in the presence of a non-unitary mixing matrix the muon
decay width gets modified as:

$\Gamma_{\mu}=\frac{G_{F}^{2}m_{\mu}^{5}}{192\pi^{3}}\sum_{i=1}^{3}\left\lvert
N_{\mu i}\right\rvert^{2}\sum_{j=1}^{3}\left\lvert
N_{ej}\right\rvert^{2}\simeq\frac{G_{F}^{2}m_{\mu}^{5}}{192\pi^{3}}\left(1-2\eta_{ee}-2\eta_{\mu\mu}\right)\equiv\frac{G_{\mu}^{2}m_{\mu}^{5}}{192\pi^{3}}\,,$
(9)

and thus $G_{F}$, as measured through muon decay, picks a non-unitarity
dependence,

$G_{F}\simeq G_{\mu}\left(1+\eta_{ee}+\eta_{\mu\mu}\right),$ (10)

which will propagate into any observable which depends on it. In the
following, we discuss in detail the dependence of these observables on the
non-unitarity parameters. The main results are already summarised in Table 1.

The SM predictions for all the observables entering in our fit will also be
necessary. A common practice is to take their SM values as given by the PDG,
however this would not be consistent since the PDG values correspond to a
global fit including the electroweak observables themselves. Hence, since our
goal is precisely to perform such a fit but in the context of non-unitarity,
it is not consistent to take values which have already been fitted. Moreover,
since the experimental data is being fitted, the resulting SM predictions
would be systematically closer to the experimental values than the values
obtained by just inserting the SM input parameters into the expression for the
observables. This may artificially diminish any preference for new physics in
the data.

For this reason, we will instead derive the SM predictions for each observable
from its parametric dependence on the SM input parameters. This implies
including the relevant loop corrections, which requires additional inputs: the
mass of the top quark ($m_{t}$), the mass of the Higgs boson ($M_{H}$), the
strong coupling constant at $M_{Z}$ ($\alpha_{s}\left(M_{Z}\right)$) and the
running of the fine-structure constant at $M_{Z}$
($\Delta\alpha\left(M_{Z}\right)$). We adopt the following values [76]:

$\displaystyle m_{t}$ $\displaystyle=172.69(30)\text{ GeV}\,,$ (11)
$\displaystyle M_{H}$ $\displaystyle=125.25(17)\text{ GeV}\,,$
$\displaystyle\alpha_{s}\left(M_{Z}\right)$ $\displaystyle=0.1185(16)\,,$
$\displaystyle\Delta\alpha\left(M_{Z}\right)$ $\displaystyle=591.05(70)\cdot
10^{-5}\,.$

Notice that we only consider SM loop contributions, which are of course
necessary for the precise evaluation of the SM prediction of each observable.
In principle one should also take into account heavy neutrino loop
contributions to the non-unitarity dependence of the observables, however it
has been shown that these can be safely neglected as they are subleading [73].

Observable | SM prediction | Experimental value
---|---|---
$M_{W}\simeq M_{W}^{\text{SM}}\left(1+0.20\left(\eta_{ee}+\eta_{\mu\mu}\right)\right)$ |  $80.356(6)$ GeV |  $80.373(11)$ GeV | -
$s^{2\hskip 4.26773pt\text{Tev}}_{\text{eff}}\simeq s^{2\hskip 4.26773pt\text{SM}}_{\text{eff}}\left(1-1.40\left(\eta_{ee}+\eta_{\mu\mu}\right)\right)$ | $0.23154(4)$ | $0.23148(33)$ | [76]
$s^{2\hskip 4.26773pt\text{LHC}}_{\text{eff}}\simeq s^{2\hskip 4.26773pt\text{SM}}_{\text{eff}}\left(1-1.40\left(\eta_{ee}+\eta_{\mu\mu}\right)\right)$ | $0.23154(4)$ | $0.23129(33)$ | [76]
$\Gamma_{\text{inv}}^{\text{LHC}}\simeq\Gamma_{\text{inv}}^{\text{SM}}\left(1-0.33\left(\eta_{ee}+\eta_{\mu\mu}\right)-1.33\eta_{\tau\tau}\right)$ |  $0.50145(5)$ GeV |  $0.523(16)$ GeV | [77]
$\Gamma_{Z}\simeq\Gamma_{Z}^{\text{SM}}\left(1+1.08\left(\eta_{ee}+\eta_{\mu\mu}\right)-0.27\eta_{\tau\tau}\right)$ |  $2.4939(9)$ GeV |  $2.4955(23)$ GeV | [76]
$\sigma^{0}_{\text{had}}\simeq\sigma^{0\hskip 4.26773pt\text{SM}}_{\text{had}}\left(1+0.50\left(\eta_{ee}+\eta_{\mu\mu}\right)+0.53\eta_{\tau\tau}\right)$ |  $41.485(8)$ nb |  $41.481(33)$ nb | [76]
$R_{e}\simeq R_{e}^{\text{SM}}\left(1+0.27\left(\eta_{ee}+\eta_{\mu\mu}\right)\right)$ | $20.733(10)$ | $20.804(50)$ | [76]
$R_{\mu}\simeq R_{\mu}^{\text{SM}}\left(1+0.27\left(\eta_{ee}+\eta_{\mu\mu}\right)\right)$ | $20.733(10)$ | $20.784(34)$ | [76]
$R_{\tau}\simeq R_{\tau}^{\text{SM}}\left(1+0.27\left(\eta_{ee}+\eta_{\mu\mu}\right)\right)$ | $20.780(10)$ | $20.764(45)$ | [76]
$R_{\mu e}^{\pi}\simeq\left(1-\left(\eta_{\mu\mu}-\eta_{ee}\right)\right)$ | $1$ | $1.0010(9)$ | [78]
$R_{\tau\mu}^{\pi}\simeq\left(1-\left(\eta_{\tau\tau}-\eta_{\mu\mu}\right)\right)$ | $1$ | $0.9964(38)$ | [78]
$R_{\mu e}^{K}\simeq\left(1-\left(\eta_{\mu\mu}-\eta_{ee}\right)\right)$ | $1$ | $0.9978(18)$ | [78]
$R_{\mu e}^{\tau}\simeq\left(1-\left(\eta_{\mu\mu}-\eta_{ee}\right)\right)$ | $1$ | $1.0018(14)$ | [78]
$R_{\tau\mu}^{\tau}\simeq\left(1-\left(\eta_{\tau\tau}-\eta_{\mu\mu}\right)\right)$ | $1$ | $1.0010(14)$ | [78]
$\left\lvert V_{ud}^{\beta}\right\rvert\simeq\sqrt{1-\left\lvert V_{us}\right\rvert^{2}}\left(1+\eta_{\mu\mu}\right)$ | $\sqrt{1-\left\lvert V_{us}\right\rvert^{2}}$ | $0.97373(31)$ | [76]
$\left\lvert V_{us}^{\tau\rightarrow K\nu}\right\rvert\simeq\left\lvert V_{us}\right\rvert\left(1+\eta_{ee}+\eta_{\mu\mu}-\eta_{\tau\tau}\right)$ | $\left\lvert V_{us}\right\rvert$ | $0.2236(15)$ | [79]
$\left\lvert V_{us}^{\tau\rightarrow K,\pi}\right\rvert\simeq\left\lvert V_{us}\right\rvert\left(1+\eta_{\mu\mu}\right)$ | $\left\lvert V_{us}\right\rvert$ | $0.2234(15)$ | [76]
$\left\lvert V_{us}^{K_{L}\rightarrow\pi e\nu}\right\rvert\simeq\left\lvert V_{us}\right\rvert\left(1+\eta_{\mu\mu}\right)$ | $\left\lvert V_{us}\right\rvert$ | $0.2229(6)$ | [76]
$\left\lvert V_{us}^{K_{L}\rightarrow\pi\mu\nu}\right\rvert\simeq\left\lvert V_{us}\right\rvert\left(1+\eta_{ee}\right)$ | $\left\lvert V_{us}\right\rvert$ | $0.2234(7)$ | [76]
$\left\lvert V_{us}^{K_{S}\rightarrow\pi e\nu}\right\rvert\simeq\left\lvert V_{us}\right\rvert\left(1+\eta_{\mu\mu}\right)$ | $\left\lvert V_{us}\right\rvert$ | $0.2220(13)$ | [76]
$\left\lvert V_{us}^{K_{S}\rightarrow\pi\mu\nu}\right\rvert\simeq\left\lvert V_{us}\right\rvert\left(1+\eta_{ee}\right)$ | $\left\lvert V_{us}\right\rvert$ | $0.2193(48)$ | [76]
$\left\lvert V_{us}^{K^{\pm}\rightarrow\pi e\nu}\right\rvert\simeq\left\lvert V_{us}\right\rvert\left(1+\eta_{\mu\mu}\right)$ | $\left\lvert V_{us}\right\rvert$ | $0.2239(10)$ | [76]
$\left\lvert V_{us}^{K^{\pm}\rightarrow\pi\mu\nu}\right\rvert\simeq\left\lvert V_{us}\right\rvert\left(1+\eta_{ee}\right)$ | $\left\lvert V_{us}\right\rvert$ | $0.2238(12)$ | [76]
$\left\lvert\dfrac{V_{us}}{V_{ud}}\right\rvert^{K,\pi\rightarrow\mu\nu}\simeq\dfrac{\left\lvert V_{us}\right\rvert}{\sqrt{1-\left\lvert V_{us}\right\rvert^{2}}}$ | $\dfrac{\left\lvert V_{us}\right\rvert}{\sqrt{1-\left\lvert V_{us}\right\rvert^{2}}}$ | $0.23131(53)$ | [76]
Table 1:  Set of precision observables used as input for our global fit. The
first column includes their dependence on the non-unitarity parameters, the
second their SM prediction and the third one their experimental value. For
$M_{W}$, we compute the current world average without including the CDF-II
measurement (see text for details).

###  2.1 Constraints from $\boldsymbol{M_{W}}$ and
$\boldsymbol{s^{2}_{\text{eff}}}$

Although $M_{W}$ and $s_{\mathrm{w}}^{2}$ do not depend directly on the non-
unitarity parameters, they provide alternative determinations of $G_{F}$ to be
compared to the one extracted from muon decay, which, as previously discussed,
does depend on $\eta$. In other words, the predicted values for $M_{W}$ and
$s_{\mathrm{w}}^{2}$ through $G_{\mu}$ will inherit the following dependence
on $\eta$

$\displaystyle M_{W}$
$\displaystyle=M_{Z}\sqrt{\dfrac{1}{2}+\sqrt{\dfrac{1}{4}-\dfrac{\pi\alpha\left(1+\eta_{ee}+\eta_{\mu\mu}\right)}{\sqrt{2}G_{\mu}M_{Z}^{2}\left(1-\Delta
r\right)}}}\,,$ (12) $\displaystyle s_{\mathrm{w}}^{2}$
$\displaystyle=\dfrac{1}{2}\left(1-\sqrt{1-\dfrac{2\sqrt{2}\pi\alpha\left(1+\eta_{ee}+\eta_{\mu\mu}\right)}{G_{\mu}M_{Z}^{2}\left(1-\Delta
r\right)}}\right)\,,$ (13)

with the SM radiative corrections included in $\Delta r=0.03657(22)$ [76].
Here, it is useful to factorize the SM contribution and keep only the leading
corrections in $\eta$, obtaining the expressions reported in Table 1. Then, we
can compute more precisely the SM predictions for $M_{W}$ and
$s_{\text{eff}}^{2}$ using, respectively, the formulas given in [80] and [81].
Notice that the “on-shell” value of $s_{\mathrm{w}}^{2}$ is related to the
effective weak angle $s_{\text{eff}}^{2}$ in the table through a
multiplicative factor $\kappa$ that includes SM loops [76]

$s_{\text{eff}}^{2}=\kappa\,s_{\mathrm{w}}^{2}\,,$ (14)

and, therefore, it has the same $\eta$-dependence than $s_{\mathrm{w}}^{2}$.

These $\eta$-dependent predictions are to be compared with the current
experimental values. For the $W$ boson mass, we compute the current world
average also considering the most recent ATLAS result [82]. For this, we
follow the PDG [76] prescription and combine the results of ATLAS, LHCb [83],
Tevatron [84] and LEP [85] following the BLUE procedure [86, 87] and taking
into account the corresponding correlations111Following the PDG [76], we will
assume a correlated uncertainty of 9 MeV for LHCb and ATLAS in order to
compute the LHC average. Then, assuming a correlation of 7 MeV between LHC and
Tevatron, we compute the hadron collider average. The world average is then
extracted by combining the hadron collider average with the LEP measurement
assuming no correlations.. In section 2.5 we will discuss the latest and
anomalous CDF-II measurement [75]. Since we find there is too much tension
with other observables even in presence of non-unitarity, we will not include
it in the global fit.

Regarding the weak mixing angle, we use two measurements of
$s^{2}_{\text{eff}}$ coming from the LHC and Tevatron. On the other hand, we
do not include low-energy determinations of $s_{\mathrm{w}}^{2}$, such as its
determination from Møller scattering. This is due to the fact that in the
electron vector coupling $g_{V}^{e}=-\frac{1}{2}+2s_{\mathrm{w}}^{2}$ there is
a partial cancellation (since $s_{\mathrm{w}}^{2}\sim 1/4$). This renders the
extraction of $s_{\mathrm{w}}^{2}$ very sensitive to higher order
contributions which would require a higher level of accuracy to include it
reliably in the fit.

###  2.2 Constraints from Z-pole observables

All of the predicitions for the Z-pole observables will be modified, as they
depend on $G_{F}$ and $s_{\mathrm{w}}^{2}$. For instance, this is the case for
the partial decay widths to charged fermions:

$\Gamma_{f}=N_{C}\dfrac{G_{\mu}M_{Z}^{3}\left(g_{V}^{f2}+g_{A}^{f2}\right)}{6\sqrt{2}\pi}\left(1+\eta_{ee}+\eta_{\mu\mu}\right),$
(15)

where $N_{C}$ is a color factor and $g_{V}^{f}$ and $g_{A}^{f}$ are the vector
and axial couplings of the $Z$ boson to fermions, respectively:

$\displaystyle g_{V}^{f}$ $\displaystyle=T_{f}-2Q_{f}s_{\mathrm{w}}^{2}\,,$
(16) $\displaystyle g_{A}^{f}$ $\displaystyle=T_{f}\,,$

with $T_{f}$ and $Q_{f}$ being the isospin and charge of the fermion.

Moreover, the $Z$ boson invisible decay width will be also modified when the
right-handed neutrinos are heavier than $M_{Z}$ and not kinematically
available:

$\Gamma_{\text{inv}}=\dfrac{G_{\mu}M_{Z}^{3}\sum_{i,j=1}^{3}\Big{|}\big{(}N^{\dagger}N\big{)}_{ij}\Big{|}^{2}}{12\sqrt{2}\pi}\left(1+\eta_{ee}+\eta_{\mu\mu}\right)\simeq\dfrac{G_{\mu}M_{Z}^{3}}{12\sqrt{2}\pi}\Big{(}3-\big{(}\eta_{ee}+\eta_{\mu\mu}+4\eta_{\tau\tau}\big{)}\Big{)}\,.$
(17)

Thus, non-unitarity effects reduce the $Z$ invisible decay width with respect
to the SM prediction which, interestingly, helped explaining the long-standing
$\sim 2\sigma$ LEP anomaly in the number of neutrinos. Nevertheless, this
anomaly is now gone after an improvement in the computation of the Bhabha
scattering cross section [88]. On top of that, a recent determination of
$\Gamma_{\text{inv}}$ from the CMS experiment [77] is $\sim 1\sigma$ away from
the SM towards larger values of $\Gamma_{\text{inv}}$, which imposes further
constraints on non-unitarity effects.

Having this in mind, the observables of our fit will be the usual combinations
used for LEP from the partial widths of the Z. Namely, we have the total width
$\Gamma_{Z}$, the hadronic cross section $\sigma^{0}_{\text{had}}$ and the
leptonic ratios $R_{\ell}$, defined as:

$\displaystyle\Gamma_{Z}$ $\displaystyle=\sum_{f\neq
t}\Gamma_{f}+\Gamma_{\text{inv}}\,,$ (18)
$\displaystyle\sigma^{0}_{\text{had}}$
$\displaystyle=\dfrac{12\pi\Gamma_{e}\Gamma_{\text{had}}}{M_{Z}^{2}\Gamma_{Z}^{2}}\,,$
(19) $\displaystyle R_{\ell}$
$\displaystyle=\dfrac{\Gamma_{\ell}}{\Gamma_{\text{had}}}\,,$ (20)

where $\Gamma_{\text{had}}=\sum_{q\neq t}\Gamma_{q}$. Moreover, we will
include the new CMS measurement of $\Gamma_{\rm inv}$ as an additional
independent observable.

The parametric dependences of the different $Z$ partial widths in order to
extract the SM predictions are taken from [89]. The expression for the singlet
vector contribution to the full hadronic width of the $Z$ is taken from [90].
It should be noted that, since all of these observables (except
$\Gamma_{\text{inv}}^{\text{LHC}}$) were measured in LEP, correlations between
them must be taken into account, which we have taken from Table B.13 of [88].

Finally, it is worth mentioning that there exist other LEP observables, such
as $R_{c}$, $R_{b}$ or the forward-backward asymmetries $A_{\rm FB}$, which
are not included here. The reason why we have chosen to leave them out is the
fact that they are not as precisely measured as other observables that depend
on the same combination of $\eta$-parameters, thus they hold less constraining
power and would end up diluting the goodness-of-fit by artificially increasing
the number of degrees of freedom.

###  2.3 Constraints from lepton flavor universality (LFU)

The tightest constraints on the universality of weak interactions among the
different lepton flavours come from ratios of meson or charged lepton decay
widths differing in the flavor of one of the leptons involved. Since they are
ratios, the $\eta$ dependence coming from $G_{F}$ will cancel out, and only
the contributions coming from the $W$ vertex will matter. Namely:

$\displaystyle\dfrac{\Gamma\left(P\rightarrow\mu\nu\right)}{\Gamma\left(P\rightarrow
e\nu\right)}$
$\displaystyle=\dfrac{\Gamma\left(P\rightarrow\mu\nu_{\mu}\right)^{\text{SM}}}{\Gamma\left(P\rightarrow
e\nu_{e}\right)^{\text{SM}}}\dfrac{\sum_{i=1}^{3}\left\lvert N_{\mu
i}\right\rvert^{2}}{\sum_{i=1}^{3}\left\lvert
N_{ei}\right\rvert^{2}}\simeq\dfrac{\Gamma\left(P\rightarrow\mu\nu_{\mu}\right)^{\text{SM}}}{\Gamma\left(P\rightarrow
e\nu_{e}\right)^{\text{SM}}}\left(R_{\mu e}^{P}\right)^{2},$ (21)
$\displaystyle\dfrac{\Gamma\left(\tau\rightarrow
P\nu\right)}{\Gamma\left(P\rightarrow\mu\nu\right)}$
$\displaystyle=\dfrac{\Gamma\left(\tau\rightarrow
P\nu_{\tau}\right)^{\text{SM}}}{\Gamma\left(P\rightarrow\mu\nu_{\mu}\right)^{\text{SM}}}\dfrac{\sum_{i=1}^{3}\left\lvert
N_{\tau i}\right\rvert^{2}}{\sum_{i=1}^{3}\left\lvert N_{\mu
i}\right\rvert^{2}}\simeq\dfrac{\Gamma\left(\tau\rightarrow
P\nu_{\tau}\right)^{\text{SM}}}{\Gamma\left(P\rightarrow\mu\nu_{\mu}\right)^{\text{SM}}}\left(R_{\tau\mu}^{P}\right)^{2},$
(22)

where $P=\pi,K$ and the ratio $R^{P}_{\alpha\beta}$ is defined as:

$R_{\alpha\beta}^{P}=1-\left(\eta_{\alpha\alpha}-\eta_{\beta\beta}\right).$
(23)

Similarly, there are also competitive bounds coming from ratios of fully
leptonic decays. In particular, the $\mu-e$ sector can be constrained via:

$\dfrac{\Gamma\left(\tau\rightarrow\mu\nu\nu\right)}{\Gamma\left(\tau\rightarrow
e\nu\nu\right)}=\dfrac{\Gamma\left(\tau\rightarrow\mu\nu_{\mu}\nu_{\tau}\right)^{\text{SM}}}{\Gamma\left(\tau\rightarrow
e\nu_{e}\nu_{\tau}\right)^{\text{SM}}}\dfrac{\sum_{i=1}^{3}\left\lvert N_{\mu
i}\right\rvert^{2}}{\sum_{i=1}^{3}\left\lvert
N_{ei}\right\rvert^{2}}\simeq\dfrac{\Gamma\left(\tau\rightarrow\mu\nu_{\mu}\nu_{\tau}\right)^{\text{SM}}}{\Gamma\left(\tau\rightarrow
e\nu_{e}\nu_{\tau}\right)^{\text{SM}}}\left(R_{\mu e}^{\tau}\right)^{2},$ (24)

whereas the $\tau-\mu$ sector is constrained by:

$\dfrac{\Gamma\left(\tau\rightarrow e\nu\nu\right)}{\Gamma\left(\mu\rightarrow
e\nu\nu\right)}=\dfrac{\Gamma\left(\tau\rightarrow
e\nu_{e}\nu_{\tau}\right)^{\text{SM}}}{\Gamma\left(\mu\rightarrow
e\nu_{e}\nu_{\mu}\right)^{\text{SM}}}\dfrac{\sum_{i=1}^{3}\left\lvert N_{\tau
i}\right\rvert^{2}}{\sum_{i=1}^{3}\left\lvert N_{\mu
i}\right\rvert^{2}}\simeq\dfrac{\Gamma\left(\tau\rightarrow
e\nu_{e}\nu_{\tau}\right)^{\text{SM}}}{\Gamma\left(\mu\rightarrow
e\nu_{e}\nu_{\mu}\right)^{\text{SM}}}\left(R_{\tau\mu}^{\tau}\right)^{2}.$
(25)

Here $R^{\tau}_{\alpha\beta}$ has the same $\eta$-dependence as in Eq. (23).
The correlations between the ratios extracted from $\tau$ decays are taken
from [78]. Notice that $R^{K}_{\tau\mu}$ has not been included in Table 1 to
avoid double counting since these decays will be included individually as
independent measurements of $V_{us}$ as described below.

###  2.4 Constraints from CKM unitarity

The unitarity of the first row of the CKM is also constrained with significant
accuracy. While violations of the unitarity of the leptonic mixing matrix
leave the CKM matrix unchanged, the processes by which the values of its
elements are extracted involve weak decays and will inherit a dependence on
the $\eta$ parameters.

Since the CKM matrix is still unitary, the following CKM unitarity relation
holds

$\left\lvert V_{ud}\right\rvert^{2}+\left\lvert
V_{us}\right\rvert^{2}+\left\lvert V_{ub}\right\rvert^{2}=1\,,$ (26)

where, given the uncertainty on the value of $\left\lvert V_{us}\right\rvert$,
the element $\left\lvert V_{ub}\right\rvert=3.82(20)\cdot 10^{-3}$ can be
safely neglected and thus one can substitute:

$\left\lvert V_{ud}\right\rvert=\sqrt{1-\left\lvert
V_{us}\right\rvert^{2}}\,.$ (27)

We will treat $\left\lvert V_{us}\right\rvert$ as a nuisance parameter which
will be marginalised over.

The most precise determination of $\left\lvert V_{ud}\right\rvert$ comes from
superallowed $\beta$ decays. As such, its extraction will be modified by the
leptonic $W$ vertex by ($1-2\eta_{ee}$) and additionally by $G_{F}^{2}$ with
($1+2\eta_{ee}+2\eta_{\mu\mu}$), which in the end amounts to:

$\left\lvert V^{\beta}_{ud}\right\rvert=\sqrt{1-\left\lvert
V_{us}\right\rvert^{2}}\left(1+\eta_{\mu\mu}\right).$ (28)

On the other hand, $\left\lvert V_{us}\right\rvert$ can be determined from
semileptonic and leptonic $K$ decays, as well as $\tau$ decays. For the case
of the $K$ semileptonic decays, similarly to the superallowed $\beta$ decays,
the $\eta$ dependence will come from the $W$ vertex and the indirect
dependence of $G_{F}$. In particular:

$\displaystyle\Big{|}V^{K\rightarrow\pi e\nu\;}_{us}\Big{|}$
$\displaystyle=\left\lvert V_{us}\right\rvert\left(1+\eta_{\mu\mu}\right),$
(29) $\displaystyle\Big{|}V^{K\rightarrow\pi\mu\nu}_{us}\Big{|}$
$\displaystyle=\left\lvert V_{us}\right\rvert\left(1+\eta_{ee}\right).$ (30)

Moreover, one can also constrain the ratio $\left\lvert
V_{us}\right\rvert/\left\lvert V_{ud}\right\rvert$ by means of the ratio
between $K$ and $\pi$ leptonic decay widths:

$\left\lvert\dfrac{V_{us}}{V_{ud}}\right\rvert^{K,\pi\rightarrow\mu\nu}=\dfrac{\left\lvert
V_{us}\right\rvert}{\sqrt{1-\left\lvert V_{us}\right\rvert^{2}}}\,.$ (31)

In this case, there is no dependence on the $\eta$ parameters, as they cancel
out due to the fact that both final states have the same flavor. However, this
ratio is still useful to constrain the nuisance parameter $V_{us}$.

Regarding $\tau$ decays, $\left\lvert V_{us}\right\rvert$ can be extracted
from its decay to a $K$:

$\left\lvert V^{\tau\rightarrow K\nu}_{us}\right\rvert=\left\lvert
V_{us}\right\rvert\left(1+\eta_{ee}+\eta_{\mu\mu}-\eta_{\tau\tau}\right),$
(32)

where the $\eta$ dependence arises from the $W$ vertex with the $\tau$ and
from $G_{F}$. Additionally, another determination can be made from the ratio
of the $\tau$ decay widths to $K$ and $\pi$. This ratio, as in the case of Eq.
(31), does not depend on $\eta$. However, as one needs to multiply it by the
determination of $\left\lvert V_{ud}\right\rvert$ to get $\left\lvert
V_{us}\right\rvert$, it inherits the dependence of Eq. (28).

$\left\lvert V^{\tau\rightarrow K,\pi}_{us}\right\rvert=\left\lvert
V_{us}\right\rvert\left(1+\eta_{\mu\mu}\right).$ (33)

It should be noted that after a reassessment of radiative corrections to the
neutrino and superallowed $\beta$ decays, the estimated value of $V_{ud}$
decreased, leading to a 4-$5\sigma$ tension with the assumption of CKM
unitarity known as the Cabibbo anomaly [91, 92, 93, 94]. More recently, the
uncertainty associated to this measurement has also been revised and increased
due to a more conservative estimate of the nuclear structure uncertainties
[76]. This has reduced the Cabibbo anomaly to a 2-$3\sigma$ effect.
Nevertheless, for a positive-definite $\eta$-matrix, this anomaly is only
worsened for $\eta_{\mu\mu}>0$, as shown by Eq. (28) and therefore will push
the global fit to very small values of this parameter.

###  2.5 About the CDF-II $\boldsymbol{M_{W}}$ measurement

Given its accuracy, one of the most important constraints on the unitarity of
the leptonic mixing matrix and heavy neutrino mixing is the comparison of
$M_{W}$ with its value obtained from $G_{F}$ from muon decay, and therefore
affected by $\eta_{ee}+\eta_{\mu\mu}$ as discussed above, see Eq. (12).

However, the most recent (and most precise) determination of $M_{W}$ by the
CDF-II collaboration [75] is around 7$\sigma$ larger than the SM and around
$3\sigma$ with respect to the prior world average. Intriguingly, this larger
value of $M_{W}$, in tension with the SM, can in principle be explained
through non-unitarity parameters for positive definite
$\eta_{ee}+\eta_{\mu\mu}$, as shown in Table 1. This possibility was explored
in Ref. [95] considering a subset of the observables studied here including in
particular the invisible width of the $Z$, LFU ratios and tests of CKM
unitarity. The conclusion of that study was that the combined fit to the CDF-
II $M_{W}$ measurement plus the other observables was in significant tension,
mainly due to the Cabibbo anomaly which, at the time, was estimated to be
$4$-$5\sigma$ prior to the revision of the uncertainty in the determination of
$V_{ud}$. Indeed, the Cabibbo anomaly would prefer negative values for $\eta$,
and it can thus only be worsened when the CDF-II $M_{W}$ anomaly is
accommodated.

The recent revision of the significance of the Cabibbo anomaly to the
$2$-$3\sigma$ level invites to reconsider the explanation of the CDF-II
anomaly through non-unitarity. Nevertheless, we find that this measurement is
still in too much tension with, not only the other determinations of $M_{W}$,
but also with other very precise measurements such as $s_{\text{eff}}^{2}$ and
the Z-pole observables. In fact, even though the tension between the CDF-II
measurement of $M_{W}$ and its SM prediction could be explained by non-
unitarity, $M_{W}$, $s_{\text{eff}}^{2}$ and the Z-pole observables all depend
on the same combination of $\eta$ parameters (i.e. $\eta_{ee}+\eta_{\mu\mu}$).
Thus, any tension between these measurements cannot be improved by the
presence of non-unitarity.

We quantify this tension through the parameter goodness-of-fit (p-g.o.f.)
[96], which is particularly suited to explore the situation in which two or
more sets of observables are in tension. In practice it amounts to splitting
the dataset into two (or more) subsets, $A_{1}$ and $A_{2}$, and computing:

$\overline{\chi}^{2}=\chi^{2}_{12}-\chi^{2}_{1}-\chi^{2}_{2}\,,$ (34)

where $\chi^{2}_{i}$ is the minimum of the $\chi^{2}$ considering the dataset
$A_{i}$, and $\chi^{2}_{12}$ is the $\chi^{2}$ minimum considering both
datasets. In the following, we take $A_{1}$ as the CDF-II $M_{W}$ measurement
alone and consider three other options for the second set of data:

  * $\bullet$

$A_{2}$ as the dataset formed by $s^{2\hskip
4.26773pt\text{LHC}}_{\text{eff}}$, $s^{2\hskip
4.26773pt\text{Tev}}_{\text{eff}}$ and the world average $M_{W}$ without CDF-
II.

  * $\bullet$

$A_{2}$ as the Z-pole observables.

  * •

$A_{2}$ as the combination of the two previous sets of data.

The results of these three p-g.o.f. are summarised in Table 2. As expected,
adding the measurements of $s_{\text{eff}}^{2}$ increases the tension already
present between the new CDF-II result and other determinations of $M_{W}$,
pushing it above the 4$\sigma$ level. Moreover, the tension of CDF-II with the
Z-pole observables alone is also above $4\sigma$, indicating a high
incompatibility due to the fact that the Z-pole observables are in good
agreement with the SM expectation. Finally, the combination of the non-CDF-II
$M_{W}$ measurements, effective weak angle and Z-pole observables results in a
tension above $5\sigma$, clearly indicating that the CDF-II measurement cannot
be reconciled with the aforementioned observables through unitarity deviations
of the leptonic mixing matrix, since they all share the same dependence on
these parameters and, therefore, we choose not to include it among the
observables of our global fit.

Set of observables | p-g.o.f. | tension | p-value
---|---|---|---
CDF-II vs $M_{W}/s^{2}_{\text{eff}}$ | $21.75/1$ | 4.7$\sigma$ | $3.1\cdot 10^{-6}$
CDF-II vs Z-pole | $21.48/1$ | 4.6$\sigma$ | $3.6\cdot 10^{-6}$
CDF-II vs $M_{W}/s^{2}_{\text{eff}}$ and Z-pole | $27.30/1$ | 5.2$\sigma$ | $1.7\cdot 10^{-7}$
Table 2: Tension between different sets of observables, quantified through the
parameter goodness-of-fit [96]. Note that the CDF-II measurement of $M_{W}$ is
in tension not only with other determinations of $M_{W}$ and $s_{\text{eff}}$,
but also with the $Z$-pole observables measured at LEP.

###  2.6 Charged Lepton Flavor Violation

In the presence of heavy neutrinos, charged lepton flavor violating processes
are no longer protected by the GIM mechanism [97], as both non-unitarity of
the leptonic mixing matrix and the scale separation between light and heavy
neutrinos prevent such cancellation. Therefore, given the strong experimental
constraints on cLFV transitions, they can be used to derive bounds on the off-
diagonal elements of the $\eta$-matrix.

Currently, the most relevant cLFV processes for probing heavy neutrinos
include [28, 38] radiative decays, three body leptonic decays and $\mu-e$
conversion in heavy nuclei, whose present bounds are summarised in Table 3.
The available list of cLFV observables is actually much longer, including for
example decays of the $Z$ [31] and Higgs [27, 33, 43] bosons, although their
current sensitivities for non-unitarity effects are lower222Nevertheless, a
future Tera-Z factory would obtain competitive bounds from LFV $Z$ decays [46,
50]..

Observable | Experimental bound
---|---
$\mu\rightarrow e\gamma$ |  $4.2\cdot 10^{-13}$ [98]
$\tau\rightarrow e\gamma$ |  $3.3\cdot 10^{-8}$ [99]
$\tau\rightarrow\mu\gamma$ |  $4.2\cdot 10^{-8}$ [100]
$\mu\rightarrow eee$ |  $1.0\cdot 10^{-12}$ [101]
$\tau\rightarrow eee$ |  $2.7\cdot 10^{-8}$ [102]
$\tau\rightarrow\mu\mu\mu$ |  $2.1\cdot 10^{-8}$ [102]
$\mu\rightarrow e$ (Ti) |  $4.3\cdot 10^{-12}$ [103]
$\mu\rightarrow e$ (Au) |  $7.0\cdot 10^{-13}$ [104]
Table 3: Summary of present 90%CL upper limits for the branching ratios of the
most important cLFV observables constraining off-diagonal elements of the
matrix $\eta$.

Radiative $\ell_{\alpha}\to\ell_{\beta}\gamma$ decays are very well-studied
and have so far dominated the constraints set by cLFV processes. Their
complete rates induced by heavy neutrinos were first reported in Refs. [16,
17, 18, 19, 20, 21]. In the limit of heavy neutrinos, i.e. heavier than
$M_{W}$, and when their potential mass differences can be neglected333In
principle, non-degenerated heavy neutrinos could lead to numerical
cancellations in some of the cLFV processes [105]. Nevertheless, we consider
such a situation very unlikely, as we explain in App. B, and therefore do not
take into account in our global analysis., these rates can be approximated by

${\rm
BR}(\ell_{\alpha}\to\ell_{\beta}\gamma)\simeq\frac{3\alpha}{2\pi}\,\big{|}\eta_{\alpha\beta}\big{|}^{2}\,,$
(35)

and therefore they are directly related to the off-diagonal entries of the
non-unitarity matrix $\eta$. Notice moreover that this relation is independent
of the heavy neutrino mass, as long as it is heavy enough, implying that the
radiative decays impose mass-independent bounds on $\eta_{\alpha\beta}$.

On the other hand, the relation between $\eta$ and other cLFV transitions is
in general more involved. The reason is that, contrary to the radiative
decays, other cLFV observables also get contributions from neutral currents
from the $Z$ and $H$ boson penguins, as well as from box-diagrams, which have
a different dependence on heavy neutrino mixings at large masses [38, 51, 52,
53, 56], and thus a different $\eta$-dependence.

As an example to illustrate this behavior, let us consider a minimal scenario
with two almost degenerated neutrinos forming a single pseudo-Dirac pair of
mass $M$ (see section 3). Then, the rate for $\mu-e$ conversion in nuclei [38]
can be approximated by

$\displaystyle{\rm CR}(\mu\to e)\simeq$
$\displaystyle\frac{\alpha^{5}m_{\mu}^{5}G_{F}^{2}\,F_{p}^{2}}{10368\pi^{4}s_{\mathrm{w}}^{2}\Gamma_{\text{capt}}}\dfrac{Z_{\text{eff}}^{4}}{Z}\,\big{|}\eta_{e\mu}\big{|}^{2}\,\Big{|}\left(A+Z\right)F_{u}+\left(2A-Z\right)F_{d}\Big{|}^{2}\,,$
(36)

with

$\displaystyle F_{u}$
$\displaystyle=-27-148s_{\mathrm{w}}^{2}-(27-64s_{\mathrm{w}}^{2})\log\frac{M^{2}}{M_{W}^{2}}-(18-48s_{\mathrm{w}}^{2})\frac{M^{2}}{M_{W}^{2}}\text{Tr}\left[\eta\right],$
(37) $\displaystyle F_{d}$
$\displaystyle=-27+\phantom{1}74s_{\mathrm{w}}^{2}+(27-32s_{\mathrm{w}}^{2})\log\frac{M^{2}}{M_{W}^{2}}+(18-24s_{\mathrm{w}}^{2})\frac{M^{2}}{M_{W}^{2}}\text{Tr}\left[\eta\right],$
(38)

and $A,Z(Z_{\rm eff}),F_{p}$ and $\Gamma_{\rm capt}$ are properties of the
nucleus: mass number, (effective) atomic number, nuclear form factor and
capture rate, respectively. Notice that we chose $\mu-e$ conversion rate as an
example, but similar expressions can be derived for other cLFV processes such
as $\mu\to eee$ or $\tau\to eee$.

Besides the mild logarithmic dependence, there is a new $M$-dependent
contribution, not present in the current with photons and thus in
$\ell_{\alpha}\to\ell_{\beta}\gamma$. This new term is suppressed by an
additional power of the small $\eta$-matrix, encoded as
$\text{Tr}\left[\eta\right]$ in this simplified scenario, but at the same time
it is enhanced by the new mass scale, so it can still be important for
sufficiently heavy neutrinos. This also happens for other cLFV processes such
as the three-body leptonic decays, with additional $M$-dependent terms
involving also the diagonal entries of $\eta$. Consequently, processes such
$\mu\to eee$ and $\mu-e$ conversion in nuclei do not impose mass-independent
bounds on $\left\lvert\eta_{e\mu}\right\rvert$ in the $M\gg M_{W}$ limit.
Moreover, the experimental limits for these processes cannot be
straightforwardly translated to bounds for
$\left\lvert\eta_{e\mu}\right\rvert$, given the additional dependence on
$\text{Tr}\left[\eta\right]$. Thus, in all generality, the
$\left\lvert\eta_{e\mu}\right\rvert$ bounds depend not only on the mass but
also on the diagonal entries of $\eta$.

Figure 1:  Upper bounds on $\left\lvert\eta_{e\mu}\right\rvert$
($\left\lvert\eta_{e\tau}\right\rvert$) imposed by current limits of $\mu\to
e\gamma$ ($\tau\to e\gamma$), $\mu\to eee$ ($\tau\to eee$) and $\mu-e$
conversion in gold (see Table 3). We considered the minimal scenario described
in section 3 with a normal ordering for light neutrino masses. The color bands
reflect how the upper bound on $\left\lvert\eta_{\alpha\beta}\right\rvert$ is
modified as we vary the other free parameters. The upper right area in red
corresponds to Yukawa couplings larger than 1, outside the perturbativity
region.

In order to understand better and to quantify these mass and $\eta$
dependencies, we show in Fig. 1 the upper bounds that can be set on
$\left\lvert\eta_{e\mu}\right\rvert$ and $\left\lvert\eta_{e\tau}\right\rvert$
from each of the cLFV rates individually. This figure corresponds to the most
minimal scenario with only a degenerate pair of heavy neutrinos, which will be
introduced and explored in detail in section 3, but it is sufficient to
exemplify and discuss the overall behavior. Besides varying the heavy neutrino
mass, we also survey the values of the additional free parameters of the model
(see Section 3) in order to also scan $\text{Tr}\left[\eta\right]$ from its
minimum444Notice that in general $\text{Tr}\left[\eta\right]$ cannot be zero
for $\eta_{\alpha\beta}>0$ due to the Schwarz inequality in Eq. (7). For
example, in the minimal model of section 3, we have
$\text{Tr}\left[\eta\right]\geq 2\left\lvert\eta_{e\mu}\right\rvert$. to its
maximum value allowed in this scenario. This defines a color band for each
observable, showing how much the bound on $\eta_{\alpha\beta}$ changes as we
vary $\text{Tr}\left[\eta\right]$. Since the radiative decays do not depend on
$\text{Tr}\left[\eta\right]$, see Eq. (35), the band is only a line and it can
be seen that they quickly saturate to a constant value as soon as the heavy
neutrino mass becomes a few times heavier than $M_{W}$.

On the other hand, the rest of observables display a more complex behaviour.
Their bounds lie in a band, thicker for heavier masses since the importance of
the $\text{Tr}\left[\eta\right]$ terms grows with $M$, and they become more
stringent for heavier neutrino masses555With the exception of a cancellation
in the $\mu-e$ conversion rate in gold for masses around 20 TeV. The exact
position of this cancellation depends however on the mass and atomic number of
the nucleus, so it will be different for each nuclei., even overcoming those
from the radiative decays. Nevertheless, whether this crossing between
observables happens within the perturbative region of the model and where
exactly depend on the values of $\text{Tr}\left[\eta\right]$. In this
particular example, $\mu-e$ conversion dominates over $\mu\to e\gamma$ for
masses above 10 TeV approximately, although the exact position varies with
$\text{Tr}\left[\eta\right]$. The situation is similar in the $\tau$-$e$
sector, where $\tau\to eee$ can dominate over $\tau\to e\gamma$ for heavy
masses above the TeV range, again depending on the value of
$\text{Tr}\left[\eta\right]$. Notice that in the $\tau$-$e$ sector this
crossing happens for lower masses than in the $\mu$-$e$ sector, since in the
former we are probing larger values of $\eta$ and, thus, the additional
$\eta$-terms are less suppressed. Nevertheless, in both sectors these effects
become relevant only close to the perturbativity limit, since they require
large masses and mixings simultaneously.

For these reasons, which cLFV observable provides the most stringent bound is
model-dependent and also changes in different regions of the parameter space.
To deal with this fact and to provide the most robust constraints possible, we
will “marginalize over” the dependencies shown in Fig. 1 by selecting a bound
that would apply to all values of the mass and $\text{Tr}\left[\eta\right]$.
In practice we simply consider the constraint stemming from radiative decays
as given in Eq. (35). This choice slightly overestimates the bounds from
radiative decays for masses close to $M_{W}$, where the GIM cancellation
starts to be recovered. Nevertheless, in the $\mu$-$e$ sector this effect is
compensated by $\mu-e$ conversion in gold, which provides a strikingly similar
constraint precisely in that region of the parameter space. While this is not
the case in the $\tau$ sector, those cLFV bounds are subdominant compared to
those derived from the Schwarz inequality, so that its effect is not relevant
in the fit. Consequently, Eq. (35) turns out to be a good approach to a
conservative and mass-independent cLFV limit, and will be added to our fit. In
any event, it should be noted that in some regions of the parameter space, for
the heaviest masses, stronger constraints than our conservative estimate might
apply as shown Fig. 1.

##  3 Global fit bounds for the 2 neutrino case (2N-SS)

In order to reproduce the two distinct mass splittings that characterize the
neutrino oscillation phenomenon, at least two of the mainly-SM neutrinos need
to become massive. Therefore the minimum number of heavy neutrinos needed in
order to have a realistic neutrino mass model is two, as we consider in this
section. We will dub this setup with two heavy neutrinos as the 2N-SS, with RH
neutrinos $N_{R}$ and $N^{\prime}_{R}$ inducing the observed light neutrino
masses and mixings.

Moreover, in order for these RH neutrinos to have sizeable mixing with the
active neutrinos while generating radiatively stable and small neutrino
masses, a lepton number protected seesaw realization is required with the
heavy neutrinos forming a pseudo-Dirac pair [64, 65, 72, 66]. The most general
neutrino mass matrix which satisfies these characteristics, in the basis
$\left(\nu_{L}\hskip 4.26773ptN^{c}_{R}\hskip
4.26773pt{N^{\prime}}_{R}^{c}\right)^{T}$, has the following form:

$\mathcal{M}_{\nu}=\begin{pmatrix}0&Yv/\sqrt{2}&\epsilon
Y^{\prime}v/\sqrt{2}\\\ Y^{T}v/\sqrt{2}&\mu^{\prime}&M\\\
\epsilon{Y^{\prime}}^{T}v/\sqrt{2}&M&\mu\end{pmatrix},$ (39)

where the $\epsilon,\mu$ and $\mu^{\prime}$ terms softly break a generalized
lepton number symmetry, $L_{e}=L_{\mu}=L_{\tau}=L_{N}=-L_{N^{\prime}}=1$.
Therefore, it is technically natural to consider
$\epsilon,\mu/M,\mu^{\prime}/M\ll 1$. Indeed, in the limit when these terms
vanish, $\mathcal{M}_{\nu}$ preserves the symmetry, yielding three massless
neutrinos ($m_{\nu}=0$) and two degenerate heavy neutrinos forming a Dirac
fermion. However, even in this limit the active-heavy mixings are already non-
zero and hence unsuppressed by the smallness of neutrino masses

$\Theta_{\alpha}=\begin{pmatrix}0&\dfrac{Y_{\alpha}v}{\sqrt{2}M}\end{pmatrix}\equiv\begin{pmatrix}0&\theta_{\alpha}\end{pmatrix},$
(40)

leading to potentially sizable non-unitarity effects given by

$\eta=\dfrac{1}{2}\begin{pmatrix}\left\lvert\theta_{e}\right\rvert^{2}&\theta_{e}\theta_{\mu}^{*}&\theta_{e}\theta_{\tau}^{*}\\\
\theta_{e}^{*}\theta_{\mu}&\left\lvert\theta_{\mu}\right\rvert^{2}&\theta_{\mu}\theta_{\tau}^{*}\\\
\theta_{e}^{*}\theta_{\tau}&\theta_{\mu}^{*}\theta_{\tau}&\left\lvert\theta_{\tau}\right\rvert^{2}\end{pmatrix}.$
(41)

The minimality of this model has two important implications for our analysis.
On the one hand, this $\eta$-matrix saturates the Schwarz inequality (i.e.
$\left\lvert\eta_{\alpha\beta}\right\rvert=\sqrt{\eta_{\alpha\alpha}\eta_{\beta\beta}}$),
meaning that the cLFV bounds will constrain not only the off-diagonal entries
of $\eta$, but also the diagonal ones. Therefore the cLFV observables need to
be added into the global fit together with the LFC ones in Table 1. On the
other hand, the flavor structure of the mixing is subject to important
constraints from the requirement of explaining the correct light neutrino mass
matrix. In fact, this flavor structure can be reconstructed from the light
neutrino masses and mixings [72], except for an overall scale $\theta$:

$\displaystyle\theta_{\alpha}$
$\displaystyle=\dfrac{\theta}{\sqrt{2}}\left(\sqrt{1+\rho}\hskip
2.84544ptU^{*}_{\alpha 3}+\sqrt{1-\rho}\hskip 2.84544ptU^{*}_{\alpha
2}\right)$ for Normal Ordering (NO), (42) $\displaystyle\theta_{\alpha}$
$\displaystyle=\dfrac{\theta}{\sqrt{2}}\left(\sqrt{1+\rho}\hskip
2.84544ptU^{*}_{\alpha 2}+\sqrt{1-\rho}\hskip 2.84544ptU^{*}_{\alpha
1}\right)$ for Inverted Ordering (IO), (43)

where

$\displaystyle\rho$ $\displaystyle=\dfrac{\sqrt{\Delta
m^{2}_{31}}-\sqrt{\Delta m^{2}_{21}}}{\sqrt{\Delta m^{2}_{31}}+\sqrt{\Delta
m^{2}_{21}}}$ for NO, (44) $\displaystyle\rho$
$\displaystyle=\dfrac{\sqrt{\Delta m^{2}_{23}}-\sqrt{\Delta m^{2}_{23}-\Delta
m^{2}_{21}}}{\sqrt{\Delta m^{2}_{23}}+\sqrt{\Delta m^{2}_{23}-\Delta
m^{2}_{21}}}$ for IO, (45)

with $\Delta m^{2}_{ij}=m^{2}_{i}-m^{2}_{j}$ and where we use for $U$ the
standard PDG parametrization with a single Majorana phase, as one neutrino
remains massless.

From oscillation experiments the values of the two light neutrino mass
splittings and the three mixing angles have been determined with good
accuracy. However, the Dirac CP phase and Majorana phases are still mostly
unconstrained. As such, in our analysis we will fix the mass splittings and
the mixing angles to their best-fit values as reported by NuFIT [106], while
treating the Dirac phase $\delta$, the Majorana phase $\phi$ and the overall
magnitude666Actually $\theta$ is a complex parameter. However, only its
modulus is relevant for our analysis, as its phase only appears in lepton
number violating processes, which are not studied in this work. of the mixing
$\theta$ as free parameters for our fit. All in all, this scenario is fully
described by 3 free parameters.

Figure 2:  Profiles from our global fit analysis of the minimal model with 2
RH neutrinos (2N-SS), for NO (upper panels) and IO (lower panels). In each
panel, we minimized over all the parameters of the fit but the corresponding
$\eta_{\alpha\alpha}$. The dashed green lines are obtained following Wilks’
theorem, while the solid black lines are the result of calibrating the test
statistic by bootstrapping it (see Appendix A).

The results of our statistical analysis for this 2N-SS setup are shown in Fig.
2, where we plot the profiles for the diagonal entries of the $\eta$-matrix
for the two neutrino mass orderings. In each panel, the dashed green line is
obtained following Wilks’ theorem, while the solid black line is the result of
calibrating our test statistic through the bootstrapping procedure (see App.
A), since, given the strong correlations among the observables implied by Eqs.
(42) and (43), the requirements for Wilks’ theorem to apply are not met.
Nevertheless, we find that in this case Wilks’ theorem still provides a good
approximation. Indeed, given the little freedom available to this rather
constrained scenario, the bound on the non-unitarity parameter $\theta$ mainly
stems from the best constrained observable in the fit, namely $\mu\to
e\gamma$, and, neglecting the contributions of the other observables, the
requirements for Wilks’ theorem to apply are approximately met. However, we
will see that this is not the case for other, less constrained, scenarios.

2N-SS | Normal Ordering | Inverted Ordering
---|---|---
|  $68\%$CL |  $95\%$CL |  $68\%$CL |  $95\%$CL
$\eta_{ee}=\dfrac{\left\lvert\theta_{e}\right\rvert^{2}}{2}$ | $6.4\cdot 10^{-6}$ | $9.4\cdot 10^{-6}$ | $[0.98,4.4]\cdot 10^{-4}$ | $5.5\cdot 10^{-4}$
$\eta_{\mu\mu}=\dfrac{\left\lvert\theta_{\mu}\right\rvert^{2}}{2}$ | $6.9\cdot 10^{-5}$ | $1.3\cdot 10^{-4}$ | $[0.20,1.0]\cdot 10^{-6}$ | $3.2\cdot 10^{-5}$
$\eta_{\tau\tau}=\dfrac{\left\lvert\theta_{\tau}\right\rvert^{2}}{2}$ | $8.6\cdot 10^{-5}$ | $2.1\cdot 10^{-4}$ | $[0.94,2.8]\cdot 10^{-5}$ | $4.5\cdot 10^{-5}$
$\text{Tr}\left[\eta\right]=\dfrac{\left\lvert\theta\right\rvert^{2}}{2}$ | $1.6\cdot 10^{-4}$ | $2.9\cdot 10^{-4}$ | $[1.1,4.8]\cdot 10^{-4}$ | $6.0\cdot 10^{-4}$
$\left\lvert\eta_{e\mu}\right\rvert=\dfrac{\left\lvert\theta_{e}\theta_{\mu}^{*}\right\rvert}{2}$ | $8.3\cdot 10^{-6}$ | $1.2\cdot 10^{-5}$ | $[0.37,1.0]\cdot 10^{-5}$ | $1.3\cdot 10^{-5}$
$\left\lvert\eta_{e\tau}\right\rvert=\dfrac{\left\lvert\theta_{e}\theta_{\tau}^{*}\right\rvert}{2}$ | $1.5\cdot 10^{-5}$ | $2.2\cdot 10^{-5}$ | $[0.25,1.2]\cdot 10^{-4}$ | $1.4\cdot 10^{-4}$
$\left\lvert\eta_{\mu\tau}\right\rvert=\dfrac{\left\lvert\theta_{\mu}\theta_{\tau}^{*}\right\rvert}{2}$ | $7.2\cdot 10^{-5}$ | $1.3\cdot 10^{-4}$ | $[0.38,3.0]\cdot 10^{-6}$ | $3.5\cdot 10^{-5}$
Table 4: Upper bounds (or preferred intervals) for the most minimal set-up
with two heavy neutrinos forming a single pseudo-Dirac pair (2N-SS), which are
obtained from the boostrapped profiles in Fig. 2 and the equivalent ones (not-
shown) for $\text{Tr}\left[\eta\right]$ and the off-diagonal elements. Note
that these results directly apply to $\eta$ and to (half of) the squared
active-sterile mixings $|\theta_{\alpha}|^{2}$. They can also be easily
translated to the $\alpha$-parametrization, as detailed in the text.

The resulting 68% and 95%CL upper limits (or preferred intervals in some
cases) are summarised in Table 4. Here and in all the RH neutrino scenarios
under study, we also provide the constraints derived for
$\text{Tr}\left[\eta\right]$, which in this simplest case amounts to
$\theta^{2}/2$, as a measurement of how large the _total_ deviation from
unitarity of the whole matrix is allowed to be regardless of its particular
flavour structure. Indeed, $\text{Tr}\left[\eta\right]$ is an invariant under
changes of basis and, given that $\eta$ is positive definite by construction,
its trace corresponds to the sum of its three eigenvalues.

It is interesting to note that in this scenario the results are rather
different between the profiles for normal and inverted orderings. Besides the
different ranges for each $\eta_{\alpha\alpha}$, as it can be seen in Fig. 2,
we find a non-unitary best-fit point for IO, while this does not happen for
NO. This qualitative difference, which also conditions the different values
obtained for the allowed ranges, can be understood as an interplay between the
preference of the data and the constrained flavor structure of $\eta$ in this
model.

Figure 3: Regions in the RH neutrino mixing flavor space consistent with
current neutrino oscillation data, fixing the mixing angles to their best fit
values, and varying the phases and absolute neutrino mass. The darker regions
correspond to the minimal scenario with 2 RH neutrinos (2N-SS) [107] described
in section 3, while the lighter regions are for the next-to-minimal scenario
with 3 RH neutrinos (3N-SS) [54] studied in section 4, where we include the
cosmological upper bound on neutrino masses. The black star corresponds to the
best-fit point for the general neutrino case (G-SS, see section 5), where the
hole triangle is allowed, indicating the actual preference of the data set
given in Table 1. Note that in the 2N-SS and 3N-SS, the axes correspond to the
normalized squared mixings since
$\eta_{\alpha\alpha}/\text{Tr}\left[\eta\right]=|\theta_{\alpha}|^{2}/|\theta|^{2}$.

On the one hand, in Table 1 a mild ($\sim 1\sigma$) preference for a non-zero
$\eta_{ee}+\eta_{\mu\mu}$ from $M_{W}$ and $s_{\mathrm{w}}^{2}$ can be seen,
as well as an also mild ($\sim 1\sigma$) preference for
$\eta_{ee}>\eta_{\mu\mu}$ from the LFU ratios. Furthermore, the CKM data
strongly disfavors non-zero values of $\eta_{\mu\mu}$, as they can only worsen
the Cabibbo anomaly, and there is also no preference for a non-zero
$\eta_{\tau\tau}$. On top of that, the strong cLFV bounds777 Notice that in
the 2N-SS we expect mixings to all three flavors, since all of them are
proportional to the overall scale $\theta$. This implies that the cLFV rates
cannot be avoided by turning just a single $\theta_{\alpha}$ off, and thus
they are specially constraining for this 2N-SS scenario. on the $\mu$-$e$
sector require either the electron or the muon mixing to be very small.
Combining all these aspects, we obtain that the data prefers a non-zero
$\eta_{ee}$ with suppressed $\eta_{\mu\mu}$ and $\eta_{\tau\tau}$. This is
represented with a black star in Fig. 3, where we display the possible flavor
patterns for both normal and inverted orderings and for the different
scenarios under consideration.

On the other hand, the flavor structure of the 2N-SS is very restricted. It is
determined up to the mass ordering and the unknown phases $\delta$ and $\phi$,
and the resulting regions in flavor space are shown as darker areas in Fig. 3.
In particular, the NO case is characterised by having a suppressed
$\left\lvert\theta_{e}\right\rvert$ with respect to
$\left\lvert\theta_{\mu}\right\rvert$ and
$\left\lvert\theta_{\tau}\right\rvert$ (for all values of $\delta$ and
$\phi$), which is precisely the opposite of what the data prefers. This is
manifest in the figure, where the dark red region is far away from the data-
preferred black star. Therefore, the NO has the best fit-point at $\eta=0$,
with stronger bounds for $\eta_{ee}$ than for the other flavors. On the
contrary, the flavor structure of the IO case is such that, for certain values
of the phases, the mixing to the electron overcomes that of the muon and tau,
and thus it can accommodate better the preference of the data. We see it again
from Fig. 3, as the dark blue region approaches more to the black star.
Consequently, the IO has a non-trivial best-fit point, with a mild ($\sim
1\sigma$) preference for non-zero $\eta$, where the best-fit for
$\eta_{\mu\mu}$ is suppressed with respect to $\eta_{ee}$ and
$\eta_{\tau\tau}$. Consistently, we find a $\Delta\chi^{2}=2.07$ in favour of
the IO minimum with respect to the NO one.

Regarding the off-diagonal elements, the bound on
$\left\lvert\eta_{e\mu}\right\rvert$ is completely dominated by the
corresponding cLFV bound. Conversely, the bounds on
$\left\lvert\eta_{e\tau}\right\rvert$ and
$\left\lvert\eta_{\mu\tau}\right\rvert$ are much stronger than the bounds
derived from their corresponding cLFV processes. This effect is related to the
saturation of the Schwarz inequality in Eq. (7) and the strong correlations
present in the flavor structure of the 2N-SS: since all three mixings are
proportional to a common scale, this also implies
$\left\lvert\eta_{e\mu}\right\rvert\propto\left\lvert\eta_{e\tau}\right\rvert\propto\left\lvert\eta_{\mu\tau}\right\rvert$
where the proportionality depends on the phases $\delta$ and $\phi$. Thus, the
very strong bound on $\left\lvert\eta_{e\mu}\right\rvert$ also induces quite
stringent bounds on $\left\lvert\eta_{e\tau}\right\rvert$ and
$\left\lvert\eta_{\mu\tau}\right\rvert$.

It should be noted that, even though the Schwarz inequality is saturated and
one can reconstruct the off-diagonal elements from the diagonal ones, the
bounds for the former cannot be inferred from the latter. This is a
consequence of the strong correlations imposed by the cLFV constraints and of
the fact that, in order to obtain each bound, the rest of the parameters are
profiled over. Therefore, it is not possible to saturate simultaneously all
bounds and it should be checked that all constraints are satisfied for a given
mixing pattern.

As an alternative to $\eta$, a lower triangular parametrization [108, 109] has
been shown to be more convenient for the study of non-unitarity in the
neutrino oscillation phenomenon [110]:

$\alpha=\begin{pmatrix}\alpha_{11}&0&0\\\ \alpha_{21}&\alpha_{22}&0\\\
\alpha_{31}&\alpha_{32}&\alpha_{33}\end{pmatrix}.$ (46)

These $\alpha$ parameters can be straightforwardly mapped to the $\eta$-matrix
[110], thus our bounds can be easily translated to this parametrization. In
particular, at $95\%$CL we obtain:

NO: $\displaystyle\hskip
14.22636pt\left\lvert\mathbb{I}-\alpha\right\rvert<\begin{pmatrix}9.4\cdot
10^{-6}&0&0\\\ 2.4\cdot 10^{-5}&1.3\cdot 10^{-4}&0\\\ 4.4\cdot
10^{-5}&2.6\cdot 10^{-4}&2.1\cdot 10^{-4}\end{pmatrix}\,,$ (47) IO:
$\displaystyle\hskip
14.22636pt\left\lvert\mathbb{I}-\alpha\right\rvert<\begin{pmatrix}5.5\cdot
10^{-4}&0&0\\\ 2.6\cdot 10^{-5}&3.2\cdot 10^{-5}&0\\\ 2.8\cdot
10^{-4}&7.0\cdot 10^{-5}&4.5\cdot 10^{-5}\end{pmatrix}.$ (48)

##  4 Global fit bounds for the 3 neutrino case (3N-SS)

The scenario in the previous section with two heavy neutrinos is the most
minimal set-up to accommodate oscillation data, however it predicts the
lightest neutrino to be massless. While this is in agreement with current
observations, the overall scale of neutrino masses still remains unknown and,
thus, it is perfectly possible that all light neutrinos are massive. In that
case, the minimum number of extra neutrinos needed to accommodate light
neutrino masses is three, in line with the three SM generations for all other
fermions. In this section we will investigate the bounds that may be derived
in this scenario that we dub 3N-SS.

Similarly to the previous section, the only way in which sizeable mixing can
be obtained while keeping neutrino masses small and stable under radiative
corrections is via a lepton number protected Lagrangian. More precisely, in
the basis $\left(\nu_{L}\hskip 4.26773ptN^{c}_{R}\hskip
4.26773pt{N^{\prime}}_{R}^{c}\hskip
4.26773pt{N^{\prime\prime}}_{R}^{c}\right)^{T}$, the mass matrix reads:

$\mathcal{M}_{\nu}=\begin{pmatrix}0&Yv/\sqrt{2}&\epsilon_{1}Y^{\prime}v/\sqrt{2}&\epsilon_{2}Y^{\prime\prime}v/\sqrt{2}\\\
Y^{T}v/\sqrt{2}&\mu_{1}&M&\mu_{3}\\\
\epsilon_{1}{Y^{\prime}}^{T}v/\sqrt{2}&M&\mu_{2}&\mu_{4}\\\
\epsilon_{2}{Y^{\prime\prime}}^{T}v/\sqrt{2}&\mu_{3}&\mu_{4}&M^{\prime}\end{pmatrix},$
(49)

where $\epsilon_{1},\epsilon_{2}\ll 1$ and the $\mu_{i}$ parameters are small
compared with $M,M^{\prime}$. Assigning
$L_{e}=L_{\mu}=L_{\tau}=L_{N}=-L_{N^{\prime}}=1$ and $L_{N^{\prime\prime}}=0$,
in the lepton number conserving limit two heavy neutrinos arrange into a Dirac
fermion with non-vanishing mixing with the active neutrinos, while the third
heavy neutrino behaves as a decoupled Majorana state. Therefore, the structure
of the $\eta$ matrix is exactly the same as in Eq. (41) and thus the Schwarz
inequality is again saturated, which requires the inclusion of cLFV bounds in
combination with the LFC observables in our global fit.

However, the relations between the mixings to the active flavors are different
from the ones present in the 2N-SS, as the increased freedom in the parameter
space implies different (looser) correlations imposed by correctly reproducing
the observed neutrino mass matrix. In particular, it has been shown in [73]
that the mixing to one flavor is determined by the other two and the entries
of the light neutrino mass matrix. For instance, if one chooses to reconstruct
the $\tau$ mixing, then:

$\displaystyle\theta_{\tau}\simeq\frac{1}{m_{e\mu}^{2}-m_{ee}m_{\mu\mu}}$
$\displaystyle\Biggl{(}\theta_{e}\left(m_{e\mu}m_{\mu\tau}-m_{e\tau}m_{\mu\mu}\right)+\theta_{\mu}(m_{e\mu}m_{e\tau}-m_{\mu\tau}m_{ee})$
$\displaystyle+$
$\displaystyle\sqrt{\theta_{e}^{2}m_{\mu\mu}-2\theta_{e}\theta_{\mu}m_{e\mu}+\theta_{\mu}^{2}m_{ee}}\,\times$
$\displaystyle\times$
$\displaystyle\sqrt{m_{e\tau}^{2}m_{\mu\mu}-2m_{e\mu}m_{e\tau}m_{\mu\tau}+m_{ee}m_{\mu\tau}^{2}+m_{e\mu}^{2}m_{\tau\tau}-m_{ee}m_{\mu\mu}m_{\tau\tau}}\Biggr{)},$
(50)

where $m_{\alpha\beta}\equiv\left(m_{\nu}\right)_{\alpha\beta}=\left(Um_{\rm
diag}\,U^{T}\right)_{\alpha\beta}$ are the entries of the light neutrino mass
matrix.

In our analysis, we fix the mass splittings and the mixing angles to their
best-fit values, leaving the Dirac phase $\delta$, the two Majorana phases
$\phi_{1}$ and $\phi_{2}$, and the lightest neutrino mass
$m_{\text{lightest}}$ as the only free parameters in $m_{\nu}$. Additionally,
we have also as free parameters the mixings $\theta_{e}$ and $\theta_{\mu}$,
both in modulus and phase. In total, this setup is characterized by 8 free
parameters. Among these parameters, the only one that is constrained by
experimental data is $m_{\text{lightest}}$, which is currently bounded from
kinematical searches at KATRIN [111] and, more strongly, from cosmology [112].
We provide our results considering the more stringent bound from Planck of
$\sum m_{\nu}<0.12$ eV (95$\%$CL). Nevertheless, we have also performed the
analysis using the looser constraint from KATRIN and obtained very similar
results.

Figure 4:  Same as Fig. 2 but for the scenario with 3 RH neutrinos (3N-SS). In
this case, however, the calibration of the test statistic is done following
the profiled bootstrapping procedure described in Appendix A. We have used the
Planck upper bound [112] on neutrino masses, but similar results are obtained
using the KATRIN bound [111]. Notice the different (log)scale for
$\eta_{\mu\mu}$.

The results of the analysis are shown in Fig. 4 and Table 5. Given the large
dimensionality of the parameter space of this setup, the bootstrapping
procedure we followed in the 2N-SS is no longer tractable, as the number of
points one has to calibrate grows exponentially with the number of free
parameters. As a consequence, we are forced to perform an approximated
bootstrap procedure, which we will dub Profiled Bootstrap, and whose details
can be found in Appendix A.

3N-SS | Normal Ordering | Inverted Ordering
---|---|---
|  $68\%$CL |  $95\%$CL |  $68\%$CL |  $95\%$CL
$\eta_{ee}=\dfrac{\left\lvert\theta_{e}\right\rvert^{2}}{2}$ | $[0.28,0.99]\cdot 10^{-3}$ | $1.3\cdot 10^{-3}$ | $[0.31,1.0]\cdot 10^{-3}$ | $1.4\cdot 10^{-3}$
$\eta_{\mu\mu}=\dfrac{\left\lvert\theta_{\mu}\right\rvert^{2}}{2}$ | $1.3\cdot 10^{-7}$ | $1.1\cdot 10^{-5}$ | $1.2\cdot 10^{-7}$ | $1.0\cdot 10^{-5}$
$\eta_{\tau\tau}=\dfrac{\left\lvert\theta_{\tau}\right\rvert^{2}}{2}$ | $[0.3,3.9]\cdot 10^{-4}$ | $1.0\cdot 10^{-3}$ | $1.7\cdot 10^{-4}$ | $8.1\cdot 10^{-4}$
$\text{Tr}\left[\eta\right]=\dfrac{\left\lvert\theta\right\rvert^{2}}{2}$ | $[0.35,1.3]\cdot 10^{-3}$ | $1.9\cdot 10^{-3}$ | $[0.33,1.0]\cdot 10^{-3}$ | $1.5\cdot 10^{-3}$
$\left\lvert\eta_{e\mu}\right\rvert=\dfrac{\left\lvert\theta_{e}\theta_{\mu}^{*}\right\rvert}{2}$ | $8.5\cdot 10^{-6}$ | $1.2\cdot 10^{-5}$ | $8.5\cdot 10^{-6}$ | $1.2\cdot 10^{-5}$
$\left\lvert\eta_{e\tau}\right\rvert=\dfrac{\left\lvert\theta_{e}\theta_{\tau}^{*}\right\rvert}{2}$ | $[1.3,5.1]\cdot 10^{-4}$ | $9.0\cdot 10^{-4}$ | $3.3\cdot 10^{-4}$ | $8.0\cdot 10^{-4}$
$\left\lvert\eta_{\mu\tau}\right\rvert=\dfrac{\left\lvert\theta_{\mu}\theta_{\tau}^{*}\right\rvert}{2}$ | $5.0\cdot 10^{-6}$ | $5.7\cdot 10^{-5}$ | $3.8\cdot 10^{-6}$ | $1.8\cdot 10^{-5}$
Table 5: Upper bounds (or preferred intervals) for the 3N-SS scenario,
considering the cosmological upper bound on neutrino masses. The results for
the diagonal entries of $\eta$ are obtained from the profiled bootstrap in
Fig. 4, while the off-diagonal ones follow Wilks’ theorem, as explained in
App. A. Similar results are obtained when considering instead the KATRIN upper
bound.

We find that assuming Wilks’ theorem does not deviate much from the results of
the approximate bootstrapping, with the main differences usually appearing
close to the physical border at $\eta_{\alpha\alpha}=0$. For the $\eta_{ee}$
profile, we find a non-zero best-fit value and a slight enhancement of the CL
with respect to Wilks’ near $\eta_{ee}=0$ due to boundary effects expected
[113]. For the $\eta_{\mu\mu}$ profile, we find a more complex behaviour. For
lower CL, since the data prefers $\eta_{ee}\neq 0$, the $\eta_{\mu\mu}$
profile is dominated by the cLFV bound, which leads to a good agreement with
Wilks’ theorem and to a very stringent $68\%$ bound. However, at higher CL,
when $\eta_{ee}$ is allowed to vanish by the rest of the observables,
$\mu\rightarrow e\gamma$ ceases to be a relevant bound, and $\eta_{\mu\mu}$ is
instead constrained by the LFC observables, which are looser. As a result, the
profile develops a plateau-like feature at the same CL at which $\eta_{ee}=0$.
Furthermore, we observe an enhancement of the CL in the plateau region, with
respect to Wilks’ expectation coincident with the same enhancement for
$\eta_{ee}=0$ which causes the plateau.

Contrary to the 2N-SS case, the results do not depend significantly on the
neutrino mass ordering. In particular, we find a non-zero best-fit now for
both mass orderings, since the additional freedom with respect to the 2N-SS
setup allows for accommodating oscillation data with larger $\eta_{ee}$ also
for NO. This can be seen in Fig. 3, where the lighter regions cover a much
bigger area of the triangle for both NO and IO. The only substantial
difference is that, for NO, there is a slight preference ($1-2\sigma$) for a
non-zero $\eta_{\tau\tau}$. This is a result of the correlation between the
mixings present in the 3N-SS. As the best-fit point stands at $\eta_{ee}\neq
0$ and $\eta_{\mu\mu}=0$, given the correlation in Eq. (50) this also induces
a best-fit point for $\eta_{\tau\tau}\neq 0$ in the NO case. However, this
does not happen for IO, since its flavor structure allows for a suppressed
$\eta_{\tau\tau}$ even when $\eta_{ee}\neq 0$ (see Fig. 3). Consequently,
since the data prefers a suppressed $\eta_{\tau\tau}$, the IO case yields
again a slightly better fit with respect to the NO, although less pronounced
compared to the 2N-SS scenario. More precisely, we find a $\Delta\chi^{2}=0.5$
in favor of the IO best-fit with respect to the NO one.

Regarding the bounds on the off-diagonal entries, it is interesting to compare
the global bounds with those imposed directly by cLFV observables888These
bounds correspond to the LFV bound column in Table 6. following Eq. (35). We
see that the constraints on $\left\lvert\eta_{e\mu}\right\rvert$ are very
similar to the ones derived directly from cLFV processes, which reflects that
in the $\mu$-$e$ sector the test statistics is dominated by the strong bounds
on $\mu\to e$ transitions. In the $e$-$\tau$ sector, we find a mild
($1-2\sigma$) preference for a non-zero $\left\lvert\eta_{e\tau}\right\rvert$,
which is induced from the correlation given in Eq. (50) that implies a non-
zero $\left\lvert\theta_{\tau}\right\rvert$ for a non-zero
$\left\lvert\theta_{e}\right\rvert$. This effect is softer for IO, since, as
explained previously, its flavor structure allows for a suppressed
$\left\lvert\theta_{\tau}\right\rvert$ even when
$\left\lvert\theta_{e}\right\rvert$ is non-zero. Nevertheless, when focusing
on the $95\%$ C.L. upper limits, we see that the bounds in both the $e$-$\tau$
and $\mu$-$\tau$ sectors are much stronger than those derived from the
corresponding cLFV processes. This is especially true for
$\left\lvert\eta_{\mu\tau}\right\rvert$, which can be understood in terms of
the $\eta_{\mu\mu}$ profiles in Fig. 4: as previously stated, $\eta_{\mu\mu}$
is constrained to be very small up to relatively high CL due to the preference
for a non-zero $\eta_{ee}$ and the very strong cLFV bound on
$\left\lvert\eta_{e\mu}\right\rvert=\sqrt{\eta_{ee}\eta_{\mu\mu}}$. This tight
bound on $\eta_{\mu\mu}$ also induces a strong bound on
$\left\lvert\eta_{\mu\tau}\right\rvert=\sqrt{\eta_{\mu\mu}\eta_{\tau\tau}}$.
Note that this interplay is a consequence of the saturation of the Schwarz
inequality.

Finally, the bounds shown in Table 5 can be again translated to the lower-
triangular $\alpha$ parametrization obtaining, at $95\%$CL:

NO: $\displaystyle\hskip
14.22636pt\left\lvert\mathbb{I}-\alpha\right\rvert<\begin{pmatrix}1.3\cdot
10^{-3}&0&0\\\ 2.4\cdot 10^{-5}&1.1\cdot 10^{-5}&0\\\ 1.8\cdot
10^{-3}&1.1\cdot 10^{-4}&1.0\cdot 10^{-3}\end{pmatrix}\,,$ (51) IO:
$\displaystyle\hskip
14.22636pt\left\lvert\mathbb{I}-\alpha\right\rvert<\begin{pmatrix}1.4\cdot
10^{-3}&0&0\\\ 2.4\cdot 10^{-5}&1.0\cdot 10^{-5}&0\\\ 1.6\cdot
10^{-3}&3.6\cdot 10^{-5}&8.1\cdot 10^{-4}\end{pmatrix}\,.$ (52)

##  5 Global fit bounds for the general neutrino case (G-SS)

When the seesaw scenario features more than three neutrinos, the dim-5 and 6
operators are in general independent (see e.g. Refs. [114, 58, 115, 116]) and
therefore there are no correlations in the flavor structure of the mixing
coming from the correct reproduction of neutrino masses and mixings. The
Schwarz inequality will in general not be saturated either, in contrast to the
previous scenarios. Thus, the 6 elements of the $\eta$-matrix are completely
independent and cLFV observables cannot be used to constrain the LFC ones.
Therefore, we can extract the bounds on the diagonal elements
$\eta_{\alpha\alpha}$ via a global fit to the observables of Table 1 and,
separately, extract the bounds on the off-diagonal entries from the bounds on
the cLFV observables using Eq. (35). Additionally, the Schwarz inequality
given by Eq. (7) also allows us to derive constraints on the off-diagonal
elements from the bounds on the diagonal ones.

Figure 5:  Same as Fig. 2 but for the general neutrino scenario with an
arbitrary number of RH neutrinos (G-SS).

The results of our global fit are shown in Fig. 5 and Table 6. As in previous
sections, Fig. 5 shows the profiles obtained following Wilks’ theorem and
those obtained after calibrating the test statistic. Note that, despite the
absence of correlations induced by reproducing the observed light neutrino
mass and mixing, the $\eta_{\alpha\alpha}\geq 0$ condition still holds, as
required by Eq. (5), which imposes a physical boundary and violates the
requirements for Wilks’ theorem to apply. Therefore, we could expect that
calibration effects will be more important close to the
$\eta_{\alpha\alpha}=0$ border. Indeed, this can be clearly appreciated in the
$\eta_{ee}$ profile, where the most substantial deviations appear near the
physical boundary at zero. As we move away from $\eta_{ee}=0$, the two
profiles converge, as naively expected. However, when moving away from the
boundary, the $\eta_{\mu\mu}$ and $\eta_{\tau\tau}$ profiles do not seem to
converge to Wilks’. This is because $\eta_{\mu\mu}$ and $\eta_{\tau\tau}$ have
a best-fit at zero and, thus, comparatively larger values would be required in
order to recover the Wilks’ behaviour.

The profiles show a slight preference ($\approx 2\sigma$) for a non-zero value
of $\eta_{ee}$, whereas $\eta_{\mu\mu}$ has quite stringent bounds. This is
again driven by the Cabibbo anomaly in the CKM sector tightly constraining the
size of $\eta_{\mu\mu}$ and several observables pushing for a non-zero
$\eta_{ee}$. Additionally, $\eta_{\tau\tau}$ is less tightly constrained than
$\eta_{\mu\mu}$, mainly because the $\tau$ sector is constrained by less
observables and does not contribute to the CKM anomaly. The preference for a
non-zero $\eta_{ee}$ induces also a preference for a non-zero
$\text{Tr}\left[\eta\right]$.

Regarding the off-diagonal elements, we show their corresponding profiles in
App. A, and collect their bounds in Table 6. The second and third columns
collect the bounds from the global fit to LFC observables including the
Schwarz inequality. In the fourth and fith columns we show the constraints on
the off-diagonal elements directly imposed by cLFV processes. We have
highlighted in bold the strongest bound for each sector. As can be seen, the
cLFV bound is the most constraining in the $e$-$\mu$ sector due to the
stringent upper limits on $\mu\to e$ transitions. Conversely, in the
$e$-$\tau$ and $\mu$-$\tau$ sectors, the indirect LFC bounds from the Schwarz
inequality clearly dominate. Notice the very strong $68\%$CL bound on
$\left\lvert\eta_{\mu\tau}\right\rvert$. This is a consequence of the strong
border effects that appear since both $\eta_{\mu\mu}$ and $\eta_{\tau\tau}$
have their best-fit at $0$. As expected, the $95\%$CL bound is, in turn, much
less tight.

G-SS | LFC Bound | LFV Bound
---|---|---
|  $68\%$CL |  $95\%$CL |  $68\%$CL |  $95\%$CL
$\eta_{ee}$ | $[0.33,1.0]\cdot 10^{-3}$ | $[0.081,1.4]\cdot 10^{-3}$ | - | -
$\eta_{\mu\mu}$ | $1.5\cdot 10^{-5}$ | $1.4\cdot 10^{-4}$ | - | -
$\eta_{\tau\tau}$ | $1.6\cdot 10^{-4}$ | $8.9\cdot 10^{-4}$ | - | -
$\text{Tr}\left[\eta\right]$ | $[0.28,1.2]\cdot 10^{-3}$ | $2.1\cdot 10^{-3}$ | - | -
$\left\lvert\eta_{e\mu}\right\rvert$ | $1.4\cdot 10^{-4}$ | $3.4\cdot 10^{-4}$ | $\mathbf{8.4\cdot 10^{-6}}$ | $\mathbf{1.2\cdot 10^{-5}}$
$\left\lvert\eta_{e\tau}\right\rvert$ | $\mathbf{4.2\cdot 10^{-4}}$ | $\mathbf{8.8\cdot 10^{-4}}$ | $5.7\cdot 10^{-3}$ | $8.1\cdot 10^{-3}$
$\left\lvert\eta_{\mu\tau}\right\rvert$ | $\mathbf{9.4\cdot 10^{-6}}$ | $\mathbf{1.8\cdot 10^{-4}}$ | $6.6\cdot 10^{-3}$ | $9.4\cdot 10^{-3}$
Table 6: Upper bounds (or preferred intervals) for the G-SS. The LFC bounds
are obtained from the global fit analysis to the observables in Table 1, in
particular from the boostrapped profiles in Fig. 6 (see also Fig. 7 in the
appendix). For off-diagonal $\eta_{\alpha\beta}$ elements, we also derive
limits from cLFV transitions and highlight the strongest bound for each flavor
sector.

As before, we can translate the bounds to the $\alpha$ parametrization, which
at $95\%$CL reads

$\left\lvert\mathbb{I}-\alpha\right\rvert=\begin{pmatrix}[0.081,1.4]\cdot
10^{-3}&0&0\\\ <2.4\cdot 10^{-5}&<1.4\cdot 10^{-4}&0\\\ <1.8\cdot
10^{-3}&<3.6\cdot 10^{-4}&<8.9\cdot 10^{-4}\end{pmatrix}.$ (53)

Finally, also in this more general setup the $\eta$ parametrization can be
connected to the mixing between active and RH neutrinos, as given in Eq. (5).
Consequently, the bounds in Table 6 can be understood as bounds on the total
mixings to a given flavor, defined as the sum of (squared) mixings of each RH
neutrinos to a given flavor:

$\eta_{\alpha\alpha}=\frac{1}{2}|\Theta_{\alpha}|^{2}\,,\qquad{\rm
with}\quad|\Theta_{\alpha}|^{2}\equiv\sum_{k=1}^{n}|\Theta_{\alpha k}|^{2}\,.$
(54)

Thus, the bounds on these total mixings $\Theta_{\alpha}$ can be compared to
those from other experimental facilities directly searching for heavy
neutrinos (see e.g. Ref. [14, 13]), bearing in mind that they depend on the RH
neutrino mass scale and apply to lighter scales than the ones considered here.
The only caveat is that most of these experimental bounds for a given flavor
are obtained switching-off the mixing with the other flavors, while we
marginalized over them. Nevertheless, we have verified than considering the
simplified scenario in which there is mixing exclusively with one given flavor
leads to very similar results to those in Table 6.

##  6 Global fit bounds for generic unitarity violation (GUV)

Finally, we consider the case in which the $\eta$-matrix is not assumed to be
positive-definite. Thus, its diagonal entries are allowed to be also negative
while the off-diagonal parameters are not required to satisfy the Schwarz’s
inequality. This generic unitarity violation (GUV) cannot be realised in a
pure Type-I Seesaw mechanism and a more elaborate particle content, with
different contributions to non-unitarity would be required [55]. This extra
freedom can lead to rather different results from the G-SS setup due to the
Cabibbo anomaly being solved by negative values of $\eta_{\mu\mu}$, which are
unattainable for a Type-I Seesaw but are possible for an unconstrained
$\eta$-matrix [55].

From the point of view of the analysis, the LFC and LFV fits can be decoupled
as the diagonal and off-diagonal elements are not correlated, similarly to the
G-SS. However, contrary to the G-SS case, the off-diagonal entries can only be
constrained from LFV bounds, as the Schwarz’s inequality no longer holds.
Furthermore, in the GUV setup there are no physical boundaries and the
$\eta_{\alpha\alpha}$ can cover both positive and negative values. This means
that Wilks’ theorem can be safely assumed, as the violations present in the
G-SS are no longer an issue.

It should be noted that, contrary to the previous cases, in this scenario we
have adopted an agnostic stance on the kind of new physics responsible for
non-unitarity. As such, the extra degrees of freedom that may contribute to
the loop level cLFV observables are unknown. We thus consider exclusively the
contribution of the light neutrinos to the cLFV processes in order to obtain a
bound on $\left\lvert\eta_{\alpha\beta}\right\rvert$, which assumes that no
substantial cancellations with the contributions from possible new particles.
In particular, the contribution of only the light neutrinos to radiative
decays is:

${\rm
BR}(\ell_{\alpha}\to\ell_{\beta}\gamma)=\frac{25\alpha}{6\pi}\,\big{|}\eta_{\alpha\beta}\big{|}^{2}\,,$
(55)

as opposed to Eq. (35) which contains also the heavy neutrino contributions.
Since we are neglecting the contributions from the new particles, the off-
diagonal bounds we will quote using Eq. (55) are only orientative and model-
dependent.

We perform our statistical analysis for the LFC observables leaving the sign
of $\eta_{\alpha\alpha}$ unconstrained. Remarkably, a much better fit to the
data than in the G-SS is found. Indeed, while the G-SS improves the fit with
respect to the SM by $\Delta\chi^{2}=3.75$, the GUV setup improves upon the SM
fit by $\Delta\chi^{2}=11.07$. As previously stated, this is mainly due to the
fact that the GUV scenario can fit the Cabibbo anomaly, whereas the G-SS
cannot.

Figure 6: Same as Fig. 2 but for the generic unitarity violation scenario
(GUV), where the $\eta$-matrix can also take negative values and Eq. (7) does
not need to be satisfied.

The resulting profiles are shown in Fig. 6, where bootstrapping provides no
deviation from Wilks’ theorem, as expected. Apart from the $2\sigma$
preference for a non-zero $\eta_{ee}$, which was already present in the G-SS,
we also find a $2\sigma$ preference for $\eta_{\mu\mu}<0$ induced by the
Cabibbo anomaly. Conversely, we do not find any significance preference for
$\eta_{\tau\tau}$ different from zero. The corresponding preferred intervals
and upper bounds are collected in Table 7.

Finally, as in the previous cases, these bounds can be translated to the
$\alpha$ parametrization. At $95\%$CL, we have:

$\displaystyle 1-\alpha_{11}$ $\displaystyle\in[0.20,1.65]\cdot 10^{-3}\,,$
$\displaystyle\left\lvert\alpha_{21}\right\rvert$ $\displaystyle<1.4\cdot
10^{-5}\,,$ $\displaystyle 1-\alpha_{22}$ $\displaystyle\in[-1.1,-0.088]\cdot
10^{-3}\,,$ $\displaystyle\left\lvert\alpha_{31}\right\rvert$
$\displaystyle<9.8\cdot 10^{-3}\,,$ (56) $\displaystyle 1-\alpha_{33}$
$\displaystyle\in[-3.1,0.56]\cdot 10^{-3}\,,$
$\displaystyle\left\lvert\alpha_{32}\right\rvert$ $\displaystyle<1.1\cdot
10^{-2}\,.$

GUV | LFC Bound |  | LFV Bound
---|---|---|---
|  $68\%$CL |  $95\%$CL |  |  $68\%$CL |  $95\%$CL
$\eta_{ee}$ | $[0.56,1.29]\cdot 10^{-3}$ | $[0.20,1.65]\cdot 10^{-3}$ | $\left\lvert\eta_{e\mu}\right\rvert$ | $5.0\cdot 10^{-6}$ | $7.2\cdot 10^{-6}$
$\eta_{\mu\mu}$ | $[-8.2,-3.3]\cdot 10^{-4}$ | $[-1.1,-0.088]\cdot 10^{-3}$ | $\left\lvert\eta_{e\tau}\right\rvert$ | $3.4\cdot 10^{-3}$ | $4.9\cdot 10^{-3}$
$\eta_{\tau\tau}$ | $[-2.2,-0.38]\cdot 10^{-3}$ | $[-3.1,0.56]\cdot 10^{-3}$ | $\left\lvert\eta_{\mu\tau}\right\rvert$ | $4.0\cdot 10^{-3}$ | $5.6\cdot 10^{-3}$
Table 7: Upper bounds (or preferred intervals) for the GUV.

##  7 Discussion and Conclusions

In this work we have updated and improved upon present constraints on the
unitarity of the leptonic mixing matrix and the mixing of heavy right-handed
neutrinos with the SM active flavours with a global fit to flavour and
electroweak precision observables. Besides updating all the experimental
constraints from the different observables under consideration, compared to
previous studies we have improved the analysis in several additional ways. In
particular, given the expected deviations from Wilks’ theorem, we have
explicitly calibrated our test statistics by boostrapping in order to properly
estimate the significance of the constraints placed. We have also explicitly
computed the SM prediction for the different observables in terms of the input
data instead of using the results of the available electroweak fits, as has
been done in the past, so as to make the global fit fully consistent.

We have provided results for four distinct scenarios. The first three
correspond to leptonic mixing unitarity deviations induced by scenarios with
2, 3 or a higher arbitrary number of extra right-handed neutrinos, dubbed
2N-SS, 3N-SS, and G-SS, respectively. The results are summarized in Tables 4,
5, and 6. The main difference among these scenarios is the level of
correlation between the parameters that describe the unitarity deviations,
provided that the correct pattern of neutrino masses and mixings is also
recovered. The deviations from unitarity are encoded in all generality through
a small Hermitian matrix $\eta$. This matrix corresponds to the coefficient of
the only $d=6$ operator obtained at tree level upon integrating out the heavy
neutrinos, and is moreover directly connected with the square of their mixing
with the active neutrino flavours. Additionally, we also considered the
possibility of lepton mixing unitarity deviations not necessarily induced by
the presence of additional right-handed neutrinos, whose results are collected
in Table 7. This scenario removes the requirement that they are described by a
positive definite matrix $\eta$. For all the scenarios considered, we also
report our results in an alternative parametrization of unitarity deviations
through a lower triangular matrix, more appropriate when studying the neutrino
oscillation phenomenon.

We derived our results from a global fit to the observables summarized in
Table 1. Thus, the set of bounds derived is valid as long as the new degrees
of freedom are heavier than the mass of the $Z$. For lighter new particles,
some of the LEP constraints at the $Z$ pole are lost, but the remaining
observables apply down to the mass of the $\tau$. Notice that we did not
include the latest CDF-II measurement on $M_{W}$ among our observables since
it is in significant tension (above $5\sigma$) with the rest of measurements
with the same parametric dependence on $\eta$ (see section 2.5).

The strongest constraint on unitarity deviations comes from the very stringent
bound from $\mu\to e$ transitions. In the 2N-SS scenario this observable
dominates the constraints on essentially all elements, given the little
freedom and strong correlations implied by having only 2 right-handed
neutrinos reproducing the measured pattern of light neutrino masses and
mixings. In the 3N-SS, $\mu\to e$ processes instead forbids mixing to both
electrons and muons simultaneously, while in the more general scenarios it
provides strong constraints on the $\eta_{e\mu}$ element. Conversely,
constraints from $\tau\to e$ and $\tau\to\mu$ transitions are generally
subleading and stronger bounds are implied through the LFC observables.
Nevertheless, when unitarity deviations are not sourced by right-handed
neutrinos exclusively, both sets of constraints would, a priori, be
uncorrelated cLFV observables are the only bounds applying on $\eta_{e\tau}$
and $\eta_{\mu\tau}$.

Regarding the lepton flavor conserving observables we find the following
behaviours: Due to the existing tension in the unitarity test of the first row
of the CKM (the so-called Cabibbo anomaly), $\eta_{\mu\mu}>0$ (or equivalently
right-handed neutrino mixing with muons $\theta_{\mu}\neq 0$, if that is the
source of $\eta$) is disfavoured, since its presence worsens the anomaly.
Moreover, even though LFU observables are generally in good agreement with the
SM, there is a mild preference ($\sim 1\sigma$) for $\eta_{ee}>\eta_{\mu\mu}$,
both in $R_{\mu e}^{\pi}$ and $R_{\mu e}^{\tau}.$ On the other hand, the
measured values of $s^{2}_{\text{eff}}$ and $M_{W}$ show a slight preference
($\sim 1$-$2\sigma$) for a deviation from unitarity either in the electron or
muon sector.

The combination of these effects leads to constraints on all the $\eta$
parameters ranging between $10^{-3}$ and $10^{-5}$ at $2\sigma$ for both the
3N-SS and G-SS, with the exception of a preference for non-unitarity at the
level of $\eta_{ee}\sim 10^{-3}$ at around $2\sigma$ (see Figs. 4 and 5). This
implies a $2\sigma$ preference for a mixing of the heavy neutrinos with the
electron at the $|\theta_{e}|\sim 10^{-2}$ level. Conversely, unitarity
deviations are very disfavoured in the case of the muons due to the Cabibbo
anomaly and the LFU constraints. Finally, the tau sector is almost exclusively
constrained by the Z-pole observables $\Gamma_{\text{inv}}$, $\Gamma_{Z}$ and
$\sigma^{0}_{\text{had}}$ which show no preference for non-unitarity.

The main difference regarding the 2N-SS is that stronger bounds on the $\eta$
elements, ranging approximately between $10^{-4}$ and $10^{-5}$ at $2\sigma$,
are found given its more constrained structure. Indeed, for normal neutrino
mass ordering there is not enough freedom to fit the preference for a dominant
mixing with the electron flavor. Thus, slightly stronger constraints and no
preference for unitarity violations are found with respect to the inverted
ordering case. Conversely, for inverted ordering (or when more than 2 right-
handed neutrinos are considered), the constraints implied by the measured
pattern of light neutrino masses and mixings are compatible with a dominant
role of the mixing to the electron flavor, as data prefers. Thus, these
scenarios provide a somewhat better fit with a preference for unitarity
violation in the electron sector.

Finally, if the unitarity deviations are not assumed to be positive definite
(as required when they are solely induced by right-handed neutrinos), negative
values in the muon sector may actually solve the Cabibbo anomaly. Thus, in
this most general, but also more complex, extension a significantly better fit
is found with a preference for unitarity deviations ($\eta_{\alpha\alpha}\neq
0$) at the $2\sigma$ level both in the electron and muon sectors and with
opposite signs.

#### Acknowledgments.

The authors thank Gonzalo Morrás for very illuminating discussions. This
project has received support from the European Union’s Horizon 2020 research
and innovation programme under the Marie Skłodowska-Curie grant agreement No
860881-HIDDeN and No 101086085 - ASYMMETRY, and from the Spanish Research
Agency (Agencia Estatal de Investigación) through the Grant IFT Centro de
Excelencia Severo Ochoa No CEX2020-001007-S and Grant PID2019-108892RB-I00
funded by MCIN/AEI/10.13039/501100011033. EFM, XM and DNT acknowledge support
from the HPC-Hydra cluster at IFT. XM acknowledges funding from the European
Union’s Horizon Europe Programme under the Marie Skłodowska-Curie grant
agreement no. 101066105-PheNUmenal. The work of DNT was supported by the
Spanish MIU through the National Program FPU (grant number FPU20/05333). JLP
also acknowledges support from Generalitat Valenciana through the plan GenT
program (CIDEGENT/2018/019) and from the Spanish Ministerio de Ciencia e
Innovacion through the project PID2020-113644GB-I00. JHG warmly thanks the
hospitality of Albert De Roeck and the EP Neutrino group during his stay at
CERN; where this project has been completed.

##  Appendix A Details of the analysis

Wilks’ theorem is a common assumption in many statistical analyses. It
establishes that the test statistic on the form

$\Delta\chi^{2}=\chi^{2}(\theta)-\min_{\theta^{\prime}}\chi^{2}(\theta^{\prime})\,,$
(57)

where

$\chi^{2}(\theta)=\sum_{d\in
D}\left(\frac{\mu_{d}(\theta)-\bar{\mu}_{d}}{\sigma_{d}}\right)^{2},$ (58)

$\theta$ is a set of model parameters, $D$ the set of data points,
$\mu_{d}(\theta)$ the expected data for $d$ given $\theta$, $\bar{\mu}_{d}$
the observed data, and $\sigma_{d}^{2}$ the variance in $\mu_{d}(\theta)$, is
everywhere distributed according to a $\chi^{2}$-distribution, with a number
of degrees of freedom that can be derived from the total number of free
parameters. When the theorem applies, it can be used to directly map the test
statistic to a confidence level.

However, for the theorem to hold, several assumptions on the behaviour of the
test statistic have to be made. One of the most important assumptions, which
is generally violated in our analysis, is that the expectation
$\mu_{d}(\theta)$ for the observables included in the fit must describe a
hyperplane in the space of possible observations as one varies the model
parameters $\theta$. A straight-forward source of violation of Wilks’ theorem
assumptions is that our parameters have a physical boundary
($\eta_{\alpha\alpha}\geq 0$) in the 2N-SS, 3N-SS and G-SS scenarios, and
therefore so will our observables. Moreover, we will also find violations when
including LFV observables, since they depend quadratically on $\eta$ unlike to
the LFC ones of Table 1, which depend linearly on $\eta$, thereby leading to a
curved expected region of the observations. Additionally, when the observables
depend on cyclic parameters, such as phases, they will describe compact
trajectories in the space of possible observations, providing and additional
source of violation of Wilks’ theorem.

Our analysis is therefore plagued with potential deviations from Wilks’
theorem, and we cannot rely on it for extracting confidence intervals.
Instead, the test statistic must be calibrated at each point of parameter
space, in order to know how it is distributed and associate the resulting
value of $\Delta\chi^{2}$ to a confidence level. This procedure is commonly
referred to as bootstrapping. We describe the bootstrapping procedure of our
analysis below.

We start by computing the $\chi^{2}(\eta_{\alpha\beta})$ where
$\eta_{\alpha\beta}$ are the parameters of the model assuming Gaussian
uncertainties for all the experimental values of the observables $O_{i}$
reported in Table 1 and taking into account correlations among them when
relevant:

$\chi^{2}(\eta_{\alpha\beta})=\sum_{i,j}\frac{\big{(}O_{i}-E_{i}(\eta_{\alpha\beta})\big{)}}{\sigma_{i}}\mathrm{corr}_{ij}\frac{\big{(}O_{j}-E_{j}(\eta_{\alpha\beta})\big{)}}{\sigma_{j}}\,,$
(59)

where $E_{i}(\eta_{\alpha\beta})$ is the expectation of $O_{i}$ given the
model parameters $\eta_{\alpha\beta}$ as given by the first two columns in
Table 1, $\sigma_{i}$ is the corresponding uncertainty reported together with
$O_{i}$ in the third column of Table 1 and $\mathrm{corr}_{ij}$ is the
correlation matrix among the different observables.

Since we are interested on showing the CL profiles for a particular parameter
of interest, for definiteness $\eta_{ee}$ but the procedure would be the same
for any other parameter, we choose the profiled $\chi^{2}$ over all the
parameters $\eta_{\alpha\beta}$ except for $\eta_{ee}$ itself as our test
statistic. Denoting the other (nuisance) parameters
$\vec{\nu}=\eta_{\alpha\beta}$ with $\alpha\beta\neq ee$, then the profiled
$\chi^{2}$ is:

$\hat{\chi}^{2}(\eta_{ee})\equiv\min_{\vec{\nu}}\chi^{2}(\eta_{ee},\vec{\nu})\,.$
(60)

If Wilks’ theorem holds, then the significance with which a point
$\eta_{ee}^{0}$ may be excluded by the observables $O_{i}$ will be given by
the square root of:

$\Delta\chi^{2}(\eta_{ee}^{0})=\hat{\chi}^{2}(\eta_{ee}^{0})-\chi^{2}_{\text{min}}\,,$
(61)

where $\chi^{2}_{\text{min}}$ is the minimum value of the test statistic,
since $\hat{\chi}^{2}(\eta_{ee})$ would follow a $\chi^{2}$ distribution with
1 degree of freedom. Note that only the value of $\eta_{ee}^{0}$ is relevant
independently of in which point of nuisance parameter space $\vec{\nu}$ the
statement is made, as Wilks’ theorem guarantees that the test statistic of Eq.
(61) is $\chi^{2}$-distributed (with one degree of freedom) in all the
parameter space. However, in the presence of violations of Wilks’ theorem,
there is no reason why this should be the case and we cannot assume a
$\chi^{2}$-distribution for the test statistic. Moreover, the distribution of
the test statistic may be different in different points of parameter space.

Instead, in order to claim an exclusion confidence level for $\eta_{ee}^{0}$,
we need to consider all points in the $\vec{\nu}$ parameter space, calibrate
the test statistic’s distribution in each of these points, and then compare
each of these calibrated distributions with the actual value $\Delta\chi^{2}$
from Eq. (61) in order to extract a confidence level. The minimum of these
confidence levels would then represent the significance of the exclusion of
$\eta_{ee}^{0}$ regardless of $\vec{\eta}$.

We therefore implement the following bootstrap procedure:

  1. 1.

The values of the parameter of interest for which we want to quote a
confidence level are fixed, $\\{\eta_{ee}^{i}\\}_{i=1}^{n}$. For each of those
values, we consider a grid of points in $\vec{\nu}$ space.

  2. 2.

For each of the $\eta_{ee}^{i}$, we consider the different $\vec{\nu}$ and
compute the corresponding predictions for the observables of the fit.

  3. 3.

Using the predictions as the central value, we generate pseudo-data for the
observables by drawing from a Gaussian distribution.

  4. 4.

For each of the sampled sets of pseudo-data, we compute the corresponding
value of the test statistic $\Delta\chi^{2}(\eta_{ee}^{i})$.

Repeating the generation of pseudo-data and consequent computation of the test
statistic many times results in a large sample of the test statistic
distribution. A confidence level is then computed by comparing the real data
with this distribution, the CL simply being the percentage of times that
$\Delta\chi^{2}(\eta_{ee}^{i})$ of the pseudo-data sets are smaller than the
value computed for the real data. Repeating this procedure for all $\vec{\nu}$
points for a certain $\eta_{ee}^{i}$ value and then picking the minimum CL, we
can finally associate an exclusion level to $\eta_{ee}^{i}$.

Since we need to generate a large amount of pseudo-data over a grid of points
in parameter space, the number of points required grows exponentially with the
number of parameters. Moreover, since at each point of the grid many numerical
minimizations of the test statistic are required in order to reconstruct its
distribution, the computational cost of this procedure quickly becomes
unfeasible. Nevertheless, for the G-SS and 2N-SS scenarios, this procedure is
still computationally tractable, since their parameter spaces are three-
dimensional. However, in the case of the 3N-SS, characterized by 8 free
parameters, this is no longer the case.

We have therefore opted for an approximation to the bootstrap procedure for
the 3N-SS scenario, which we dub Profiled Bootstrap. As it name suggests, it
consists on only performing the calibration in the points that profile the
test statistic along the direction of some parameter of interest $\eta_{ee}$.
As such, for each value $\eta_{ee}^{i}$, instead of considering a grid of
points in nuisance parameter space $\vec{\nu}$, we will only consider the
point:

$\hat{\vec{\nu}}(\eta_{ee}^{i})\equiv\underset{\vec{\nu}}{\text{argmin
}}\chi^{2}(\eta_{ee}^{i},\vec{\nu})\,.$ (62)

We are thus reducing the task of exploring a grid of points in $\vec{\nu}$
space to just a single point. The rest of the procedure follows steps 2-4 as
before.

The requirement that needs to be satisfied in order for the profiled bootstrap
to be an acceptable approximation to the complete bootstrap is that the point
$\hat{\vec{\nu}}(\eta_{ee}^{i})$ needs to be close to the point in $\vec{\nu}$
space that yields the minimum CL for a fixed $\eta_{ee}^{i}$. While this can
seem to intuitively always be the case, there may be scenarios in which points
in the parameter space with a value of the test statistic far from its minimum
may yield a smaller CL if the distribution of the test statistic leans toward
even larger values in that point of the parameter space. We have thus compared
the original bootstrap and the profiled bootstrap in the G-SS case, which is a
scenario in which both procedures are tractable, to test its accuracy. The
results are shown in the upper panels of Fig. 7 where we can see that, in
general, both the solid and dashed black lines are in good agreement.

Figure 7: Comparison between the bootstrap (solid black) and profiled
bootstrap (dashed black) for the G-SS case. We also show the expectations of
Wilks’ theorem (dashed green) for reference. In general we see quite good
agreement between both procedures for the diagonal elements in the upper
panels, but not for the off-diagonal combinations in the lower panels.

The reason why this behaviour is expected to also hold for the 3N-SS is that
both the G-SS and 3N-SS showcase very similar freedom when it comes to fitting
the data: even though the 3N-SS case presents non trivial correlations between
the $\eta_{\alpha\alpha}$ parameters (see Eq. (50)), they do not strongly
constrain the flavor pattern999This is not the case of the 2N-SS, which
showcases a quite predictive flavor pattern. However, it should be noted the
flavor pattern of the 3N-SS becomes more predictive for $m_{\text{lightest}}$
values much smaller than the existing bounds [54]., as it can be seen in Fig.
3.

However, while the agreement between these two procedures is good for the
profiles of the diagonal elements $\eta_{\alpha\alpha}$, this is no longer the
case for the off-diagonal combinations
$\sqrt{\eta_{\alpha\alpha}\eta_{\beta\beta}}$, as can be seen in the lower
panels of Fig. 7. Therefore, assuming the same discrepancy will be present for
the 3N-SS, we adopted a more conservative approach and the off-diagonal bounds
quoted in Table 5 are extracted assuming Wilks’ theorem, which provides the
weaker constraint between the two options. Note that the latter only applies
for the 3N-SS, since in the 2N-SS and G-SS scenarios we are able to follow the
full bootstrap procedure also for the off-diagonal elements.

##  Appendix B About cancellations in the cLFV rates

In section 2.6 we neglected potential differences between the masses of all
heavy neutrinos, so we could simplify the discussion of cLFV transitions. The
bounds we derived after including these rates should still apply to more
generic scenarios, nevertheless in principle having several mass scales could
lead to potential cancellations in some of the cLFV rates, as advocated in
Ref. [105]. We devote this appendix to such a discussion.

First of all, let us emphasize that this discussion does not apply to the
minimal scenarios in sections 3 and 4, where the Schwarz inequality in Eq. (7)
is saturated, meaning that the cLFV bounds must be included in our global fit
and that they impose bounds also for the diagonal entries of the
$\eta$-matrix. The reason is that there is only a single relevant heavy scale
and, moreover, the flavor mixing pattern is quite defined by oscillation data.
Thus, there is not enough freedom to accommodate the potential cancellations
we discuss here, and the bounds derived for those scenarios are robust.

The question about potential cancellations in the cLFV rates becomes relevant
for models with more heavy neutrinos, for which the general bounds in section
5 should be respected, and when having several heavy neutrino mass scales. In
such a scenario, and taking the radiative decays as an example, we have

$\Gamma(\ell_{\alpha}\to\ell_{\beta}\gamma)\propto\left|\sum_{k}\Theta_{\alpha
k}\Theta^{*}_{\beta
k}\,F_{\gamma}\Big{(}M_{N_{k}}^{2}/M_{W}^{2}\Big{)}\right|^{2}\,,$ (63)

where the sum goes over all the heavy neutrinos, and the missing factors and
explicit form of the loop-function $F_{\gamma}$ can be found for instance in
Ref. [38].

If all heavy neutrinos are degenerate or, more generally, if the value of
$F_{\gamma}\Big{(}M_{N_{k}}^{2}/M_{W}^{2}\Big{)}$ is common for all
contributions, the sum of $\Theta_{\alpha k}\Theta^{*}_{\beta k}$ gives the
$\eta_{\alpha\beta}$ dependence in Eq. (5). Otherwise, additional dependence
on other parameters beyond $\eta_{\alpha\beta}$ will be present, potentially
even leading to a suppression of the radiative decay for non-vanishing
$\eta_{\alpha\beta}$. Nevertheless, the loop function monotonically increases
from $1/8$ to $1/2$ as the heavy neutrino mass varies from $M_{W}$ to
infinity. Thus, deviations from a constant value for
$F_{\gamma}\Big{(}M_{N_{k}}^{2}/M_{W}^{2}\Big{)}$ cannot be large.
Furthermore, a cancellation requires opposite-sign contributions from the
mixing terms $\Theta_{\alpha k}\Theta^{*}_{\beta k}$, which would also
(partially) suppress $\eta_{\alpha\beta}$. Therefore, having a large
$\left\lvert\eta_{\alpha\beta}\right\rvert$ with suppressed
$\ell_{\alpha}\to\ell_{\beta}\gamma$ is possible only for rather large
$\Theta_{\alpha k}\Theta^{*}_{\beta k}$ mixings and for very specific values
of the heavy neutrino masses, so that is only a partial cancellation in
$\left\lvert\eta_{\alpha\beta}\right\rvert$ but the values of $F_{\gamma}$ are
such that the radiative decays are strongly suppressed.

Furthermore, this kind of configurations would lead to a numerical
cancellation only for the radiative decays, but in general they will not
suppress other cLFV observables such as the 3-body decays or $\mu-e$
conversion in nuclei, which can actually be more restrictive, as we have seen
in section 2.6. While adding more neutrinos enhances the number of parameters
and the freedom to choose very specific values of mixings and masses so that
all the cLFV channels are suppressed (but not the respective
$\left\lvert\eta_{\alpha\beta}\right\rvert$), this seems an extremely unlikely
configuration.

Nevertheless, to account also for this possibility, our results for the G-SS
in section 5 are given considering separately the lepton flavor conserving and
violating channels, see Table 6. While for
$\left\lvert\eta_{e\mu}\right\rvert$ the bounds from the cLFV observables are
stronger, they may be regarded as more model dependent, as discussed above.
Conversely, LFC bounds, which are also the most stringent for
$\left\lvert\eta_{e\tau}\right\rvert$ and
$\left\lvert\eta_{\mu\tau}\right\rvert$ cannot be avoided by these kind of
cancellations.

## References

  * [1] S. Dodelson and L. M. Widrow, Sterile-neutrinos as dark matter, Phys. Rev. Lett. 72 (1994) 17–20, [hep-ph/9303287].
  * [2] X.-D. Shi and G. M. Fuller, A New dark matter candidate: Nonthermal sterile neutrinos, Phys. Rev. Lett. 82 (1999) 2832–2835, [astro-ph/9810076].
  * [3] K. Abazajian, G. M. Fuller, and M. Patel, Sterile neutrino hot, warm, and cold dark matter, Phys. Rev. D 64 (2001) 023501, [astro-ph/0101524].
  * [4] T. Asaka, S. Blanchet, and M. Shaposhnikov, The nuMSM, dark matter and neutrino masses, Phys. Lett. B 631 (2005) 151–156, [hep-ph/0503065].
  * [5] A. Falkowski, J. Juknevich, and J. Shelton, Dark Matter Through the Neutrino Portal, arXiv:0908.1790.
  * [6] M. Lindner, A. Merle, and V. Niro, Enhancing Dark Matter Annihilation into Neutrinos, Phys. Rev. D 82 (2010) 123529, [arXiv:1005.3116].
  * [7] V. Gonzalez Macias and J. Wudka, Effective theories for Dark Matter interactions and the neutrino portal paradigm, JHEP 07 (2015) 161, [arXiv:1506.03825].
  * [8] M. Blennow, E. Fernandez-Martinez, A. Olivares-Del Campo, S. Pascoli, S. Rosauro-Alcaraz, and A. V. Titov, Neutrino Portals to Dark Matter, Eur. Phys. J. C 79 (2019), no. 7 555, [arXiv:1903.00006].
  * [9] M. Fukugita and T. Yanagida, Baryogenesis Without Grand Unification, Phys. Lett. B 174 (1986) 45–47.
  * [10] E. K. Akhmedov, V. A. Rubakov, and A. Y. Smirnov, Baryogenesis via neutrino oscillations, Phys. Rev. Lett. 81 (1998) 1359–1362, [hep-ph/9803255].
  * [11] T. Asaka and M. Shaposhnikov, The $\nu$MSM, dark matter and baryon asymmetry of the universe, Phys. Lett. B 620 (2005) 17–26, [hep-ph/0505013].
  * [12] M. Drewes, B. Garbrecht, P. Hernandez, M. Kekic, J. Lopez-Pavon, J. Racker, N. Rius, J. Salvado, and D. Teresi, ARS Leptogenesis, Int. J. Mod. Phys. A 33 (2018), no. 05n06 1842002, [arXiv:1711.02862].
  * [13] C. Antel et al., Feebly Interacting Particles: FIPs 2022 workshop report, in Workshop on Feebly-Interacting Particles, 5, 2023.  arXiv:2305.01715.
  * [14] E. Fernández-Martínez, M. González-López, J. Hernández-García, M. Hostert, and J. López-Pavón, Effective portals to heavy neutral leptons, arXiv:2304.06772.
  * [15] M. Hostert. Heavy Neutrino Limits GitHub repository, https://github.com/mhostert/Heavy-Neutrino-Limits.
  * [16] S. T. Petcov, The Processes $\mu\rightarrow e+\gamma,\mu\rightarrow e+\overline{e},\nu^{\prime}\rightarrow\nu+\gamma$ in the Weinberg-Salam Model with Neutrino Mixing, Sov. J. Nucl. Phys. 25 (1977) 340. [Erratum: Sov.J.Nucl.Phys. 25, 698 (1977), Erratum: Yad.Fiz. 25, 1336 (1977)].
  * [17] S. M. Bilenky, S. T. Petcov, and B. Pontecorvo, Lepton Mixing, mu –$>$ e + gamma Decay and Neutrino Oscillations, Phys. Lett. B 67 (1977) 309.
  * [18] T.-P. Cheng and L.-F. Li, Muon Number Nonconservation in Gauge Theories, Stud. Nat. Sci. 12 (1977) 659–681.
  * [19] W. J. Marciano and A. I. Sanda, Exotic Decays of the Muon and Heavy Leptons in Gauge Theories, Phys. Lett. B 67 (1977) 303–305.
  * [20] B. W. Lee, S. Pakvasa, R. E. Shrock, and H. Sugawara, Muon and Electron Number Nonconservation in a V-A Gauge Model, Phys. Rev. Lett. 38 (1977) 937. [Erratum: Phys.Rev.Lett. 38, 1230 (1977)].
  * [21] B. W. Lee and R. E. Shrock, Natural Suppression of Symmetry Violation in Gauge Theories: Muon - Lepton and Electron Lepton Number Nonconservation, Phys. Rev. D 16 (1977) 1444.
  * [22] R. E. Shrock, New Tests For, and Bounds On, Neutrino Masses and Lepton Mixing, Phys. Lett. B 96 (1980) 159–164.
  * [23] J. Schechter and J. W. F. Valle, Neutrino Masses in SU(2) x U(1) Theories, Phys. Rev. D 22 (1980) 2227.
  * [24] R. E. Shrock, General Theory of Weak Leptonic and Semileptonic Decays. 1. Leptonic Pseudoscalar Meson Decays, with Associated Tests For, and Bounds on, Neutrino Masses and Lepton Mixing, Phys. Rev. D 24 (1981) 1232\.
  * [25] R. E. Shrock, General Theory of Weak Processes Involving Neutrinos. 2. Pure Leptonic Decays, Phys. Rev. D 24 (1981) 1275.
  * [26] P. Langacker and D. London, Mixing Between Ordinary and Exotic Fermions, Phys. Rev. D 38 (1988) 886.
  * [27] A. Pilaftsis, Lepton flavor nonconservation in H0 decays, Phys. Lett. B 285 (1992) 68–74.
  * [28] A. Ilakovac and A. Pilaftsis, Flavor violating charged lepton decays in seesaw-type models, Nucl. Phys. B 437 (1995) 491, [hep-ph/9403398].
  * [29] E. Nardi, E. Roulet, and D. Tommasini, Limits on neutrino mixing with new heavy particles, Phys. Lett. B 327 (1994) 319–326, [hep-ph/9402224].
  * [30] D. Tommasini, G. Barenboim, J. Bernabeu, and C. Jarlskog, Nondecoupling of heavy neutrinos and lepton flavor violation, Nucl. Phys. B 444 (1995) 451–467, [hep-ph/9503228].
  * [31] J. I. Illana and T. Riemann, Charged lepton flavor violation from massive neutrinos in Z decays, Phys. Rev. D 63 (2001) 053004, [hep-ph/0010193].
  * [32] W. Loinaz, N. Okamura, S. Rayyan, T. Takeuchi, and L. C. R. Wijewardhana, Quark lepton unification and lepton flavor nonconservation from a TeV scale seesaw neutrino mass texture, Phys. Rev. D 68 (2003) 073001, [hep-ph/0304004].
  * [33] E. Arganda, A. M. Curiel, M. J. Herrero, and D. Temes, Lepton flavor violating Higgs boson decays from massive seesaw neutrinos, Phys. Rev. D 71 (2005) 035011, [hep-ph/0407302].
  * [34] W. Loinaz, N. Okamura, S. Rayyan, T. Takeuchi, and L. C. R. Wijewardhana, The NuTeV anomaly, lepton universality, and nonuniversal neutrino gauge couplings, Phys. Rev. D 70 (2004) 113004, [hep-ph/0403306].
  * [35] S. Antusch, C. Biggio, E. Fernandez-Martinez, M. B. Gavela, and J. Lopez-Pavon, Unitarity of the Leptonic Mixing Matrix, JHEP 10 (2006) 084, [hep-ph/0607020].
  * [36] S. Antusch, J. P. Baumann, and E. Fernandez-Martinez, Non-Standard Neutrino Interactions with Matter from Physics Beyond the Standard Model, Nucl. Phys. B 810 (2009) 369–388, [arXiv:0807.1003].
  * [37] C. Biggio, The Contribution of fermionic seesaws to the anomalous magnetic moment of leptons, Phys. Lett. B 668 (2008) 378–384, [arXiv:0806.2558].
  * [38] R. Alonso, M. Dhen, M. B. Gavela, and T. Hambye, Muon conversion to electron in nuclei in type-I seesaw models, JHEP 01 (2013) 118, [arXiv:1209.2679].
  * [39] A. Abada, D. Das, A. M. Teixeira, A. Vicente, and C. Weiland, Tree-level lepton universality violation in the presence of sterile neutrinos: impact for $R_{K}$ and $R_{\pi}$, JHEP 02 (2013) 048, [arXiv:1211.3052].
  * [40] E. Akhmedov, A. Kartavtsev, M. Lindner, L. Michaels, and J. Smirnov, Improving Electro-Weak Fits with TeV-scale Sterile Neutrinos, JHEP 05 (2013) 081, [arXiv:1302.1872].
  * [41] L. Basso, O. Fischer, and J. J. van der Bij, Precision tests of unitarity in leptonic mixing, EPL 105 (2014), no. 1 11001, [arXiv:1310.2057].
  * [42] A. Abada, A. M. Teixeira, A. Vicente, and C. Weiland, Sterile neutrinos in leptonic and semileptonic decays, JHEP 02 (2014) 091, [arXiv:1311.2830].
  * [43] E. Arganda, M. J. Herrero, X. Marcano, and C. Weiland, Imprints of massive inverse seesaw model neutrinos in lepton flavor violating Higgs boson decays, Phys. Rev. D 91 (2015), no. 1 015001, [arXiv:1405.4300].
  * [44] S. Antusch and O. Fischer, Non-unitarity of the leptonic mixing matrix: Present bounds and future sensitivities, JHEP 10 (2014) 094, [arXiv:1407.6607].
  * [45] S. Antusch and O. Fischer, Testing sterile neutrino extensions of the Standard Model at future lepton colliders, JHEP 05 (2015) 053, [arXiv:1502.05915].
  * [46] A. Abada, V. De Romeri, S. Monteil, J. Orloff, and A. M. Teixeira, Indirect searches for sterile neutrinos at a high-luminosity Z-factory, JHEP 04 (2015) 051, [arXiv:1412.6322].
  * [47] A. Abada, V. De Romeri, and A. M. Teixeira, Impact of sterile neutrinos on nuclear-assisted cLFV processes, JHEP 02 (2016) 083, [arXiv:1510.06657].
  * [48] A. Abada and T. Toma, Electric Dipole Moments of Charged Leptons with Sterile Fermions, JHEP 02 (2016) 174, [arXiv:1511.03265].
  * [49] A. Abada and T. Toma, Electron electric dipole moment in Inverse Seesaw models, JHEP 08 (2016) 079, [arXiv:1605.07643].
  * [50] V. De Romeri, M. J. Herrero, X. Marcano, and F. Scarcella, Lepton flavor violating Z decays: A promising window to low scale seesaw neutrinos, Phys. Rev. D 95 (2017), no. 7 075028, [arXiv:1607.05257].
  * [51] E. Arganda, M. J. Herrero, X. Marcano, R. Morales, and A. Szynkman, Effective lepton flavor violating H$\ell$i$\ell$j vertex from right-handed neutrinos within the mass insertion approximation, Phys. Rev. D 95 (2017), no. 9 095029, [arXiv:1612.09290].
  * [52] M. J. Herrero, X. Marcano, R. Morales, and A. Szynkman, One-loop effective LFV $Zl_{k}l_{m}$ vertex from heavy neutrinos within the mass insertion approximation, Eur. Phys. J. C 78 (2018), no. 10 815, [arXiv:1807.01698].
  * [53] X. Marcano and R. A. Morales, Flavor techniques for LFV processes: Higgs decays in a general seesaw model, Front. in Phys. 7 (2020) 228, [arXiv:1909.05888].
  * [54] M. Chrzaszcz, M. Drewes, T. E. Gonzalo, J. Harz, S. Krishnamurthy, and C. Weniger, A frequentist analysis of three right-handed neutrinos with GAMBIT, Eur. Phys. J. C 80 (2020), no. 6 569, [arXiv:1908.02302].
  * [55] A. M. Coutinho, A. Crivellin, and C. A. Manzari, Global Fit to Modified Neutrino Couplings and the Cabibbo-Angle Anomaly, Phys. Rev. Lett. 125 (2020), no. 7 071802, [arXiv:1912.08823].
  * [56] K. A. U. Calderón, I. Timiryasov, and O. Ruchayskiy, Improved constraints and the prospects of detecting TeV to PeV scale Heavy Neutral Leptons, arXiv:2206.04540.
  * [57] S. Weinberg, Baryon and Lepton Nonconserving Processes, Phys. Rev. Lett. 43 (1979) 1566–1570.
  * [58] A. Broncano, M. B. Gavela, and E. E. Jenkins, The Effective Lagrangian for the seesaw model of neutrino mass and leptogenesis, Phys. Lett. B 552 (2003) 177–184, [hep-ph/0210271]. [Erratum: Phys.Lett.B 636, 332 (2006)].
  * [59] P. Minkowski, $\mu\to e\gamma$ at a Rate of One Out of $10^{9}$ Muon Decays?, Phys. Lett. B 67 (1977) 421–428.
  * [60] R. N. Mohapatra and G. Senjanovic, Neutrino Mass and Spontaneous Parity Nonconservation, Phys. Rev. Lett. 44 (1980) 912.
  * [61] T. Yanagida, Horizontal gauge symmetry and masses of neutrinos, Conf. Proc. C 7902131 (1979) 95–99.
  * [62] M. Gell-Mann, P. Ramond, and R. Slansky, Complex Spinors and Unified Theories, Conf. Proc. C 790927 (1979) 315–321, [arXiv:1306.4669].
  * [63] G. C. Branco, W. Grimus, and L. Lavoura, The Seesaw Mechanism in the Presence of a Conserved Lepton Number, Nucl. Phys. B 312 (1989) 492–508.
  * [64] J. Kersten and A. Y. Smirnov, Right-Handed Neutrinos at CERN LHC and the Mechanism of Neutrino Mass Generation, Phys. Rev. D 76 (2007) 073005, [arXiv:0705.3221].
  * [65] A. Abada, C. Biggio, F. Bonnet, M. B. Gavela, and T. Hambye, Low energy effects of neutrino masses, JHEP 12 (2007) 061, [arXiv:0707.4058].
  * [66] K. Moffat, S. Pascoli, and C. Weiland, Equivalence between massless neutrinos and lepton number conservation in fermionic singlet extensions of the Standard Model, arXiv:1712.07611.
  * [67] R. N. Mohapatra, Mechanism for Understanding Small Neutrino Mass in Superstring Theories, Phys. Rev. Lett. 56 (1986) 561–563.
  * [68] R. N. Mohapatra and J. W. F. Valle, Neutrino Mass and Baryon Number Nonconservation in Superstring Models, Phys. Rev. D 34 (1986) 1642\.
  * [69] E. K. Akhmedov, M. Lindner, E. Schnapka, and J. W. F. Valle, Left-right symmetry breaking in NJL approach, Phys. Lett. B 368 (1996) 270–280, [hep-ph/9507275].
  * [70] M. Malinsky, J. C. Romao, and J. W. F. Valle, Novel supersymmetric SO(10) seesaw mechanism, Phys. Rev. Lett. 95 (2005) 161801, [hep-ph/0506296].
  * [71] E. Fernandez-Martinez, M. B. Gavela, J. Lopez-Pavon, and O. Yasuda, CP-violation from non-unitary leptonic mixing, Phys. Lett. B 649 (2007) 427–435, [hep-ph/0703098].
  * [72] M. B. Gavela, T. Hambye, D. Hernandez, and P. Hernandez, Minimal Flavour Seesaw Models, JHEP 09 (2009) 038, [arXiv:0906.1461].
  * [73] E. Fernandez-Martinez, J. Hernandez-Garcia, J. Lopez-Pavon, and M. Lucente, Loop level constraints on Seesaw neutrino mixing, JHEP 10 (2015) 130, [arXiv:1508.03051].
  * [74] E. Fernandez-Martinez, J. Hernandez-Garcia, and J. Lopez-Pavon, Global constraints on heavy neutrino mixing, JHEP 08 (2016) 033, [arXiv:1605.08774].
  * [75] CDF Collaboration, T. Aaltonen et al., High-precision measurement of the $W$ boson mass with the CDF II detector, Science 376 (2022), no. 6589 170–176.
  * [76] Particle Data Group Collaboration, R. L. Workman et al., Review of Particle Physics, PTEP 2022 (2022) 083C01.
  * [77] CMS Collaboration, Precision measurement of the Z boson invisible width in pp collisions at $\sqrt{s}$ = 13 TeV, arXiv:2206.07110.
  * [78] D. Bryman, V. Cirigliano, A. Crivellin, and G. Inguglia, Testing Lepton Flavor Universality with Pion, Kaon, Tau, and Beta Decays, arXiv:2111.05338.
  * [79] HFLAV Collaboration, Y. S. Amhis et al., Averages of b-hadron, c-hadron, and $\tau$-lepton properties as of 2018, Eur. Phys. J. C 81 (2021), no. 3 226, [arXiv:1909.12524].
  * [80] M. Awramik, M. Czakon, A. Freitas, and G. Weiglein, Precise prediction for the W boson mass in the standard model, Phys. Rev. D 69 (2004) 053006, [hep-ph/0311148].
  * [81] M. Awramik, M. Czakon, and A. Freitas, Electroweak two-loop corrections to the effective weak mixing angle, JHEP 11 (2006) 048, [hep-ph/0608099].
  * [82] ATLAS Collaboration, Improved W boson Mass Measurement using 7 TeV Proton-Proton Collisions with the ATLAS Detector, .
  * [83] LHCb Collaboration, R. Aaij et al., Measurement of the W boson mass, JHEP 01 (2022) 036, [arXiv:2109.01113].
  * [84] CDF, D0 Collaboration, T. A. Aaltonen et al., Combination of CDF and D0 $W$-Boson Mass Measurements, Phys. Rev. D 88 (2013), no. 5 052018, [arXiv:1307.7627].
  * [85] ALEPH, DELPHI, L3, OPAL, LEP Electroweak Collaboration, S. Schael et al., Electroweak Measurements in Electron-Positron Collisions at W-Boson-Pair Energies at LEP, Phys. Rept. 532 (2013) 119–244, [arXiv:1302.3415].
  * [86] L. Lyons, D. Gibaut, and P. Clifford, How to combine correlated estimates of a single physical quantity, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 270 (1988), no. 1 110–117.
  * [87] A. Valassi, Combining correlated measurements of several different physical quantities, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 500 (2003), no. 1 391–405. NIMA Vol 500.
  * [88] P. Janot and S. Jadach, Improved Bhabha cross section at LEP and the number of light neutrino species, Phys. Lett. B 803 (2020) 135319, [arXiv:1912.02067].
  * [89] A. Freitas, Higher-order electroweak corrections to the partial widths and branching ratios of the Z boson, JHEP 04 (2014) 070, [arXiv:1401.2447].
  * [90] P. A. Baikov, K. G. Chetyrkin, J. H. Kuhn, and J. Rittinger, Complete ${\cal O}(\alpha_{s}^{4})$ QCD Corrections to Hadronic $Z$-Decays, Phys. Rev. Lett. 108 (2012) 222003, [arXiv:1201.5804].
  * [91] C.-Y. Seng, M. Gorchtein, H. H. Patel, and M. J. Ramsey-Musolf, Reduced Hadronic Uncertainty in the Determination of $V_{ud}$, Phys. Rev. Lett. 121 (2018), no. 24 241804, [arXiv:1807.10197].
  * [92] A. Czarnecki, W. J. Marciano, and A. Sirlin, Radiative Corrections to Neutron and Nuclear Beta Decays Revisited, Phys. Rev. D 100 (2019), no. 7 073008, [arXiv:1907.06737].
  * [93] C.-Y. Seng, X. Feng, M. Gorchtein, and L.-C. Jin, Joint lattice QCD–dispersion theory analysis confirms the quark-mixing top-row unitarity deficit, Phys. Rev. D 101 (2020), no. 11 111301, [arXiv:2003.11264].
  * [94] Y. Grossman, E. Passemar, and S. Schacht, On the Statistical Treatment of the Cabibbo Angle Anomaly, JHEP 07 (2020) 068, [arXiv:1911.07821].
  * [95] M. Blennow, P. Coloma, E. Fernández-Martínez, and M. González-López, Right-handed neutrinos and the CDF II anomaly, Phys. Rev. D 106 (2022), no. 7 073005, [arXiv:2204.04559].
  * [96] M. Maltoni and T. Schwetz, Testing the statistical compatibility of independent data sets, Phys. Rev. D 68 (2003) 033020, [hep-ph/0304176].
  * [97] S. L. Glashow, J. Iliopoulos, and L. Maiani, Weak Interactions with Lepton-Hadron Symmetry, Phys. Rev. D 2 (1970) 1285–1292.