A frequentist analysis of three right-handed neutrinos with GAMBIT

TL;DR

This paper conducts the first comprehensive frequentist analysis of a Standard Model extension with three right-handed neutrinos, covering masses from 60 MeV to 500 GeV, integrating diverse experimental constraints.

Contribution

It provides the most complete assessment to date of the three right-handed neutrino model below the TeV scale using GAMBIT, including parameter constraints and flavor mixing insights.

Findings

Derived combined profile likelihoods for key parameters.

Identified statistically preferred regions in parameter space.

Explored flavor mixing patterns for different lightest neutrino masses.

Abstract

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

Tables

Table 1. Table 1 : Electroweak precision observables measurements and uncertainties, taken from Ref.  PDG17 .

Observable	Value
Input parameters
$G_{\mu }$ [GeV-2]	$1.1663787\left(6\right)\times {10}^{-5}$
$m_Z$ [GeV]	$91.1875\left(21\right)$
Constraints
$m_W$ [GeV]	$80.385\left(15\right)$
$s_{eff}^2$	$0.23155\pm 0.00005$
${\Gamma }_{\mathrm{inv}}$ [MeV]	$499.0\pm 1.6$
${\Gamma }_{W\to e{\overline{\nu }}_e}$ [MeV]	$223\pm 6$
${\Gamma }_{W\to \mu {\overline{\nu }}_{\mu }}$ [MeV]	$222\pm 5$
${\Gamma }_{W\to \tau {\overline{\nu }}_{\tau }}$ [MeV]	$237\pm 6$

Table 2. Table 2 : Experimental upper bounds on LFV processes, along with the experiments that provided that bound. When more than one experiment is cited, the HFLAV average is used  Amhis:2016xyh . All upper bounds are given at the 90 % percent 90 90\% C.L..

Process	Branch. Frac.	Reference
LFV decay
${\mu }^{-}\to e^{-}\gamma$	$4.2\times {10}^{-13}$	MEG TheMEG:2016wtm
${\tau }^{-}\to e^{-}\gamma$	$5.4\times {10}^{-8}$	BaBar Aubert:2009ag ,Belle Hayasaka:2007vc
${\tau }^{-}\to {\mu }^{-}\gamma$	$5.0\times {10}^{-8}$	BaBar Aubert:2009ag ,Belle Hayasaka:2007vc
${\mu }^{-}\to e^{-}{e}^{-}{e}^{+}$	$1.0\times {10}^{-12}$	SINDRUM Bellgardt:1987du
${\tau }^{-}\to e^{-}{e}^{-}{e}^{+}$	$1.4\times {10}^{-8}$	BaBar Lees:2010ez ,Belle Hayasaka:2010np
${\tau }^{-}\to {\mu }^{-}{\mu }^{-}{\mu }^{+}$	$1.2\times {10}^{-8}$	ATLAS Aad:2016wce ,BaBar Lees:2010ez
		Belle Hayasaka:2010np ,LHCb Aaij:2014azz
${\tau }^{-}\to {\mu }^{-}{e}^{-}{e}^{+}$	$1.1\times {10}^{-8}$	BaBar Lees:2010ez ,Belle Hayasaka:2010np
${\tau }^{-}\to e^{-}{e}^{-}{\mu }^{+}$	$0.84\times {10}^{-8}$	BaBar Lees:2010ez ,Belle Hayasaka:2010np
${\tau }^{-}\to e^{-}{\mu }^{-}{\mu }^{+}$	$1.6\times {10}^{-8}$	BaBar Lees:2010ez ,Belle Hayasaka:2010np
${\tau }^{-}\to {\mu }^{-}{\mu }^{-}{e}^{+}$	$0.98\times {10}^{-8}$	BaBar Lees:2010ez ,Belle Hayasaka:2010np
LFV conversion
$\mu -e$ (Ti)	$4.3\times {10}^{-12}$	SINDRUM II Dohmen:1993mp
$\mu -e$ (Au)	$7\times {10}^{-13}$	SINDRUM II Bertl:2006up
$\mu -e$ (Pb)	$4.6\times {10}^{-11}$	SINDRUM II Honecker:1996zf

Table 3. Table 3 : Experimental measurements for all tests of lepton universality.

Obs.	Measured
$R_{e\mu }^{\pi }$	$\left(1.2327\pm 0.0023\right)\times {10}^{-4}$  Tanabashi:2018oca
$R_{e\mu }^K$	$\left(2.488\pm 0.010\right)\times {10}^{-5}$  Lazzeroni:2012cx
$R_{\mu e}^{\tau }$	$0.9762\pm 0.0028$  Amhis:2016xyh
$R_{\mu e}^W$	$0.980\pm 0.018$  Aaij:2016qqz
$R_{\tau e}^W$	$1.063\pm 0.027$  Schael:2013ita
$R_{\tau \mu }^W$	$1.070\pm 0.026$  Schael:2013ita
$R_K^B$	$0.745\pm 0.089$  Aaij:2014ora
$R_{K^{\ast }}^B$ (1)	$0.66\pm 0.09$  Aaij:2017vbb
$R_{K^{\ast }}^B$ (2)	$0.69\pm 0.10$  Aaij:2017vbb

Table 4. Table 4 : Experimental values of ( V C ​ K ​ M ) u ​ s subscript subscript 𝑉 𝐶 𝐾 𝑀 𝑢 𝑠 (V_{CKM})_{us} and the average value of ( V C ​ K ​ M ) u ​ d subscript subscript 𝑉 𝐶 𝐾 𝑀 𝑢 𝑑 (V_{CKM})_{ud} used in the calculation of the CKM likelihood. The factor f + ​ ( 0 ) = 0.959 ± 0.005 subscript 𝑓 0 plus-or-minus 0.959 0.005 f_{+}(0)=0.959\pm 0.005 is taken from Aoki:2016frl .

Parameter	Process	Value	Ref.
	$K_L\to \pi e\nu$	$0.2163\left(6\right)$
	$K_L\to \pi \mu \nu$	$0.2166\left(6\right)$
$\left|{\left(V_{\text{CKM}}^{\text{exp}}\right)}_{us}\right|f_{+}\left(0\right)$	$K_S\to \pi e\nu$	$0.2155\left(13\right)$	Antonelli:2010yf ; Aoki:2016frl
	$K^{\pm }\to {\pi }^{0}e\nu$	$0.2160\left(11\right)$
	$K^{\pm }\to {\pi }^0\mu \nu$	$0.2158\left(14\right)$
	$\frac{\text{BR}\left(\tau \to K\nu \right)}{\text{BR}\left(\tau \to \pi \nu \right)}$	$0.2262\left(13\right)$
$\left|{\left(V_{\text{CKM}}^{\text{exp}}\right)}_{us}\right|$	$\tau \to K\nu$	$0.2214\left(22\right)$	Follana:2007uv ; Amhis:2012bh
	$\tau \to l,\tau \to s$	$0.2173\left(22\right)$
$\left|{\left(V_{\text{CKM}}^{\text{exp}}\right)}_{ud}\right|$	Average	$0.97417\left(21\right)$	Patrignani:2016xqp

Table 5. Table 5 : Parameter ranges adopted for the full model scans, with + + ( − - ) for normal (inverted) hierarchy of the active neutrino masses.

Parameter	Value/Range	Prior
Active neutrino parameters
${\theta }_{12}$ [rad]	$[0.547684,0.628144]$	flat
${\theta }_{23}$ [rad]	$[0.670206,0.925025]$	flat
${\theta }_{13}$ [rad]	$[0.139452,0.155509]$	flat
$m_{{\nu }_0}$ [eV]	$[{10}^{-7},0.23]$	log
$\Delta m_{21}^2$ $\left[{10}^{-5}{\text{eV}}^2\right]$	$[6,9]$	flat
$\Delta m_{3l}^2$ $\left[{10}^{-3}{\text{eV}}^2\right]$	$[\pm 2,\pm 3]$	flat
${\alpha }_1$, ${\alpha }_2$ [rad]	$[0,2\pi ]$	flat
Sterile neutrino parameters
$\delta$ [rad]	$[0,2\pi ]$	flat
Re ${\omega }_{ij}$ [rad]	$[0,2\pi ]$	flat
Im ${\omega }_{ij}$	$[-15,15]$	flat
$M_I$ [GeV]	$[0.06,500]$	log
$R_{\mathrm{order}}$	[1,6]	flat
Nuisance parameters
$m_H$ [GeV]	$[124.1,127.3]$	flat

Table 6. Table 6 : Set of targeted scans performed for normal (N) and inverted (I) hierarchy in addition to the full parameter scans. Parameters not shown in this table are taken as in Table  5 .

$M_1$ [GeV]	$\Delta M_{21}$ [GeV]	$m_{{\nu }_0}$ [eV]	$\alpha$	$(s,m)$	Hierarchy	Other
$[0.1,0.3162]$	$[{10}^{-10},0.1]$	$[{10}^{-7},0.23]$	$(e,\mu ,\tau )$	$(0.5,-0.5)$	N, I
$[0.1,0.4217]$	$[{10}^{-10},0.1]$	$[{10}^{-7},0.23]$	$(e,\mu ,\tau )$	$(0.5,-0.5)$	N, I
$[0.3162,2.0]$	$[{10}^{-10},0.1]$	$[{10}^{-7},0.23]$	$(e,\mu ,\tau )$	$(0.5,0.5)$	N, I
$[2.0,60]$	$[{10}^{-20},{10}^{-10}]$	$[{10}^{-6},0.23]$	$(e,\mu ,\tau )$	$(0.5,0)$	N, I	, flat $m_{{\nu }_0}$ prior
$[2.0,60]$	$[{10}^{-20},{10}^{-10}]$	${10}^{-4},{10}^{-5},{10}^{-6}$	$(e,\mu ,\tau )$	$(0.5,0)$	N, I	fixed $m_{{\nu }_0}$
$[60,500]$	$[{10}^{-20},{10}^{-10}]$	$[{10}^{-6},0.23]$	$(e,\mu ,\tau )$	$(0.7,0.25)$	N, I	, flat $m_{{\nu }_0}$ prior
$[60,500]$	$[{10}^{-20},{10}^{-10}]$	${10}^{-4},{10}^{-5},{10}^{-6}$	$(e,\mu ,\tau )$	$(0.7,0.25)$	N, I	fixed $m_{{\nu }_0}$
$[0.06,0.14]$	$[{10}^{-10},0.1]$	$[{10}^{-7},0.23]$	$(e,\mu )$	$(0.5,-0.5)$	N, I	, flat $m_{{\nu }_0}$ prior
$[60,500]$	$[{10}^{-20},{10}^{-10}]$	$[{10}^{-6},0.23]$	$(e,\mu )$	$(0.7,0.25)$	I	flat $m_{{\nu }_0}$ prior
$[0.14,0.2]$	$[{10}^{-10},0.1]$	$[{10}^{-7},0.23]$	$\left(e\right)$	$(0.5,-0.5)$	N,I	, flat $m_{{\nu }_0}$ prior
$[0.2,0.4217]$	$[{10}^{-10},0.1]$	$[{10}^{-7},0.23]$	$\left(e\right)$	$(0.5,-0.5)$	N,I	, flat $m_{{\nu }_0}$ prior
$[0.14,0.3162]$	$[{10}^{-10},0.1]$	$[{10}^{-7},0.23]$	$\left(\mu \right)$	$(0.5,-0.5)$	N,I	, flat $m_{{\nu }_0}$ prior
$[0.1,0.3162]$	$[{10}^{-10},0.1]$	$[{10}^{-7},0.23]$	$\left(\tau \right)$	$(0.5,-0.5)$	N, I	-
$[0.1,0.4217]$	$[{10}^{-10},0.1]$	$[{10}^{-7},0.23]$	$\left(\tau \right)$	$(0.5,-0.5)$	N, I	-
$[0.175,0.3611]$	$[{10}^{-20},{10}^{-10}]$	$[{10}^{-2},0.23]$	$\left(\tau \right)$	$(0.5,0.5)$	N, I	-
$[0.25,0.3611]$	$[{10}^{-20},{10}^{-10}]$	$[{10}^{-2},0.23]$	$\left(\tau \right)$	$(0.5,0.5)$	N, I	-
$[0.25,0.4]$	$[{10}^{-20},{10}^{-10}]$	$[{10}^{-2},0.23]$	$\left(\tau \right)$	$(1.0,0)$	N, I	-
$[0.3611,0.4492]$	$[{10}^{-10},0.1]$	$[{10}^{-2},0.23]$	$\left(\tau \right)$	$(0.5,-0.5)$	N, I	-
$[0.3611,0.4492]$	$[{10}^{-20},{10}^{-10}]$	$[{10}^{-2},0.23]$	$\left(\tau \right)$	$(0.5,0.5)$	N, I	-
$[0.4,0.5]$	$[{10}^{-10},0.1]$	$[{10}^{-7},0.23]$	$\left(\tau \right)$	$(0.5,-0.5)$	N, I	-
$[0.3162,2.0]$	$[{10}^{-10},0.1]$	$[{10}^{-7},0.23]$	$\left(\tau \right)$	$(0.5,0.5)$	N, I	-
$[0.3162,1.4]$	$[{10}^{-10},0.1]$	$[0.03,0.23]$	$\left(\tau \right)$	$(0.5,0.5)$	I
$[1.0,1.5]$	$[{10}^{-7},0.01]$	$[0.03,0.23]$	$\left(\tau \right)$	$(0.5,0.5)$	I
$[1.25,1.45]$	$[{10}^{-20},{10}^{-10}]$	$[0.01,0.23]$	$\left(\tau \right)$	$(0.5,0.5)$	I	-
$[1.4,1.78]$	$[{10}^{-7},0.01]$	$[0.03,0.23]$	$\left(\tau \right)$	$(0.5,0.5)$	I
$[1.65,1.85]$	$[{10}^{-20},{10}^{-10}]$	$[0.01,0.23]$	$\left(\tau \right)$	$(0.5,0.5)$	I	-
$[1.25,1.45]$	$[{10}^{-20},{10}^{-10}]$	$[0.01,0.23]$	$\left(\tau \right)$	$(0.5,0.5)$	I	-
$[2.0,60]$	$[{10}^{-20},{10}^{-10}]$	$[{10}^{-6},0.23]$	$\left(\tau \right)$	$(0.5,0)$	N, I	flat $m_{{\nu }_0}$ prior
$[60,500]$	$[{10}^{-20},{10}^{-10}]$	$[{10}^{-6},0.23]$	$\left(\tau \right)$	$(0.7,0.25)$	N, I	flat $m_{{\nu }_0}$ prior

Table 7. Table 7 : Upper limits on U α ​ I 2 superscript subscript 𝑈 𝛼 𝐼 2 U_{\alpha I}^{2} and U α ​ I ​ U β ​ I subscript 𝑈 𝛼 𝐼 subscript 𝑈 𝛽 𝐼 U_{\alpha I}U_{\beta I} within 2 σ 𝜎 \sigma in the high mass region M I ≳ 80 greater-than-or-equivalent-to subscript 𝑀 𝐼 80 M_{I}\gtrsim 80 GeV, for normal (N) and inverted (I) hierarchy.

Hierarchy	Coupling	Upper limit ($\mathrm{𝟐}\sigma$)
N	$U_{eI}^2$	$4.92\times {10}^{-4}$
N	$U_{\mu I}^2$	$2.42\times {10}^{-4}$
N	$U_{\tau I}^2$	$9.59\times {10}^{-3}$
N	$U_{eI}{U}_{\mu I}$	$3.49\times {10}^{-5}$
N	$U_{eI}{U}_{\tau I}$	$1.37\times {10}^{-3}$
N	$U_{\mu I}{U}_{\tau I}$	$1.25\times {10}^{-3}$
I	$U_{eI}^2$	$8.15\times {10}^{-4}$
I	$U_{\mu I}^2$	$3.46\times {10}^{-4}$
I	$U_{\tau I}^2$	$9.91\times {10}^{-3}$
I	$U_{eI}{U}_{\mu I}$	$3.34\times {10}^{-5}$
I	$U_{eI}{U}_{\tau I}$	$2.19\times {10}^{-3}$
I	$U_{\mu I}{U}_{\tau I}$	$1.43\times {10}^{-3}$

Table 8. Table 8 : Capabilities and module functions implemented for active and sterile neutrino masses and mixings.

Capability	Function (Return Type): Brief description	Dependencies	Options
ordering	ordering(bool): Calculates the hierarchy type.	None
m_nu	M_nu(Eigen::Matrix3cd): Calculates the diagonalised LHN mass matrix.	ordering	ordering(bool)
md21	md21(double): Calculates the square mass splitting $\Delta m_{21}^2$.	m_nu
md31	md31(double): Calculates the square mass splitting $\Delta m_{31}^2$.	m_nu
md32	md32(double): Calculates the square mass splitting $\Delta m_{32}^2$.	m_nu
min_mass	min_mass(double): Calculates the minimal neutrino mass.	ordering, m_nu
UPMNS	UPMNS(Eigen::Matrix3cd): Calculates the PMNS matrix.	None
SeesawI_Theta	CI_Theta(Eigen::Matrix3cd): Calculates the active-sterile mixing matrix in seesaw type-I using the C-I parametrisation.	m_nu, UPMNS, SMINPUTS
SeesawI_Vnu	Vnu(Eigen::Matrix3cd): Calculates the active mixing matrix in seesaw type-I.	UPMNS, SeesawI_Theta
Unitarity	Unitarity_UPMNS(bool): Checks for unitarity in the SM neutrino mixing matrix.	m_nu, UPMNS
	Unitarity_SeesawI(bool): Checks for unitarity in the full neutrino mixing matrix in seesaw type-I.	m_nu, SeesawI_Theta, SeesawI_Vnu

Table 9. Table 9 : Capabilities and module functions implemented that calculate log-likelihoods for the active neutrino parameters.

Capability	Function (Return Type): Brief description	Dependencies
md21_lnL	md21_lnL(double): Computes the log-likehood for $\Delta m_{21}^2$.	ordering, md21
md3l_lnL	md3l_lnL(double): Computes the log-likehood for $\Delta m_{31}^2$ for normal hierarchy or $\Delta m_{32}^2$ for inverted.	ordering, md31, md32
deltaCP_lnL	deltaCP_lnL(double): Computes the log-likehood for ${\delta }_{CP}$.	ordering, deltaCP
theta12_lnL	theta12_lnL(double): Computes the log-likehood for ${\theta }_{12}$.	ordering, theta12
theta23_lnL	theta23_lnL(double): Computes the log-likehood for ${\theta }_{23}$.	ordering, theta23
theta13_lnL	theta13_lnL(double): Computes the log-likehood for ${\theta }_{13}$.	ordering, theta13
sum_mnu_lnL	sum_mnu_lnL(double): Computes the log-likehood for $\sum m_{\nu }$	None

Table 10. Table 10 : Capabilities and module functions that calculate electroweak precision observables and their likelihoods.

Capability	Function (Return Type): Brief description	Dependencies
prec_sinW2_eff	RHN_sinW2_eff(triplet<double>): Calculates $s_{eff}^2$.	SeesawI_Theta
mW	RHN_mW(triplet<double>): Calculates $m_W$.	sinW2, SeesawI_Theta
Z_gamma_nu	Z_gamma_nu_2l(triplet<double>): Calculates the decay width of $Z$ to neutrinos.	SM_spectrum, SeesawI_Theta, SeesawI_Vnu
W_to_l_decays	RHN_W_to_l_decays(vector<double>): Calculates the decay width of the processes $W\to l\nu$.	SMINPUTS, mw, SeesawI_Theta
lnL_sinW2_eff	lnL_sinW2_eff_chi2(double): Computes the log-likehood for $s_{eff}^2$.	prec_sinW2_eff
lnL_W_mass	lnL_W_mass_chi2(double): Computes the log-likehood for $m_W$ .	mW
lnL_Z_inv	lnL_Z_inv(double): Computes the log-likehood for ${\Gamma }_{inv}$.	Z_gamma_nu
lnL_W_decays	lnL_W_decays_chi2(double): Computes the log-likehood for ${\Gamma }_{W\to l\nu }$.	W_to_l_decays, W_plus_decay_rates

Table 11. Table 11 : Capabilities and module functions to calculate LFV observables.

Capability	Function (Return Type): Brief description	Dependencies
muegamma	RHN_muegamma(double): Calculates $BR\left({\mu }^{-}\to e^{-}\gamma \right)$.	SMINPUTS, m_nu, SeesawI_Vnu, SeesawI_Theta, mu_minus_decay_rates
tauegamma	RHN_tauegamma(double): Calculates $BR\left({\tau }^{-}\to e^{-}\gamma \right)$.	SMINPUTS, m_nu, SeesawI_Vnu, SeesawI_Thet, tau_minus_decay_rates
taumugamma	RHN_taumugamma(double): Calculates $BR\left({\tau }^{-}\to {\mu }^{-}\gamma \right)$.	SMINPUTS, m_nu, SeesawI_Vnu, SeesawI_Theta, tau_minus_decay_rates
mueee	RHN_mueee(double): Calculates $BR\left({\mu }^{-}\to e^{-}{e}^{-}{e}^{+}\right)$.	SMINPUTS, m_nu, SeesawI_Vnu, SeesawI_Theta, mu_minus_decay_rates
taueee	RHN_taueee(double): Calculates $BR\left({\tau }^{-}\to e^{-}{e}^{-}{e}^{+}\right)$.	SMINPUTS, m_nu, SeesawI_Vnu, SeesawI_Theta, tau_minus_decay_rates
taumumumu	RHN_taumumumu(double): Calculates $BR\left({\tau }^{-}\to {\mu }^{-}{\mu }^{-}{\mu }^{+}\right)$.	SMINPUTS, m_nu, SeesawI_Vnu, SeesawI_Theta, tau_minus_decay_rates
taumuee	RHN_taumuee(double): Calculates $BR\left({\tau }^{-}\to {\mu }^{-}{e}^{-}{e}^{+}\right)$.	SMINPUTS, m_nu, SeesawI_Vnu, SeesawI_Theta, tau_minus_decay_rates
taueemu	RHN_taueemu(double): Calculates $BR\left({\tau }^{-}\to e^{-}{e}^{-}{\mu }^{+}\right)$.	SMINPUTS, m_nu, SeesawI_Vnu, SeesawI_Theta, tau_minus_decay_rates
tauemumu	RHN_tauemumu(double): Calculates $BR\left({\tau }^{-}\to e^{-}{\mu }^{-}{\mu }^{+}\right)$.	SMINPUTS, m_nu, SeesawI_Vnu, SeesawI_Theta, tau_minus_decay_rates
taumumue	RHN_taumumue(double): Calculates $BR\left({\tau }^{-}\to {\mu }^{-}{\mu }^{-}{e}^{+}\right)$.	SMINPUTS, m_nu, SeesawI_Vnu, SeesawI_Theta, tau_minus_decay_rates
mueTi	RHN_mueTi(double): Calculates $R\left(\mu -e\right)$ in a Ti nucleus.	SMINPUTS, m_nu, SeesawI_Vnu, SeesawI_Theta
mueAu	RHN_mueAu(double): Calculates $R\left(\mu -e\right)$ in a Au nucleus.	SMINPUTS, m_nu, SeesawI_Vnu, SeesawI_Theta
muePb	RHN_muePb(double): Calculates $R\left(\mu -e\right)$ in a Pb nucleus.	SMINPUTS, m_nu, SeesawI_Vnu, SeesawI_Theta

Table 12. Table 12 : Capabilities and module functions for the likelihoods computed for the LFV observables.

Capability	Function (Return Type): Brief description	Dependencies
l2lgamma_lnL	l2lgamma_likelihood(double): Computes the log-likelihood for $l^{-}\to l^{-}\gamma$.	muegamma, tauegamma, taumugamma
l2lll_lnL	l2lll_likelihood(double): Computes the log-likelihood for $l^{-}\to l^{-}{l}^{-}{l}^{+}$.	mueee, taueee, taumumumu, taumumue, tauemumu
mu2e_lnL	mu2e_likelihood(double): Computes the log-likelihood associated with $\mu -e$ conversion.	mueTi, muAu, muePb

Table 13. Table 13 : Capabilities and module functions that calculate lepton universality observables.

Capability	Function (Return Type): Brief description	Dependencies
R_pi	RHN_R_pi(double): Calculates the test of lepton universality $R_{e\mu }^{\pi }$.	SMINPUTS, SeesawI_Theta, SeesawI_Vnu
R_K	RHN_R_K(double): Calculates the test of lepton universality $R_{e\mu }^K$..	SMINPUTS, SeesawI_Theta, SeesawI_Vnu
R_tau	RHN_R_tau(double): Calculates the test of lepton universality $R_{\mu e}^{\tau }$.	SMINPUTS, SeesawI_Theta
R_W	RHN_R_W(vector<double>): Calculates the test of lepton universality $R_{\alpha \beta }^W$.	W_to_l_decays
RK	RHN_RK(double): Calculates the test of lepton universality $R_K$.	SMINPUTS, SeesawI_Theta
RKstar_0045_11	RHN_RKstar_0045_11(double): Calculates the test of lepton universality $R_{K^{\ast }}$ for the range .	SMINPUTS, SeesawI_Theta
RKstar_11_60	RHN_RKstar_11_60(double): Calculates the test of lepton universality $R_{K^{\ast }}$ for the range .	SMINPUTS, SeesawI_Theta

Table 14. Table 14 : Capabilities and module functions for the likelihoods computed for lepton universality tests.

Capability	Function (Return Type): Brief description	Dependencies
lnL_R_pi	lnL_R_pi(double): Calculates the total log-likelihood for lepton universality tests on leptonic decays of $\pi$ mesons.	R_pi
lnL_R_K	lnL_R_K(double): Calculates the total log-likelihood for lepton universality tests on leptonic decays of $K$ mesons.	R_K
lnL_R_tau	lnL_R_tau(double): Calculates the total log-likelihood for lepton universality tests on leptonic decays of $\tau$ leptons.	R_tau
lnL_R_W	lnL_R_W(double): Calculates the total log-likelihood for lepton universality tests on leptonic decays of $W$ bosons.	R_W
LUV_M	LUV_measurements(): Collates all measurements of semi-leptonic tests of lepton universality in $B$ meson decays.	RK, RKstar_0045_11, RKstar_11_60
LUV_LL	lnL_lepuniv(double): Calculates the total log-likelihood for semi-leptonic tests of lepton universality in $B$ meson decays.	LUV_M

Table 15. Table 15 : Capabilities and module functions implemented to calculate CKM unitarity and its likelihood.

Capability	Function (Return Type): Brief description	Dependencies
calc_Vus	calc_Vus(double): Calculates the profiling value of $V_{us}$ for a particular $\Theta$.	SMINPUTS, SeesawI_Theta
lnLckm_Vus	lnL_ckm_Vus(double): Computes the total log-likelihood from CKM unitarity for a given parameter $V_{us}$	SMINPUTS, SeesawI_Theta
lnLckm_Vusmin	lnL_ckm_Vusmin(double): Computes the total log-likelihood from CKM unitarity profiling over $V_{us}$.	SMINPUTS, SeesawI_Theta, calc_Vus

Table 16. Table 16 : Capabilities and module functions implemented to calculate neutrinoless double-beta decay observables and likelihood.

Capability	Function (Return Type): Brief description	Dependencies
Thalf_0nubb_Xe	RHN_Thalf_0nubb_Xe(double): Calculates the half-life for Xe detector.	m_nu, UPMNS, SeesawI_Theta
Thalf_0nubb_Ge	RHN_Thalf_0nubb_Ge(double): Calculates the half-life for Ge detector.	m_nu, UPMNS, SeesawI_Theta
mbb_0nubb_Xe	RHN_mbb_0nubb_Xe(double): Calculates $m_{\beta \beta }$ for Xe detector.	m_nu, UPMNS, SeesawI_Theta
mbb_0nubb_Ge	RHN_mbb_0nubb_Ge(double): Calculates $m_{\beta \beta }$ for Ge detector.	m_nu, UPMNS, SeesawI_Theta
lnL_0nubb_KamLAND_Zen	lnL_0nubb_KamLAND_Zen(double): Calculates KamLAND-Zen log-likelihood based on half-life.	Thalf_0nubb_Xe
lnL_0nubb_GERDA	lnL_0nubb_GERDA(double): Calculates GERDA log-likelihood based on half-life.	Thalf_0nubb_Ge
lnL_0nubb	lnL_0nubb(double): Calculates the total log-likelihood based on half-life.	lnL_0nubb_KamLAND_Zen, lnL_0nubb_GERDA
lnL_mbb_0nubb_KamLAND_Zen	lnL_mbb_0nubb_KamLAND_Zen(double): Calculates KamLAND-Zen log-likelihood based on $m_{\beta \beta }$.	mbb_0nubb_Xe
lnL_mbb_0nubb_GERDA	lnL_mbb_0nubb_GERDA(double): Calculates GERDA log-likelihood based on $m_{\beta \beta }$.	mbb_0nubb_Ge
lnL_mbb_0nubb	lnL_mbb_0nubb(double): Calculates the total log-likelihood based on $m_{\beta \beta }$.	lnL_mbb_0nubb_KamLAND_Zen, lnL_mbb_0nubb_GERDA

Table 17. Table 17 : Capabilities and module functions implemented in NeutrinoBit to calculate BBN observables and likelihood for sterile neutrino models.

Capability	Function (Return Type): Brief description	Dependencies
Gamma_RHN2piplusl	Gamma_RHN2piplusl(std::vector<double>): Calculates $\Gamma \left(N_I\to {\pi }^{+}{l}_{\alpha }^{-}\right)$.	SMINPUTS, SeesawI_Theta
Gamma_RHN2Kplusl	Gamma_RHN2Kplusl(std::vector<double>): Calculates $\Gamma \left(N_I\to K^{+}{l}_{\alpha }^{-}\right)$.	SMINPUTS, SeesawI_Theta
Gamma_RHN2Dplusl	Gamma_RHN2Dplusl(std::vector<double>): Calculates $\Gamma \left(N_I\to D^{+}{l}_{\alpha }^{-}\right)$.	SMINPUTS, SeesawI_Theta
Gamma_RHN2Dsl	Gamma_RHN2Dsl(std::vector<double>): Calculates $\Gamma \left(N_I\to D_s{l}_{\alpha }^{-}\right)$.	SMINPUTS, SeesawI_Theta
Gamma_RHN2Bplusl	Gamma_RHN2Bplusl(std::vector<double>): Calculates $\Gamma \left(N_I\to B^{+}{l}_{\alpha }^{-}\right)$.	SMINPUTS, SeesawI_Theta
Gamma_RHN2Bcl	Gamma_RHN2Bcl(std::vector<double>): Calculates $\Gamma \left(N_I\to B_c{l}_{\alpha }^{-}\right)$.	SMINPUTS, SeesawI_Theta
Gamma_RHN2pi0nu	Gamma_RHN2pi0nu(std::vector<double>): Calculates $\Gamma \left(N_I\to {\pi }^0{\nu }_{\alpha }\right)$.	SMINPUTS, SeesawI_Theta
Gamma_RHN2etanu	Gamma_RHN2etanu(std::vector<double>): Calculates $\Gamma \left(N_I\to \eta {\nu }_{\alpha }\right)$.	SMINPUTS, SeesawI_Theta
Gamma_RHN2etaprimenu	Gamma_RHN2etaprimenu(std::vector<double>): Calculates $\Gamma \left(N_I\to {\eta }^{\prime }{\nu }_{\alpha }\right)$.	SMINPUTS, SeesawI_Theta
Gamma_RHN2etacnu	Gamma_RHN2etacnu(std::vector<double>): Calculates $\Gamma \left(N_I\to {\eta }_c{\nu }_{\alpha }\right)$.	SMINPUTS, SeesawI_Theta
Gamma_RHN2rhoplusl	Gamma_RHN2rhoplusl(std::vector<double>): Calculates $\Gamma \left(N_I\to {\rho }^{+}{l}_{\alpha }^{-}\right)$.	SMINPUTS, SeesawI_Theta
Gamma_RHN2Dstarplusl	Gamma_RHN2Dstarplusl(std::vector<double>): Calculates $\Gamma \left(N_I\to D^{\ast +}{l}_{\alpha }^{-}\right)$.	SMINPUTS, SeesawI_Theta
Gamma_RHN2Dstarsl	Gamma_RHN2Dstarsl(std::vector<double>): Calculates $\Gamma \left(N_I\to D_s^{\ast }{l}_{\alpha }^{-}\right)$.	SMINPUTS, SeesawI_Theta
Gamma_RHN2rho0nu	Gamma_RHN2rho0nu(std::vector<double>): Calculates $\Gamma \left(N_I\to {\rho }^0{\nu }_{\alpha }\right)$.	SMINPUTS, SeesawI_Theta
Gamma_RHN2omeganu	Gamma_RHN2omeganu(std::vector<double>): Calculates $\Gamma \left(N_I\to \omega {\nu }_{\alpha }\right)$.	SMINPUTS, SeesawI_Theta
Gamma_RHN2phinu	Gamma_RHN2phinu(std::vector<double>): Calculates $\Gamma \left(N_I\to \varphi {\nu }_{\alpha }\right)$.	SMINPUTS, SeesawI_Theta
Gamma_RHN2Jpsinu	Gamma_RHN2Jpsinu(std::vector<double>): Calculates $\Gamma \left(N_I\to J/\psi {\nu }_{\alpha }\right)$.	SMINPUTS, SeesawI_Theta
Gamma_RHN23nu	Gamma_RHN23nu(std::vector<double>): Calculates $\Gamma \left(N_I\to {\nu }_{\alpha }\overline{{\nu }_{\beta }}\nu \beta \right)$.	SMINPUTS, SeesawI_Theta
Gamma_RHN2llnu	Gamma_RHN2llnu(std::vector<double>): Calculates $\Gamma \left(N_I\to l_{\alpha \ne \beta }^{-}{l}_{\beta }^{+}{\nu }_{\beta }\right)$.	SMINPUTS, SeesawI_Theta
Gamma_RHN2null	Gamma_RHN2Kplusl(std::vector<double>): Calculates $\Gamma \left(N_I\to {\nu }_{\alpha }{l}_{\beta }^{+}{l}_{\beta }^{-}\right)$.	SMINPUTS, SeesawI_Theta
Gamma_RHN2nuuubar	Gamma_RHN2nuuubar(std::vector<double>): Calculates $\Gamma (N_I\to {\nu }_{\alpha }{q}^u\overline{q^u}$.	SMINPUTS, SeesawI_Theta
Gamma_RHN2nuddbar	Gamma_RHN2nuddbar(std::vector<double>): Calculates $\Gamma (N_I\to {\nu }_{\alpha }{q}^d\overline{q^d}$.	SMINPUTS, SeesawI_Theta
Gamma_RHN2ludbar	Gamma_RHN2ludbar(std::vector<double>): Calculates $\Gamma (N_I\to l_{\alpha }{q}_{\beta }^u\overline{q_{\gamma }^d}$.	SMINPUTS, SeesawI_Theta
Gamma_BBN	Gamma_BBN(std::vector<double>): Calculates the total decay width for each RHN.	Gamma_*, as listed above
lnL_bbn	lnL_bbn(double): Calculates the log-likelihood.	Gamma_BBN

Table 18. Table 18 : Capabilities and module functions that calculate the magnitude and argument of the matrix elements of Θ Θ \Theta ( I = 1,2,3 ).

Capability	Function (Return Type): Brief description	Dependencies	Options
UeI	UeI(double): Magnitude of the matrix element $U_{eI}=\left|{\Theta }_{eI}\right|$.	SeesawI_Theta	upper_limit(double), lower_limit(double)
UmuI	UmuI(double): Magnitude of the matrix element $U_{\mu I}=\left|{\Theta }_{\mu I}\right|$.	SeesawI_Theta	upper_limit(double), lower_limit(double)
UtauI	UtauI(double): Magnitude of the matrix element $U_{\tau I}=\left|{\Theta }_{\tau I}\right|$.	SeesawI_Theta	upper_limit(double), lower_limit(double)
UeI_phase	UeI_phase(double): Argument of the matrix element ${\Theta }_{eI}$.	SeesawI_Theta
UmuI_phase	UmuI_phase(double): Argument of the matrix element ${\Theta }_{\mu I}$.	SeesawI_Theta
UtauI_phase	UtauI_phase(double): Argument of the matrix element ${\Theta }_{\tau I}$.	SeesawI_Theta

Table 19. Table 19 : Capabilities and module functions implemented in NeutrinoBit to calculate direct search likelihoods for sterile neutrino models.

Capability	Function (Return Type): Brief description	Dependencies
lnLpienu	lnL_pienu(double): Calculates the log-likelihood for PIENU.	Ue1, Ue2, Ue3
lnLps191e	lnL_ps191_e(double): Calculates the log-likelihood for PS-191 in the electron sector.	Ue1, Ue2, Ue3, Um1, Um2, Um3, Ut1, Ut2, Ut3
lnLps191mu	lnL_ps191_mu(double): Calculates the log-likelihood for PS-191 in the muon sector.	Ue1, Ue2, Ue3, Um1, Um2, Um3, Ut1, Ut2, Ut3
lnLcharme	lnL_charm_e(double): Calculates the log-likelihood for CHARM in the electron sector.	Ue1, Ue2, Ue3, Um1, Um2, Um3, Ut1, Ut2, Ut3
lnLcharmmu	lnL_charm_mu(double): Calculates the log-likelihood for CHARM in the muon sector.	Ue1, Ue2, Ue3, Um1, Um2, Um3, Ut1, Ut2, Ut3
lnLcharmtau	lnL_charm_tau(double): Calculates the log-likelihood for CHARM in the tau sector.	Ut1, Ut2, Ut3
lnLdelphi_shortlived	lnL_delphi_short_lived(double): Calculates the log-likelihood for DELPHI’s short-lived RHN analyses.	Ue1, Ue2, Ue3, Um1, Um2, Um3, Ut1, Ut2, Ut3
lnLdelphi_longlived	lnL_delphi_long_lived(double): Calculates the log-likelihood for DELPHI’s long-lived RHN analyses.	Ue1, Ue2, Ue3, Um1, Um2, Um3, Ut1, Ut2, Ut3
lnLatlase	lnL_atlas_e(double): Calculates the log-likelihood for ATLAS in the electron sector.	Ue1, Ue2, Ue3
lnLatlasmu	lnL_atlas_mu(double): Calculates the log-likelihood for ATLAS in the muon sector.	Um1, Um2, Um3
lnLlhce	lnL_lhc_e(double): Calculates the log-likelihood for CMS in the electron sector.	Ue1, Ue2, Ue3
lnLlhcmu	lnL_lhc_mu(double): Calculates the log-likelihood for CMS in the muon sector.	Um1, Um2, Um3
lnLe949	lnL_e949(double): Calculates the log-likelihood for E949.	Um1, Um2, Um3
lnLnutev	lnL_nutev(double): Calculates the log-likelihood for NuTeV.	Um1, Um2, Um3

Table 20. Table 20 : Capabilities and module functions for perturbativity constraints and coupling slide.

Capability	Function (Return Type): Brief description	Dependencies
perturbativity_lnL	perturbativity_likelihood(double): Calculates the log-likelihood for the perturbativity of Yukawa couplings.	SMINPUTS, SeesawI_Theta
RHN_coupling_slide	coupling_slide(double): Calculates the log-likelihood for the coupling slide.	SeesawI_Theta, Ut1, Ut2, Ut3

Table 21. Table 21 : Effective atomic number Z eff subscript 𝑍 eff Z_{\text{eff}} , nuclear form factor F p subscript 𝐹 𝑝 F_{p} and capture rate Γ capt subscript Γ capt \Gamma_{\text{capt}} for the relevant nuclei.

${𝐍}_{𝐙}^{𝐙+𝐍}$	${𝐙}_{\text{eff}}$	${𝐅}_{𝐩}$	${𝚪}_{\text{capt}}\left({10}^6{s}^{-1}\right)$
${\text{Ti}}_{22}^{48}$	17.6	0.54	2.59
${\text{Au}}_{79}^{197}$	33.5	0.16	13.07
${\text{Pb}}_{82}^{208}$	34.0	0.15	13.45

Table 22. Table 22 : Decay constants (in MeV) of pseudoscalar charged (top left), pseudoscalar neutral (top right), vector charged (bottom left) and vector neutral (bottom right) mesons.

$f_{{\pi }^{+}}$	130.2 Rosner:2015wva	$f_{{\pi }^0}$	130.2 Rosner:2015wva
$f_{K^{+}}$	155.7 Rosner:2015wva	$f_{\eta }$	81.7 Bondarenko:2018ptm
$f_{D^{+}}$	212.6 Rosner:2015wva	$f_{{\eta }^{\prime }}$	-94.7 Bondarenko:2018ptm
$f_{D_s}$	249.9 Rosner:2015wva	$f_{{\eta }_c}$	237 Bondarenko:2018ptm
$f_{B^{+}}$	190 Rosner:2015wva
$f_{B_c}$	434 Colquhoun:2015oha
$f_{{\rho }^{+}}$	209 Ebert:2006hj	$f_{{\rho }^0}$	209 Ebert:2006hj
$f_{D^{\ast +}}$	246.75 Dhiman:2017urn	$f_{\omega }$	195 Atre:2009rg
$f_{D_s^{\ast }}$	284 Dhiman:2017urn	$f_{\varphi }$	229 Ebert:2006hj
		$f_{J/\psi }$	418 Becirevic:2013bsa

Equations

$$\mathcal{L}$$

$$-\frac{1}{2}\left(\overline{{\nu }_R^c}{M}_M{\nu }_R+\overline{{\nu }_R}{M}_M^{†}{\nu }_R^c\right).$$

$$\frac{1}{2}(\overline{{\nu }_L}\overline{{\nu }_R^c})\mathcal{M}\left(\begin{array}{c}{\nu }_L^c\\ {\nu }_R\end{array}\right),$$

$$\mathcal{M}=\left(\begin{array}{cc}\delta m_{\nu }^{1loop}& M_D\\ M_D^T& M_M+\delta M_N^{1loop}\end{array}\right),$$

$$\mathcal{U}=\left(\begin{array}{cc}\cos (\theta )& \sin (\theta )\\ -\sin ({\theta }^{†})& \cos ({\theta }^{†})\end{array}\right)\left(\begin{array}{cc}{U}_{\nu }& \\ & U_N^{\ast }\end{array}\right)$$

$$\cos (\theta )$$

$$\sin (\theta )$$

$${\mathcal{U}}^{†}\mathcal{M}{\mathcal{U}}^{\ast }=\left(\begin{array}{cc}{m}_{\nu }^{\mathrm{diag}}& \\ & M_N^{\mathrm{diag}}\end{array}\right)$$

$$M_N^{\mathrm{diag}}$$

$$m_{\nu }^{\mathrm{diag}}$$

$$\theta =M_D{M}_M^{-1}=vF{M}_M^{-1}$$

$$\mathcal{U}=\Bigg{[}\begin{pmatrix}\mathbb{I}-\frac{1}{2}\theta\theta^{\dagger}&\theta\\ -\theta^{\dagger}&\mathbb{I}-\frac{1}{2}\theta^{\dagger}\theta\end{pmatrix}+\mathcal{O}(\theta^{3})\Bigg{]}\begin{pmatrix}U_{\nu}&\\ &U_{N}^{\ast}\end{pmatrix},$$

$$m_{\nu }=m_{\nu }^{\mathrm{tree}}+\delta m_{\nu }^{1loop}$$

$$m_{\nu }^{\mathrm{tree}}$$

$$M_N=M_M+\frac{1}{2}\left({\theta }^{†}\theta M_M+M_M^T{\theta }^T{\theta }^{\ast }\right)+\delta M_N^{1loop}.$$

$${\left(\delta m_{\nu }^{1loop}\right)}_{\alpha \beta }=\sum _I{F}_{\alpha I}{M}_I{F}_{I\beta }^{T}l(M_I),$$

$$\begin{gathered} l(M_I)=\frac{1}{(4\pi {)}^2}[\left(\frac{3\text{ln}[(M_I/m_Z{)}^2]}{(M_I/m_Z{)}^2-1}\right) \\ +\left(\frac{\text{ln}[(M_I/m_H{)}^2]}{(M_I/m_H{)}^2-1}\right)]. \end{gathered}$$

$$\nu =V_{\nu }^{†}{\nu }_L-U_{\nu }^{†}\theta {\nu }_R^c+V_{\nu }^T{\nu }_L^c-U_{\nu }^T{\theta }^{\ast }{\nu }_R$$

$$N=V_N^{†}{\nu }_R+{\Theta }^T{\nu }_L^c+V_N^T{\nu }_R^c+{\Theta }^{†}{\nu }_L,$$

$$V_{\nu }$$

$$V_N$$

$$\Theta =\theta U_N^{\ast }.$$

$$-\frac{g}{\sqrt{2}}{\overline{N}}_I{\Theta }_{I\alpha }^{†}{\gamma }^{\mu }{e}_{L\alpha }{W}_{\mu }^{+}-\frac{g}{\sqrt{2}}\overline{e_{L\alpha }}{\gamma }^{\mu }{\Theta }_{\alpha I}{N}_I{W}_{\mu }^{-}$$

$$-\frac{g}{2\cos {\theta }_W}\overline{N_I}{\Theta }_{I\alpha }^{†}{\gamma }^{\mu }{\nu }_{L\alpha }{Z}_{\mu }-\frac{g}{2\cos {\theta }_W}\overline{{\nu }_{L\alpha }}{\gamma }^{\mu }{\Theta }_{\alpha I}{N}_i{Z}_{\mu }$$

$$\displaystyle-\frac{g}{\sqrt{2}}\frac{M_{I}}{m_{W}}\Theta_{\alpha i}h\overline{\nu_{L\alpha}}N_{I}-\frac{g}{\sqrt{2}}\frac{M_{I}}{m_{W}}\Theta^{\dagger}_{I\alpha}h\overline{N_{I}}\nu_{L\alpha}\$$

$$U_{\alpha I}^2$$

$$U_I^2$$

$$U_{\alpha }^2$$

$$(\delta m_{\nu }^{1loop}{)}_{\alpha \alpha }$$

$$m_{\nu }=-\theta \tilde{M}{\theta }^T$$

$$\tilde{M}=\big{[}1-\frac{1}{v^{2}}M_{M}M_{N}^{\rm diag}l(M_{N}^{\rm diag})\big{]}M_{M}\,.$$

$${\tilde{M}}_{IJ}\simeq {\tilde{M}}_{IJ}^{\mathrm{diag}}=M_I{\delta }_{IJ}\left(1-\frac{M_I^2}{v^2}l(M_I)\right).$$

$$\Theta =i{U}_{\nu }\sqrt{m_{\nu }^{\text{diag}}}\mathcal{R}{\sqrt{{\tilde{M}}^{\text{diag}}}}^{-1},$$

$$\mathcal{R}={\mathcal{R}}^{23}{\mathcal{R}}^{13}{\mathcal{R}}^{12},$$

Full text

††journal: Eur. Phys. J. C

1. 1 Introduction

2. 1.1 Motivation

3. 1.2 Main improvements compared to previous studies

4. 2 Right-handed neutrino physics

5. 2.1 Basic definitions

6. 2.2 The seesaw limit

7. 2.3 The role of the matrix $U_{N}$

8. 2.4 Casas-Ibarra parametrisation

9. 2.5 The symmetry protected scenario

10. 2.6 Connection to the model with $n=2$

11. 3 Observables, experiments and likelihoods

12. 3.1 Active neutrino mixing

13. 3.2 Indirect constraints

14. 3.2.1 Electroweak precision observables

15. 3.2.2 Lepton flavour violation

16. 3.2.3 Lepton universality

17. 3.2.4 CKM unitarity

18. 3.2.5 Neutrinoless double-beta decay

19. 3.2.6 Big Bang Nucleosynthesis

20. 3.3 Direct RHN searches

21. 3.3.1 PIENU

22. 3.3.2 PS-191

23. 3.3.3 CHARM

24. 3.3.4 E949

25. 3.3.5 NuTeV

26. 3.3.6 DELPHI

27. 3.3.7 ATLAS

28. 3.3.8 CMS

29. 3.3.9 LHCb

30. 3.3.10 Other experiments

31. 4 Scanning strategy and parameter ranges

32. 4.1 Parameters and priors

33. 4.2 Targeted scans

34. 4.3 Scanning framework

35. 4.4 Data postprocessing

36. 4.5 Capped likelihood

37. 5 Results and discussion

38. 5.1 General constraints on the RHN mass and mixing

39. 5.2 Discussion of individual bounds

40. 5.3 Lightest neutrino mass and flavour mixing

41. 5.4 Discussion of excesses likelihoods

42. 6 Conclusions & Outlook

43. A GAMBIT Implementation

44. A.1 Neutrino models

45. A.2 Right-handed neutrino likelihood functions

46. B Full expressions for the relevant observables

47. B.1 Decay widths and form factors for LFV observables

48. B.2 Decay widths relevant for Big Bang Nucleosynthesis

49. C Distinguishing symmetry protected from tuned parameter choices

50. D Partial likelihoods

1 Introduction

1.1 Motivation

The observation of neutrino flavour oscillations is one of the strongest hints for the existence of particle physics beyond the Standard Model (SM). The oscillations imply that neutrinos have small masses, while the minimal SM predicts that they are massless. At the same time neutrinos are the only elementary fermions that are only known to exist with left handed chirality $\nu_{L}$. If right handed neutrinos $\nu_{R}$ exist, one could immediately add a Dirac mass term $\bar{\nu_{L}}M_{D}\nu_{R}$ to the SM Lagrangian in analogy to all other known fermions. The fact that the $\nu_{R}$ have not been seen yet could easily be explained because they are "sterile", i.e., not charged under any known gauge interactions. The same property also makes it possible for them to have a Majorana mass term $\bar{\nu_{R}}M_{M}\nu_{R}^{c}$ in addition to the Dirac mass. For eigenvalues of $M_{M}$ that are much larger than the observed light neutrino masses, the smallness of the neutrino masses can be explained via the seesaw mechanism Minkowski:1977sc ; Mohapatra:1979ia ; Mohapatra:1980yp ; GellMann:1980vs ; Yanagida:1980xy . Neutrino oscillation data is, however, not sufficient to pin down the value of $M_{M}$, known as seesaw scale, because it is primarily sensitive to the combination $M_{D}M_{M}^{-1}M_{D}^{T}$. The range of allowed values spans from a few eV deGouvea:2005er up to the scale of Grand Unification Asaka:2015eda . For specific choices of their Majorana mass the $\nu_{R}$ could in addition solve a number of long standing problems in cosmology. For instance, they could explain the baryon asymmetry of our Universe via leptogenesis during the decay Fukugita:1986hr or production Akhmedov:1998qx ; Asaka:2005pn of the heavy neutrinos or provide a viable dark matter candidate Dodelson:1993je ; Shi:1998km . An overview of the cosmological implications of different choices of $M_{M}$ can e.g. be found in Ref. Drewes:2013gca .

Experiments can directly search for heavy neutrinos if $M_{M}$ is below the TeV scale. Such searches have been performed in various different facilities, including high energy colliders and fixed target experiments. This is the mass range we consider in the present article. In addition, the $\nu_{R}$ would indirectly affect precision observables or searches for rare processes. A summary of different existing constraints can be found in the reviews Atre:2009rg ; Boyarsky:2009ix ; Drewes:2013gca ; Deppisch:2015qwa ; Cai:2017mow . For the future a wide range of different searches have been proposed, an overview can be found in Refs. Antusch:2016ejd ; Beacham:2019nyx ; Alimena:2019zri . In order to decide about the best possible search strategy is it important to understand which parameter region is already ruled out by past experiments. This is in fact a non-trivial question because different observables are correlated in the seesaw model, and the requirement to simultaneously respect all known experimental results imposes stronger constrains on the model parameter space than superimposing individual bounds. Such global constraints can only be derived within a given model. An important quantity in this context is the unknown number $n$ of right handed neutrino flavours. The minimal number that is required to explain the light neutrino oscillation data is $n=2$, which would necessarily require the lightest SM neutrino to be massless. The minimal number that is required to generate masses for all three SM neutrinos is $n=3$. This choice is also somewhat appealing in view of the fact that there are three fermion generations in the SM, and it is mandatory for anomaly freedom in many gauge extensions of the SM. The goal of the present work is to impose global constraints on the parameter space of the model with $n=3$, based on the combination of direct, indirect and cosmological constraints summarised in section 3.

Several authors have previously imposed global constraints on the properties of right handed neutrinos. Here we exclusively focus on models in which the right handed neutrinos can explain the light neutrino oscillation data.111The authors of Ref. deGouvea:2015euy considered a single heavy neutrino, but made the conservative assumption that this particle may predominantly decay into a dark sector via new interactions. This e.g. excludes most sterile neutrino Dark Matter models because the feeble coupling of such particles that is required to ensure their longevity implies that its contribution to the light neutrino mass generation can be neglected Boyarsky:2006jm .222We refer the reader to Refs. Adhikari:2016bei ; Boyarsky:2018tvu for recent reviews on sterile neutrino Dark Matter. One of the most complete studies of indirect constraints on the parameter space for $n=2$ in the last few years was presented in Ref. Antusch:2014woa , where multiple electroweak precision observables and flavour-violating decays were included, along with tests of lepton universality and the unitarity of the CKM matrix. Loop corrections to some of these relations were considered in Ref. Fernandez-Martinez:2015hxa . The authors of Drewes:2016jae included direct search constraints and those from big bang nucleosynthesis. The model with $n=3$ is much less studied. Recent analyses of indirect constraints include Refs. Escrihuela:2015wra ; Fernandez-Martinez:2016lgt , direct search constraints and BBN have been added to this in Ref. Drewes:2015iva .

1.2 Main improvements compared to previous studies

In this paper, we present the first full frequentist analysis of the $n=3$ right-handed neutrino (RHN) extension of the SM, for a wide range of RHN masses from about 60 MeV to 500 GeV. We opted for a frequentist analysis rather than a Bayesian analysis since this is best suited to fully explore the valid parameter space while avoiding prior dependence and volume effects of the parameter space (however, we emphasize that we do not perform a full sampling-based goodness-of-fit analysis and instead resort for practical reasons to an approximate treatment of likelihood and their sampling statistics). We improve on different aspects of earlier analyses by combining all the strongest limits exerted by experiments as well as indirect signatures in a statistically consistent manner. Previous studies that examined the parameter space for $n=3$ either used a subset of the constraints included here Escrihuela:2015wra ; Fernandez-Martinez:2016lgt or used less rigorous statistical methods Drewes:2015iva and focused on specific regions of the parameter space Gorbunov:2014ypa .

• •

While most previous studies fixed the mixing angles and mass differences in the active neutrino sector to the best fit values as presented in Esteban:2016qun , we take into account likelihoods for the active neutrino observables.

• •

Electroweak observables require precise calculations for its comparison with the extremely accurate measurements. We therefore use the calculation of the SM prediction for $\sin\theta_{w}^{eff}$ up to two-loop order Ferroglia:2012ir .

• •

Most studies of lepton flavour violation in neutrino models focus exclusively on the most constraining processes, such as $\mu\to e\gamma$ and $\mu\to eee$ Antusch:2014woa ; Drewes:2015iva . In this work we include all lepton flavour violating processes, in particular all leptonic $\tau$ decays, for which we use the most recent average of experimental results provided by HFLAV Amhis:2016xyh , as well as $\mu-e$ conversion in nuclei (Pb, Au and Ti).

• •

For neutrinoless double-beta decay, in comparison with Drewes:2015iva , we opt to carry out our analysis conservatively; in addition, the upper limit on the effective Majorana mass and hence the mixing is encoded in the form of a (one-sided) Gaussian likelihood, not as a strict cut.

• •

Lepton universality tests are often centered on leptonic decays of mesons, $K$ and $\pi$, $\tau$-leptons and $W$-bosons Antusch:2014woa . We supplement these tests of universality with the recently observed semileptonic decays of B-mesons Aaij:2014ora ; Aaij:2015yra ; Aaij:2017vbb .

• •

We improve the treatment of CKM unitarity with respect to the discussion in Ref. Drewes:2015iva .

• •

Concerning direct searches, previous studies have used only a subset of the experiments considered here Ruchayskiy:2011aa ; deGouvea:2015euy , or chose to place a hard cut at the upper limits presented in the individual papers Drewes:2015iva ; Drewes:2016jae . We implement the strongest constraints over the mass range as likelihoods. The statistical combination of these likelihoods also leads to more accurate profile likelihood contours in comparison to simply overlaying individual limits.

• •

We study in detail the flavour mixing pattern of the three RHN, for different values of the lightest neutrino mass. We discuss the limit where the lightest neutrino is massless and the connection to the $n=2$ case.

We use here the open-source software package GAMBIT gambit . It includes an interface to Diver ScannerBit , a differential evolution-based scanner that provides efficient sampling performance for frequentist scans.

This paper is organised as follows. In section 2, the model, parametrisation used and essential quantities are defined. All the observables and experiments that are considered are subsequently discussed in detail in section 3. Our scanning strategy, parameter ranges and applied priors are mentioned in section 4. The results are presented in section 5 and we discuss the implications of the combined constraints for future searches in section 6. In Appendix A we comment on the details of the implementation in GAMBIT, in Appendix B we explicitly give the expressions for the different observables, in Appendix C we provide details on how we interpret our results in view of the criterion of technical naturalness, and in Appendix D we show the different partial likelihoods.

2 Right-handed neutrino physics

2.1 Basic definitions

The addition of three RHNs to the particle content of the Standard Model introduces in total 18 new parameters. In this section we summarise basic relations in the seesaw model and define our notation, following Ref. Drewes:2015iva .

The most general renormalisable Lagrangian that can be constructed from SM fields and the $\nu_{R}$ has the following form:

[TABLE]

Hereby, $\ell_{L}=(\nu_{L},e_{L})^{T}$ indicate the left-handed leptons333Throughout this article we use four component spinor notation, where the chiral spinors $\nu_{R}$ and $\ell_{L}$ have only two non-zero components ($P_{R}\nu_{R}=\nu_{R}$ and $P_{L}\ell_{L}=\ell_{L}$). As a result, no explicit chiral projectors are necessary in the weak interaction term (21). of the SM and $\Phi$ is the Higgs doublet with $\tilde{\Phi}=\epsilon\Phi^{\ast}$ and $\epsilon$ being the Levi-Civita tensor. $M_{M}$ is the Majorana mass matrix for $\nu_{R}$ and $F$ is the Yukawa coupling matrix. We work in a flavour basis where $M_{M}=\text{diag}(M_{1},M_{2},M_{3})$.

After electroweak symmetry breaking (EWSB), the complete neutrino mass term reads

[TABLE]

with

[TABLE]

where $M_{D}=Fv$, $v$ being the Higgs vacuum expectation value ($v=174$ GeV in the ground state). We include the one loop corrections $\delta m_{\nu}^{1loop}$ and ${\delta}M_{N}^{1loop}$ as we aim for performing an analysis to be consistent at second order in the Yukawa couplings $F$. The mass matrix (3) can be diagonalised by a matrix of the form Fernandez-Martinez:2015hxa

[TABLE]

with

[TABLE]

Hereby, $\theta$ indicates the matrix that mediates the mixing between the active neutrinos $\nu_{L}$ and the sterile neutrinos $\nu_{R}$. We can generally write

[TABLE]

with

[TABLE]

The additional complex conjugation of $U_{N}$ ensures that the relation among mass and flavour eigenstates will be analogous for left-handed neutrinos (LHNs) and RHNs within the notation. In the second relation in eq. (8) we have neglected the difference between the eigenvalues of $M_{M}$ and $M_{N}$, which is of second order in $\theta$. This is justified for the present purpose because of the experimental constraints on the magnitude of the elements $\theta_{\alpha I}$, which we discuss further below.

2.2 The seesaw limit

The limit of small $\theta_{\alpha I}$ is usually referred to as the seesaw limit, it corresponds to $M_{D}\ll M_{M}$ (in terms of eigenvalues). It allows the approximation

[TABLE]

and

[TABLE]

leading to

[TABLE]

with

[TABLE]

and

[TABLE]

The loop correction to the light neutrino mixing matrix is given by Pilaftsis:1991ug :

[TABLE]

where $l(M_{I})$ is a loop function given by

[TABLE]

The light and heavy neutrino mass eigenstates are described by the flavour vectors

[TABLE]

and

[TABLE]

respectively. We can define the matrices $V_{\nu}$ and $V_{N}$ that represent the mixing between mass and interaction eigenstates in the respective sectors as

[TABLE]

while mixing between the two sectors is encoded in the matrix

[TABLE]

This quantity is of primary interest because it controls the interactions of the heavy neutrinos with the physical Higgs field $h$ and the gauge bosons $W$ and $Z$,

[TABLE]

Here $g$ is the weak gauge coupling constant and $\theta_{W}$ the Weinberg angle. For convenience, we introduce the notation

[TABLE]

From the relations (3) and (7) it is straightforward to derive the relation

[TABLE]

2.3 The role of the matrix $U_{N}$

In our numerical scan we approximate $U_{N}$ by unity.444Note that the approximation $U_{N}=\mathbb{I}$ also allows to neglect $\delta M_{N}^{1loop}$ because it only amounts to a change in the matrix $U_{N}$ Drewes:2019mhg . For generic parameter choices this can be justified because we work in a basis where $M_{N}$ is diagonal, and the physical mass matrix (13) is also diagonal up to corrections of second order in $\theta$. These corrections can lead to a large deviation of $U_{N}$ from unity only if the eigenvalues of $M_{M}$ are quasi-degenerate, so that the $\mathcal{O}[\theta^{2}]$ terms in the matrix (13) are relevant.

If a degeneracy between only two of the RHNs is caused by a symmetry, cf. sec. 2.5, then it can be shown that the effect of $U_{N}$ on the $U_{\alpha I}^{2}$ is small even if individual entries of $U_{N}$ are larger than the $U_{\alpha I}^{2}$ Drewes:2019mhg . This means that the production cross sections for heavy neutrinos are not affected. However, the branching ratio between lepton number violating (LNV) and lepton number conserving heavy neutrino decays is affected by $U_{N}$ Drewes:2019byd . This has no big effect on our scan because constraints from searches for LNV are sub-dominant in almost the entire mass range that we consider, but it may have important implications for future searches.

$U_{N}$ can have a big impact on the individual mixings $U_{\alpha I}^{2}$ of each heavy neutrino if all three Majorana masses are degenerate. This can be accommodated in technically natural scenarios discussed in the following section 2.5, cf. in particular footnote 6. The practical impact on experimental searches is, however, limited because most experiments are not able to kinematically resolve small mass splittings and therefore only probe $U_{\alpha}^{2}$ in this regime (rather than the couplings $U_{\alpha I}^{2}$ of individual heavy neutrino flavours). Also in this case observables that are sensitive to LNV are the only ones that are likely to be affected.

Finally, if the degeneracy between the heavy neutrino masses is accidental, then the proof in Ref. Drewes:2019mhg does not apply, and $U_{N}$ can have a significant effect on the $U_{\alpha I}^{2}$ even if only two heavy neutrinos have degenerate masses. Our results contain a significant number of points of this kind because we performed several scans with "agnostic" parameter ranges that do not suppress fine-tuned points, cf.  table 5. However, the fact that experiments are unlikely to resolve the individual resonances in this case implies that they are only sensitive to the quantities $U_{a}^{2}$, where the summation is to be taken over the mass degenerate heavy neutrino flavours only. As in the previous two cases, the effect of $U_{N}$ on the total production rate is minor because the matrix mainly re-distributes coupling between the mass degenerate states. The main affect would again be on LNV observables.

In summary, if any heavy neutrinos are discovered in the future, a comparison between the branching ratios of lepton number violating and lepton number conserving decays will give important insight into the mechanism of neutrino mass generation and will be crucial to identify any underlying symmetries.

2.4 Casas-Ibarra parametrisation

In the current work, we use the Casas-Ibarra (C-I) parametrisation Casas:2001sr , generalised to include the 1-loop correction to the left-handed neutrino mass matrix Lopez-Pavon:2015cga . This provides a simple way to impose constraints from light neutrino oscillation data in our scan. This parametrisation is based on the observation that $m_{\nu}$ in eq. (12) can be expressed as

[TABLE]

with

[TABLE]

Since the loop function is smooth we can neglect the difference in the eigenvalues of $M_{M}$ and $M_{N}$,

[TABLE]

In this scheme the sterile neutrino mixing matrix, i.e. the matrix encoding the mixing among LHNs and RHNs (20) can be written as

[TABLE]

where $U_{\nu}$ is the PMNS matrix introduced above, $m_{\nu}^{\text{diag}}$ is the diagonalised, one-loop-corrected LHN mass matrix and $\tilde{M}^{\text{diag}}$ is the analogous RHN mass matrix, given by (28). Furthermore, $\mathcal{R}$ is a complex, orthogonal matrix that is parametrised by complex angles $\omega_{ij}$

[TABLE]

where $\mathcal{R}^{ij}$ has the non-zero elements

[TABLE]

Since we work in the flavour basis in which the Yukawa couplings of the charged leptons are diagonal, $U_{\nu}$ can be parametrised as

[TABLE]

where $U_{\pm\delta}=\mathrm{diag}(e^{{\mp}i\delta/2},1,e^{{\pm}i\delta/2})$ and $V^{ij}$, parametrised by the LHN mixing angles $\theta_{ij}$, has non-zero elements analogous to $\mathcal{R}$. Furthermore, $\alpha_{1}$, $\alpha_{2}$ and $\delta$ are CP-violating phases.

The C-I parametrisation scheme generates by construction Yukawa couplings and mixing angles $\Theta$ that are consistent with light neutrino oscillation data up to second order in $\theta$. This has two disadvantages. First, one may find it unsatisfactory that we treat light neutrino oscillation data differently from other constraints. Second, the C-I is a “bottom up” parametrisation. There is usually no simple relation between the C-I parameters and parameters that may be well-motivated from a model building viewpoint, and any theory-motivated prior on the RHNs’ mixings and masses would acquire a rather convoluted form in the C-I parametrisation. In particular, there is no simple way to distinguish “natural” from “fine tuned” parameter choices. Hence, we refrain from performing Bayesian scans in the current work, and instead concentrate on a likelihood-based frequentist treatment. In view of the high dimensionality of the parameter space and the complicated functional form of the different constraints, the disadvantages of the C-I parametrisation are, however, compensated for by the numerical advantage that one gains.

2.5 The symmetry protected scenario

The smallness of the light neutrino masses $m_{i}$ can be explained in different ways by the seesaw relation (26). One possibility is that the $N_{I}$ are very heavy, i.e., $M_{I}\gg v$, in which case the smallness of $m_{i}$ is due to the smallness of the ratio $v/M_{I}$. This choice for the mass scale(s) $M_{I}$ is well-motivated by Grand Unified Theories,555See Fukugita:2003en ; Croon:2019kpe for a review on neutrino masses in the context of Grand Unified Theories. but raises the question of radiative corrections to the Higgs potential from the Yukawa couplings of the RHNs Vissani:1997ys .

This “hierarchy problem” can be avoided in low scale seesaw scenarios. Low values of $M_{I}$ are natural because in the limit $M_{I}\to 0$ the $B-L$ symmetry in the SM is restored. In this case, however, the smallness of $m_{i}$ can no longer be explained efficiently by the suppression of $v/M_{I}$, as it typically requires couplings

[TABLE]

that are very small, in particular for seesaw scales as low as 100 MeV.

Such small values for fundamental parameters are considered ’unnatural’ by many theorists Giudice:2008bi , though some possible explanations have been proposed Froggatt:1978nt . However, this estimate relies on the underlying assumption that there are no cancellations (accidental or otherwise) in the seesaw relation (26), which would allow for much larger $U_{\alpha I}^{2}=|\Theta_{\alpha I}|^{2}$ than the naive estimate (35) suggests while keeping the eigenvalues $m_{i}^{2}$ of $m_{\nu}^{\dagger}m_{\nu}$ small.

Hence, a technically natural tHooft:1979rat way to obtain small neutrino masses $m_{i}$ can be realised if the Lagrangian (2.1) approximately respects a $B-\bar{L}$ symmetry Shaposhnikov:2006nn ; Kersten:2007vk (cf. also Gluza:2002vs ), where $\bar{L}$ is a generalised lepton number under which combinations of the $\nu_{Ri}$ are charged. Such $B-\bar{L}$ symmetry is exact if the Yukawa coupling and mass matrix take the form Moffat:2017feq

[TABLE]

in which case the light neutrinos are exactly massless $m_{i}=0$. In order to generate non-zero light neutrino masses this symmetry has to be slightly broken, i.e.,

[TABLE]

where the entries of the matrices $\upmu$ and $\upepsilon$ are small symmetry breaking parameters.

If the symmetry is not exact $M_{M}$ can have off-diagonal elements, see for example Ref. Abada:2007ux . Throughout this work we use a basis in which $M_{M}$ is diagonal. The diagonalisation affects the form of the Yukawa matrix $F$, but as long as the off diagonal elements of $\upmu$ are small, this only leads to a small modification of the flavour structure. For the following discussion we will therefore adapt the simpler form Abada:2018oly 666 An important exception is the case $\mu\ll 1$, $\bar{M}^{\prime}\simeq\bar{M}$. In that situation even small off-diagonal elements $\upmu_{ij}$ can lead to a comparably large misalignment between the basis in which $F$ has the form (2.5) and the heavy neutrino mass basis, which means which that all heavy neutrinos have unsuppressed Yukawa couplings $\sim F_{a}$ in spite of the fact that $\epsilon^{\prime}_{a}\ll 1$, cf. ref. Drewes:2019byd for a discussion. However, in this case all three mass eigenstate $N_{I}$ have approximately the same mass $\bar{M}$ and cannot be distinguished kinematically. In this case the experimentally relevant mixing is $U_{a}^{2}$, the magnitude of which is controlled by the large entries $F_{a}$. Heavy neutrino oscillations in the detector Boyanovsky:2014una ; Cvetic:2015ura ; Anamiati:2016uxp ; Dib:2016wge ; Das:2017hmg ; Antusch:2017ebe ; Antusch:2017pkq ; Cvetic:2018elt ; Hernandez:2018cgc ; Cvetic:2019rms could provide an indirect way to access the small mass splitting and phenomenologically study this specific case.

[TABLE]

with $\epsilon_{\alpha}^{\prime},\epsilon_{\alpha},\mu,\ll 1$ being small symmetry breaking parameters and $F_{\alpha}$ being of the order of one. This means that one heavy neutrino practically decouples while the other two approximately form a Dirac spinor with mass $\bar{M}$.

In this symmetry protected scenario there is no upper limit on $U_{\alpha I}^{2}$ from neutrino oscillation data. In the mass range considered here the upper limit comes from the experimental constraints, while for larger masses there is a theoretical bound $U_{\alpha I}^{2}<4\pi(n-1)(v/\bar{M})^{2}$ from the requirement that the Yukawa couplings remain perturbative Asaka:2015eda . This provides a theoretical motivation for a low scale seesaw with experimentally accessible mixings $U_{\alpha I}^{2}$. Specific examples that motivate this limit include “inverse seesaw”  Wyler:1982dd ; Mohapatra:1986aw ; Mohapatra:1986bd ; Bernabeu:1987gr , “linear seesaw” Akhmedov:1995ip ; Akhmedov:1995vm , scale invariant Khoze:2013oga and some technicolour-type models Appelquist:2002me ; Appelquist:2003uu and also the $\nu$MSM Asaka:2005pn ; Shaposhnikov:2006nn .

2.6 Connection to the model with $n=2$

The parametrisation (2.5) suggests that the $B-\bar{L}$ symmetric limit for the model with $n=3$ should contain the model with $n=2$, as the third heavy neutrino decouples for $\epsilon^{\prime}_{a}\to 0$. This is, for example, observed in the $\nu$MSM. However, some care is required when taking this limit if one wants to be consistent with neutrino oscillation data.

First, it is clear that not all seven symmetry breaking parameters $\epsilon_{a},\epsilon^{\prime}_{a},\mu$ can be set to zero because this would give exactly massless light neutrinos. Which of these parameters are non-zero and how small they are with respect to each other depends on the way how the symmetry is broken and thus on the particle physics model in which the Lagrangian (2.1) is embedded. It is not possible to make a model independent statement about the relative size of the $\epsilon^{\prime}_{a}$ in relation to other model parameters.

Second, the parametrisation (2.5) is not the most general one: If we allow for small off diagonal elements in the general form (37), then all three heavy neutrinos can have unsuppressed interactions if $\bar{M}^{\prime}\simeq\bar{M}$, cf. footnote 6. Hence, if $\bar{M}^{\prime}$ is degenerate with $\bar{M}$, one cannot expect to recover the $n=2$ model even if $\epsilon_{\alpha}\ll 1$.

Finally, as discussed in more detail in Appendix C, there are Casas-Ibarra parameter choices that yield small values of $m_{\nu_{0}}$, but correspond to highly fine-tuned scenarios where this smallness is due to accidental cancellations. These solutions can imitate the symmetry protected scenario and can also circumvent the seesaw upper limit and thus reach high values of $U_{\alpha I}^{2}$.

3 Observables, experiments and likelihoods

Models with heavy right-handed neutrinos, as described above, will alter the SM predictions for different observables that are already significantly constrained by experimental results. In this analysis, we implemented all relevant constraints such as active neutrino likelihoods (3.1) and direct detection experiments which currently exert the strongest bounds over the considered mass range (3.3); these include beam dump and peak search experiments, which looked for RHNs in meson, tau and gauge boson decays. Besides, we similarly include the most relevant indirect constraints: electroweak precision observables (3.2.1), lepton flavour violating processes (3.2.2), lepton universality constraints (3.2.3), BBN (3.2.6), neutrinoless double-beta decay (3.2.5) and CKM unitarity (3.2.4).

In this section, we will focus on the physics and statistics aspects of our likelihood functions. The corresponding implementation of GAMBIT capabilities and module functions associated with the various observables are discussed in detail in Appendix A.

3.1 Active neutrino mixing

In contrast to previous studies, we include likelihoods for the active neutrino mixing observables in our analysis: the three mixing angles $\theta_{12},\theta_{13},\theta_{23}$, the mass splittings $\Delta m^{2}_{21}$ and $\Delta m^{2}_{3\ell}$ with $\ell=1$ for normal ordering and $\ell=2$ for inverted ordering, as well as the CP-phase $\delta_{\mathrm{CP}}$. We use the most recent publically available results of the global analysis of solar, atmospheric, reactor and accelerator neutrino data in the framework of three neutrino oscillations provided by the NuFIT collaboration (as of January 2018) Esteban:2016qun ; NuFit , including

• •

the solar neutrino experiments Homestake chlorine Cleveland:1998nv , Gallex/GNO Kaether:2010ag and SAGE Abdurashitov:2009tn , SNO Aharmim:2011vm , the four phases of Super-Kamiokande Hosaka:2005um ; Cravens:2008aa ; Abe:2010hy and two phases of Borexino Bellini:2011rx ; Bellini:2008mr ; Bellini:2014uqa ,

• •

the atmospheric experiments IceCube/DeepCore Aartsen:2014yll ,

• •

the reactor experiments KamLAND Gando:2013nba , Double-Chooz An:2016srz , Daya-Bay An:2016ses and Reno reno ,

• •

the accelerator experiments MINOS Adamson:2013whj ; Adamson:2013ue , T2K t2k and NO$\nu$A nova ,

• •

the cosmic microwave background measurement Planck Ade:2015xua

For our global fit, we take the provided one-dimensional $\Delta\chi^{2}$ tables for both orderings of the NuFIT collaboration NuFit . For more detailed information, we refer to Esteban:2016qun and references therein. We emphasize that using higher dimensional tables that account for correlations would in general lead to (slightly) more stringent results on the RHN parameter space, hence our treatment can be considered as conservative.

3.2 Indirect constraints

3.2.1 Electroweak precision observables

The leptonic charge currents are modified by the RHNs, and hence the value of $G_{\mu}$ that is measured via the muon decay will differ from the actual Fermi constant $G_{F}$ which is defined in terms of the fine structure constant and mass of the $Z$ boson. The correction can be written as Drewes:2015iva

[TABLE]

and is caused by the non-unitarity of the flavour mixing matrix $V_{\nu}$, see Eq. (18), which leads to a slight suppression of the muon decay.

Both the weak mixing angle $\theta_{w}$ and the mass of the $W$ boson $m_{W}$ depend on $G_{\mu}$ at one loop, which means they also get a correction from the active-sterile mixing matrix $\Theta$, which is given by Antusch:2014woa

[TABLE]

where $s_{w}^{2}=\sin^{2}\theta_{w}$. Since experiments typically measure the effective Weinberg angle $s^{2}_{eff}$, and assuming the QCD corrections factorize from the leptonic corrections Antusch:2015mia , we use for the SM prediction the highly accurate calculation, including corrections up to two-loops, from Ferroglia:2012ir

[TABLE]

Other electroweak precision observables affected by the presence of the heavy neutrinos are the decays of the $Z$ and $W$ bosons, in particular the invisible decay width of the $Z$ boson, $\Gamma_{\rm{inv}}$, and the leptonic decays of $W$. Under the assumption that the radiative corrections factorize from the heavy neutrino contribution, at least up to order $\theta^{2}$ Antusch:2015mia ; Fernandez-Martinez:2015hxa , one can write the invisible decay width of the Z as Abada:2013aba

[TABLE]

where we have neglected the contribution from $Z\to N_{i}N_{j}$ due to being of order $\theta^{4}$, and for the SM decay $Z\to\nu_{i}\nu_{j}$ we use the 2-loop calculation from Dubovyk:2018rlg .

The contribution of heavy neutrinos to the $W$ decay widths to leptons can be written as Antusch:2014woa

[TABLE]

where we defined $x_{\alpha}\equiv m_{l_{\alpha}}^{2}/m_{W}^{2}$.

We construct Gaussian likelihoods for these observables using the experimental measurements and uncertainties displayed in Table 1. All these observables depend on $G_{\mu}$ (eq.(39)) either directly or through another observable ($s_{w}$ or $m_{W}$). Since the experimental measurements of these quantities are independent of each other, we assume them to be uncorrelated.

3.2.2 Lepton flavour violation

Flavour changing neutral processes, such as lepton flavour violation (LFV), are strongly suppressed in the Standard Model at one loop due to the GIM mechanism Glashow:1970gm . Hence, any non-trivial contribution to these processes from physics beyond the Standard Model would dominate over the SM contribution, which in turn makes the experimental determination of these observables a smoking gun of new physics. Several experiments have attempted to measure LFV processes with outstanding precision and they have imposed a set of upper limits on their branching fractions. In Table 2 we list the most significant of these observables, along with the experimental upper bound on their branching ratios and the experiment that provided it.

The experimental upper bounds for LFV $\mu$ and $\tau$ decays in Table 2 are given as branching fractions with respect to the total decay width of the respective lepton PDG17 ; Ilakovac:1994kj ,

[TABLE]

In the model with three heavy neutrinos the leading contributions to these observables arise from dipole and box diagrams with mixing between the active and sterile neutrinos, given by the active-sterile mixing matrix $\Theta$. The relevant LFV processes containing these diagrams are of the form $l_{\alpha}^{-}\to l_{\beta}^{-}\gamma$, $l_{\alpha}^{-}\to l_{\beta}^{-}l_{\beta}^{-}l_{\beta}^{+}$, $l_{\alpha}^{-}\to l_{\beta}^{-}l_{\gamma}^{-}l_{\gamma}^{+}$ and $l_{\alpha}^{-}\to l_{\gamma}^{-}l_{\gamma}^{-}l_{\beta}^{+}$. The associated decay widths can be found in Appendix B.1.

Lastly, LFV processes can result in a neutrinoless $\mu-e$ conversion inside a nucleus. Muons captured by a nucleus typically decay in orbit providing a continuous spectrum of energy for the electron in the final state. In coherent flavour violating conversion, $\mu^{-}N\to e^{-}N$, final state electrons have a discrete energy spectrum, corresponding to the mass of the decaying muon. Consequently experiments measure the rate at which this conversion happens, with respect to the rate of capture by the nucleus,

[TABLE]

The corresponding expressions for the conversion ratio, as well as the nuclear parameters for the two nuclei studied, $\text{Ti}_{22}^{48}$, $\text{Au}_{79}^{197}$ and $\text{Pb}_{82}^{208}$, can be found in Appendix B.1.

The likelihoods for these LFV observables are all Gaussian upper limit likelihoods. They are computed as

[TABLE]

using the experimental data from Table 2. More specifically, we assume a measured value of $x_{0}$ for all observables777In the cases where the experiments do not provide a measured value we take $x_{0}=0$., and set $\sigma=v/1.64$ for full Gaussians and $\sigma=v/1.28$ for one-sided Gaussians, where $v$ is the quoted upper 90% C.L. limit.

3.2.3 Lepton universality

Recent measurements of meson decays Aaij:2014ora ; Aaij:2015yra ; Aaij:2017vbb have put into question the flavour-independence of leptonic charged currents, as predicted by the SM. Previous tests of lepton universality performed by LEP and SLC, using lifetime measurements of the tau and muon as well as the partial decay widths of the Z boson, showed no such deviation. This has lead to the formulation of many BSM theories attempting to explain the deviation shown in meson decays with sterile neutrinos Abada:2013aba ; Boucenna:2015raa ; Bryman:2019ssi .

The presence of right-handed neutrinos modifies the leptonic currents and thus triggers a contribution to processes testing lepton universality such as in the fully leptonic decays of charged mesons, $X^{+}\to l^{+}\nu$, or the semileptonic decays of $B$ mesons $B^{0/\pm}\to X^{0/\pm}l^{+}l^{-}$.

In order to cancel the considerable hadronic uncertainties present in the decays of pseudoscalar mesons, lepton universality tests are best formulated using ratios between lepton species. For fully leptonic and semileptonic decays of mesons, these ratios are expressed as

[TABLE]

respectively.

In case of fully leptonic decays, one can express the test of lepton universality in terms of deviations from the SM prediction as

[TABLE]

where the sterile neutrino contribution can be calculated from the active-sterile mixing matrix $\Theta$ as Shrock:1980ct ; Drewes:2015iva

[TABLE]

where we used

[TABLE]

with $\vartheta$ being the Heaviside step function, $r_{\alpha}\equiv m_{l_{\alpha}}^{2}/m_{X}^{2}$ and $r_{I}\equiv M_{I}^{2}/m_{X}^{2}$. The SM predictions used in eq. (49) for the tests of lepton universality for pions and kaons are $R^{\pi}_{e\mu,SM}=1.2354\times 10^{-4}$ and $R^{K}_{e\mu,SM}=2.477\times 10^{-5}$, respectively Cirigliano:2007xi .

The contribution from heavy right-handed neutrinos to the semileptonic decays of $B$ mesons is much less significant than to the leptonic decays. As argued in Ref. Abada:2013aba , the effect on B decays to charmed mesons, $B^{\pm}\to Dl\nu$, is completely negligible. Semileptonic decays to $K$ mesons are more affected, particularly the decays $B^{+}\to K^{+}l^{+}l^{-}$ and $B^{0}\to K^{*0}l^{+}l^{-}$. Assuming that $m_{l}\ll m_{K^{(*)}}$ and that the Wilson coefficient $C_{7}\ll C_{9},C_{10}$, one can approximate the ratios $R_{K}$ and $R_{K^{*}}$ as Ghosh:2014awa

[TABLE]

and the BSM contributions to the Wilson coefficients $\Delta C_{9}^{\alpha}$ and $\Delta C_{10}^{\alpha}$ can be expressed as He:2017osj

[TABLE]

with $x_{t}=m_{t}^{2}/m_{W}^{2}$, $x_{I}=M_{I}^{2}/m_{W}^{2}$ and the loop function

[TABLE]

NNL calculations for the Standard Model contribution to the Wilson coefficients $C_{9}$ and $C_{10}$ used in Eq. (52) gives $C_{9}^{SM}=4.211$ and $C_{10}^{SM}=-4.103$ Altmannshofer:2008dz ; Choudhury:2017ijp .

In addition to meson decays, other common tests of lepton universality include the decays of the $W$ boson to leptons as well as $\tau$ decays. The ratio of decay widths of $W$ to charged leptons $l_{\alpha}$ and $l_{\beta}$ can be written as Antusch:2014woa

[TABLE]

Deviations from the SM for the lepton universality test in $\tau$ decays follow the same form as in Eq. (50) and the SM prediction is $R^{\tau}_{\mu e,SM}=0.973$ Pich:2009zza .

These tests of lepton universality are implemented as Gaussian likelihoods centered on the experimentally measured value. The experimental measurements, with their corresponding uncertainties888The experimental uncertainties for $R^{B}_{K^{(*)}}$ are obtained as the sum in quadrature of the statistical and systematic uncertainties provided by Aaij:2014ora and Aaij:2017vbb ., are shown in Table 3. The measurements of $R^{\pi}_{e\mu}$ include subleading decays with $\gamma$’s, hence the upper limit shown is the PDG average of the ratios of $\Gamma(\pi^{+}\to l_{\alpha}^{+}\nu_{\alpha})+\Gamma(\pi^{+}\to l_{\alpha}^{+}\nu_{\alpha}\gamma)$, based on the measurements in  Britton:1992pg ; Czapek:1993kc ; Aguilar-Arevalo:2015cdf . Two experimental measurements are shown for $R_{K^{*}}$ corresponding to two regions of the dilepton invariant mass $0.045<q^{2}<1.1(\text{GeV}^{2}/c^{4})$ for (1) and $1.1<q^{2}<6.0(\text{GeV}^{2}/c^{4})$ for (2).

3.2.4 CKM unitarity

The determination of the CKM matrix elements $(V_{CKM}^{exp})_{ab}^{i}$ is usually done under the implicit assumption of a zero active-sterile mixing matrix, $\Theta=0$. The measurements of the $(V_{CKM}^{exp})_{ab}^{i}$ therefore need to be adjusted to take into account effects of RHNs.

Firstly, the smallest element of the CKM matrix, $(V_{CKM})_{ub}$, can be neglected in our study as its absolute value $|(V_{CKM})_{ub}|^{2}\sim 10^{-5}$ is much smaller than our sensitivity to the $\Theta$ parameter. Hence, under the assumption of the unitary of the CKM matrix, one can derive the following relation:

[TABLE]

Thus, we use the various experimental measurements of $(V_{CKM}^{exp})_{us}$ Antonelli:2010yf ; Follana:2007uv ; Amhis:2012bh and $(V_{CKM}^{exp})_{ud}$ Patrignani:2016xqp to simultaneously constrain the true value of $|(V_{CKM})_{us}|$ and active-sterile mixing matrix $\Theta$.

Following Refs. Antusch:2014woa ; Drewes:2015iva , the experimental measurements and true value of CKM matrix element $(V_{CKM})_{us,ud}$ are related via

[TABLE]

where we defined the functions $f^{i}$ to encode the contribution of RHNs to the process considered in each experiment. The decay processes considered to extract the value of $|(V_{CKM}^{exp})_{us}|$, and the $f(\Theta)$ functions, are given by Antusch:2014woa

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The situation is simpler in the determination of the $|(V_{CKM}^{exp})_{ud}|$ element as the uncertainty is dominated by the superallowed $0^{+}\to 0^{+}$ nuclear beta transitions measurements, which need to be modified accordingly to:

[TABLE]

The experimentally measured values of $|(V_{CKM}^{exp})^{i}_{us}|$ in each of the decay processes above are listed in Tab. 4, and the value of $|(V_{CKM}^{exp})_{ud}|=0.97417\pm 0.00021$ is taken from the world average Patrignani:2016xqp .

We thus construct the likelihood for this constraint from a chi-squared function, $2\ln\mathcal{L}=-\chi^{2}$, where the discriminant measures the deviation of the true value $(V_{CKM})_{us,ud}$ and the experimental measurements $(V_{CKM}^{exp})^{i}_{us,ud}$, and is given by

[TABLE]

Due to the unitarity relation in Eq. 56, the value $(V_{CKM})_{ud}$ is obtained from $(V_{CKM})_{us}$ for every parameter point, and thus the only free floating parameters are the value of $(V_{CKM})_{us}$ and the active-sterile mixing matrix, $\Theta$. For simplicity, and since this is the only constraint to depend strongly on the value of $(V_{CKM})_{us}$, we optimise on its value for each $\Theta$, which removes the necessity of making $(V_{CKM})_{us}$ part of the scanning model. This approach is similar to the discussion in Drewes:2015iva , but we improve upon it by optimising on the true value $(V_{CKM})_{us}$, including the $\Theta$ corrections, for each parameter point, rather than the value measured experimentally.

3.2.5 Neutrinoless double-beta decay

Double-beta decay refers to the decay of two neutrons into two protons while emitting two electrons and two anti-neutrinos. In case of neutrinos having a Majorana nature, lepton number would be violated and neutrinoless double-beta decay ($0\nu\beta\beta$) induced. Besides the exchange the light neutrinos, the exchange of RHNs is similarly possible and would alter the expected effective neutrino mass $m_{\beta\beta}$. The effective mass is constrained by half life measurements of $0\nu\beta\beta$ decay. The most stringent limits are currently set by the GERDA experiment (Germanium) Agostini:2017iyd with $m_{\beta\beta}<0.15-0.33\;\text{eV}$ (90% CL), and KamLAND-Zen (Xenon) KamLAND-Zen:2016pfg , $m_{\beta\beta}<0.061-0.165\;\text{eV}$ (90% CL). The effective mass $m_{\beta\beta}$, can be theoretically evaluated in term of the mixings and masses of the light and right handed neutrinos Drewes:2016lqo

[TABLE]

Hereby, the first term denotes the contribution from LHNs, the second the one from RHNs. With a typical momentum exchange of around 100 MeV in $0\nu\beta\beta$ decay, RHNs with a mass above this threshold participate in the process only virtually. This suppression is taken into account by the factor Drewes:2016lqo

[TABLE]

The typical momentum exchange $p^{2}$ depends not only on the specific isotope in consideration but is also subject to the theoretical model in which the constraints are derived and the value of the nucleon axial-vector constant. An overview is given in Faessler:2014kka : For our analysis, we use the “Argonne” model and the lower of the two values for $p^{2}$ (quenched), which yields the most conservative constraints: $\sqrt{\langle p^{2}\rangle}=178\;\text{MeV}$ for xenon, and $\sqrt{\langle p^{2}\rangle}=159\;\text{MeV}$ for germanium. A more dedicated analysis of the impact of different limits due to nuclear uncertainties is beyond the scope of this work. Since we are focusing on profile likelihood for our results, this approach is largely equivalent to profiling over systematic uncertainties assuming a flat prior that spans the entire range of values $\langle p^{2}\rangle$ in ref. Faessler:2014kka . For our analysis we use the experimental values, as stated above, as one-sided Gaussian likelihoods, choosing the higher of the two values in order to remain conservative.

3.2.6 Big Bang Nucleosynthesis

If RHNs decay shortly before or during BBN, the typical energy of decay products, here $\sim M_{I}\geq 50\rm\;MeV$, is significantly higher than the plasma temperature at that time, $\sim 100$ keV. Therefore, either by dissociating formed nuclei, or by causing deviations from thermal equilibrium, they will affect the abundances of primordial elements, which are however observationally well constrained. The requirement that the RHN decay happens sufficiently early enough before BBN implies an upper limit on the lifetime ($\tau_{I}$) of RHNs, or equivalently, a lower bound on the mixing $U_{I}^{2}$ Dolgov:2003sg . However, in the presence of multiple RHN species, BBN cannot constrain individual mixing angles $U_{\alpha I}^{2}$ (22) but only the total mixing $U_{I}^{2}~{}\eqref{UI}$.

We consider leptonic decay channels for all RHNs masses, when kinematically allowed, as well as hadronic decays to mesons and leptons. As shown in Bondarenko:2018ptm , for low masses the hadronic decay width is dominated by channels with a single meson and a lepton, while for masses above the hadronisation scale, $\Lambda_{\textrm{had}}\sim 1$ GeV, it can be approximated by computing the decay to free quarks. The decay width for each topology is listed in Appendix B.2, with expressions and values for the decay constants taken from  Atre:2009rg , Gorbunov:2007ak , Canetti:2012kh , Bondarenko:2018ptm and Ballett:2019bgd , along with a detailed comparison of the various expressions.

In the current study, we require the lifetime of each RHN to be less than $0.1s$ Ruchayskiy:2012si , which is implemented in the likelihood as a step function. In principle, this limit can be weakened if the lightest active neutrino has a mass $<\mathcal{O}(10^{-3})$ eV, since the RHNs do not necessarily thermalize in this case Hernandez:2014fha . We leave, however, the implementation of refined BBN constraints in GAMBIT for future work. Note that a lifetime bound that is stronger by a factor of two would lead to proportionally stronger constraints on the total mixing $U_{I}^{2}$.

3.3 Direct RHN searches

Different experiments search with various approaches directly for RHNs. One can distinguish between three types: peak searches (PIENU), searches at beam dump experiments (PS-191, CHARM, E949, NuTeV), and searches at $e^{+}e^{-}$ or $pp$ colliders (DELPHI, ATLAS, CMS).

One possibility to look for RHN, is to search for peaks in the lepton energy spectrum of a meson decay. If, for example, a meson of mass $m_{X}$ decays into an RHN of mass $M_{I}$ and an electron/muon with mass $m_{l_{\alpha}}$, this peak will be approximately at

[TABLE]

Even in situations where backgrounds are sizeable, a peak search can hence be used to impose constraints on the mixing.

In beam dump experiments, the large background signal that is usually present near the target hinders the detection of charged particles that are produced along with the RHNs. On the other hand, RHNs with mass below the D meson scale can be long-lived enough to travel macroscopic distances. Looking for their charged decay products some distance away from the target leads to (almost) background-free experimental situations.

In collision experiments ($e^{+}e^{-}$ or $pp$), vector bosons or mesons get produced that subsequently can decay leptonically. The bounds on these processes are then able to constrain the corresponding active-sterile mixing angles in a certain mass range.

To implement the direct detection constraints as likelihoods, we follow two different approaches, depending on the information that is provided in each study. Firstly, some of the experiments found no signal events and had no background counts after cuts (DELPHI, CHARM, PS191 and NuTeV). In this case, since the processes in the experiments are essentially Poissonian, we construct the likelihood (to observe $n$ events) as a Poisson distribution. The number of expected counts, $\mu$, is a function of the RHN masses and mixings, i.e. $\mu=\mu(M_{I},U_{{\alpha}I}^{4})$ (assuming the experiment does so as well, the fourth power takes both production and decay of RHNs into account). For expected $\mu$ events and background $b$, the likelihood is:

[TABLE]

With no reported detections ($n=0$) and background cuts reducing $b$ to approximately zero,

[TABLE]

To connect $\mu$ with our model parameters, we use the fact that the expected signal counts are proportional to the LHN-RHN mixing, $\mu\propto U_{{\alpha}I}^{4}$. The factor of proportionality is set to reproduce the results from the experimental papers (assuming that these limits are based on the common Feldman-Cousins procedure Feldman:1997qc , where e.g. a 95% CL upper limit would correspond to an expected number of signal counts of $\mu=3.09$).

On the other hand, for the experiments which either quote non-zero signal events and/or backgrounds, or if this information is ambiguous (CHARM ($\nu_{\tau}$ re-interpretation), PIENU, ATLAS and E949), we model the constraint likelihood as Gaussian upper limits, i.e. we model them as half-Gaussians with zero mean and error set according to the confidence level at which the results are presented. For example, in the case of an experiment that presents limits at $90\%$ CL, for a half Gaussian, this lies within $1.28\sigma$ of the mean.

It is worth noting that collider experiments often use simplified model assumptions to compute the confidence level intervals presented in their results. Since we use these to construct our likelihoods, we are incorporating these assumptions as well, in spite of the fact that our confidence intervals are computed by profiling over the multidimensional parameter space. Given that a full collider simulation is beyond the scope of this study, we employ the provided simplified model limits as given. We acknowledge, however, that the true limits may be slightly weaker due to, e.g a reduction of the production cross-section, and we defer the exploration of the differences between the collider predictions of simplified and full models to future work.

3.3.1 PIENU

The PIENU experiment PIENU:2011aa sought to detect RHNs in the mass range of $68-129$ MeV by searching for peaks in the energy spectrum of the decay process $\pi^{+}\rightarrow e^{+}\nu$. It was, hence, sensitive to the mixing $|\Theta_{eI}|^{2}\equiv U_{eI}^{2}$ and $\mu$ in eq. (72) is also taken to scale as $U_{eI}^{2}$ in our analysis. Although no peaks were found, exact information on the number of background events is unavailable. Further, production processes in peak searches are, in general, unaffected by the Majorana/Dirac nature of the RHNs; hence, no correction is necessary here.

The constraints on $U_{eI}^{2}$ are at 90% CL, so it is implemented in GAMBIT as a half-Gaussian with zero mean and error set at 1.28$\sigma$.

After our analysis was complete we became aware of the slightly stronger updated constraints presented in Ref. Aguilar-Arevalo:2017vlf , which are not included in our scan.

3.3.2 PS-191

This experiment Bernardi:1987ek was designed for the purpose of detecting neutrino decays. RHNs would be produced via either of the following mechanisms: $\pi^{+}/K^{+}\rightarrow e^{+}\nu_{e}$, or $\pi^{+}/K^{+}\rightarrow\mu^{+}\nu_{\mu}$, and would then decay via $\nu_{R}\rightarrow\mu^{-}e^{+}\nu$, $\nu_{R}\rightarrow e^{-}\mu^{+}\nu$, $\nu_{R}\rightarrow e^{-}\pi^{+}$, $\nu_{R}\rightarrow\mu^{-}\mu^{+}\nu$, $\nu_{R}\rightarrow\mu^{-}\pi^{+}$ or $\nu_{R}\rightarrow e^{-}\pi^{+}\pi^{0}$. Thus, PS-191 could constrain the quantities $U_{eI}^{4}$ and $U_{\mu I}^{4}$ for RHNs with a mass between $20-450$ MeV.

Having found no signal or background events, it placed constraints on these quantities at 90% CL. We deviate from the original analysis in two ways. The first is necessitated by the fact that in the original analysis, the constraints were derived under the assumption that the RHNs interact only through the charged current. In Ruchayskiy:2011aa , these limits were re-interpreted with the inclusion of neutral current interactions. Thus, instead of the signal count being proportional to the fourth power of the relevant flavour mixing, it is proportional to $U_{e/\mu I}^{2}\times\sum_{\alpha}c_{\alpha}U_{\alpha I}^{2}$, with the coefficients given by

[TABLE]

We use these revised bounds here. The limits are encoded in likelihood form as in eqn. (72), with the aforementioned proportionality factor being $2.44$.

3.3.3 CHARM

RHNs were searched for in CHARM Bergsma:1985is using two strategies, one with a neutrino beam from dumping protons on copper (BD) and another using a wide-band neutrino beam (WBB) from primary protons.

In BD, the production of RHNs was assumed to occur through the decay of D mesons. They would then decay via $\nu_{R}\rightarrow e^{+}e^{-}\nu_{e}$, $\nu_{R}\rightarrow\mu^{+}\mu^{-}\nu_{\mu}$ or $\nu_{R}\rightarrow e^{+}\mu^{-}\nu_{e}$, $\mu^{+}e^{-}\nu_{\mu}$ (and the anti-particle counterparts) and the decay products were looked for.

In WBB, RHN production was assumed to occur via neutrino-nucleus neutral current scattering $\nu_{\mu}N\rightarrow\nu_{R}X$. The subsequent decay $\nu_{R}\rightarrow\mu R$, R representing hadrons, was then searched for. The limits from the WBB analysis are, however, weaker than those exerted by other experiments in the same mass range, and are not considered here.

The BD analysis yielded no candidate events or background and hence placed limits on $U_{eI}$ and $U_{\mu I}$ at 90% CL. Further, the original analysis assumed the possibility of RHNs interacting solely via the charged current; we use the results re-interpreted after the inclusion of neutral current interactions Ruchayskiy:2011aa as discussed in section 3.3.2, i.e. the signal count is proportional to $U_{e/\mu I}^{2}\times\sum_{\alpha}c_{\alpha}U_{\alpha I}^{2}$ and once again use eqn. (72) to represent the likelihood, with the proportionality factor being $2.44$.

In Orloff:2002de , the data from the CHARM experiment was re-analyzed assuming that RHNs mix solely with tau-flavoured leptons, and was able to place limits at 90% CL on $U_{\tau I}$, which we implement as a half-Gaussian with zero mean and error set at 1.28$\sigma$.

Dirac RHNs were assumed in both the original and tau-specific analyses, so the limits presented are also re-scaled by dividing them by $\sqrt{2}$.

3.3.4 E949

In this experiment Shaykhiev:2011zz ; Artamonov:2014urb , RHNs were searched for in the decay of kaons produced in a beam dump: $K^{+}\rightarrow\mu^{+}\nu_{R}$. Constraints on $U_{\mu I}$ were placed at $90\%$ CL in the mass range $175-300$ MeV; we also divide the limits by a factor of $\sqrt{2}$ to account for the Majorana nature of RHNs in our model.

The likelihood is modeled as a half-Gaussian with zero mean, error set at 1.28$\sigma$ and $\mu\propto U_{\mu I}^{2}$.

3.3.5 NuTeV

The NuTeV experiment Vaitaitis:1999wq searched for RHNs through their decay into the following final states: $\mu e\nu$, $\mu\mu\nu$, $\mu\pi$ and $\mu\rho$. They were assumed to be produced in the decay of mesons. $90\%$ CL limits on $U_{\mu I}$ were placed for RHNs with a mass between $0.25-2$ GeV.

Information about the assumed Dirac or Majorana nature of the RHNs is not present, so we take the conservative route and presume Majorana RHNs were considered in the analysis. No candidate events or background were detected, so the likelihood is modeled as in eq. (72), with a proportionality factor of $2.44$ and $\mu$ scaling as $U_{\alpha I}^{4}$.

3.3.6 DELPHI

At DELPHI Abreu:1997uq , $e^{+}e^{-}\rightarrow Z^{0}\rightarrow\nu_{R}\bar{\nu}$ was the dominant RHN production mechanism; the process $Z^{0}\rightarrow\nu_{R}\bar{\nu_{R}}$ would be suppressed due to the additional $U^{2}$ factor. The products of the RHN decaying via the weak and neutral current were then searched for, according to: $\nu_{R}\rightarrow\nu Z^{\ast}$, $Z^{\ast}\rightarrow\nu\bar{\nu}$, $l\bar{l}$, $q\bar{q}$ or $\nu_{R}\rightarrow l^{\prime}W^{\ast}$, $W^{\ast}\rightarrow\nu\bar{l},q\bar{q^{\prime}}$. DELPHI could constrain $\Theta_{eI}$, $\Theta_{\mu I}$ and $\Theta_{\tau I}$ for RHNs having a mass between $0.5-80$ GeV.

Since the RHNs could have existed long enough to travel macroscopic distances of upto $100$ cm, different signatures had to be considered and the analysis was split to tackle the short- and long-lived cases separately.

In the short-lived RHN case, depending on the particle mass, two signatures were looked for. For masses less than about 30 GeV, due to the large boost received by the RHNs, the signature would be a monojet. Background coming from leptonic Z boson decays or $\gamma\gamma$ processes were accounted for. Higher masses open the decay channel into $q\bar{q}$ (and a lepton, depending on the channel), and the signature in this case would be two acollinear jets which are also acoplanar with respect to the beam axis. Most of the background in this scenario came from hadronic Z decays with missing energy; a neural network was used to remove all of them from the final data.

Longer-lived RHNs were looked for using displaced vertices and calorimeter clusters. The former was useful in tracking RHNs with an intermediate lifetime; however, a cluster finding algorithm along with vertex reconstruction did not find any signals. Calorimeter clusters were used to detect the longest-lived RHNs, whose decay products would interact with the outermost layers/components of the experimental setup: the signature would be a cluster of hits in a small angular region coincident with the beam collision, which could be traced back to the initial interaction point.

The analysis was carried out assuming Majorana RHNs and yielded one candidate event and no background events. In our analysis, this means the proportionality factor is $3.09$ and $\mu$ scales as $U_{\alpha I}^{4}$.

A caveat must be mentioned here: the DELPHI analysis presented bounds on the mixing in a flavour-independent manner: the limit on $U^{2}$, as presented in the paper, applies equally to $U_{e}^{2}$, $U_{\mu}^{2}$ and $U_{\tau}^{2}$, as they mention. In the mass range under consideration, the mass of the tauon will, of course, influence the strength of the limit and, as they quote, the presented bounds become weaker for masses below $\sim 4$ GeV. However, the extent of the kinematic suppression due to the tauon mass is not quantitatively discussed; we use the limits as is, noting that it is highly likely that NA62 will subsume these bounds in the near future Lurkin:2017tmu .

3.3.7 ATLAS

The process relevant for RHN production in ATLAS Aad:2015xaa is $pp\rightarrow(W^{\pm})^{\ast}\rightarrow l^{\pm}\nu_{R}$. The RHNs were taken to be heavier than the W boson, allowing it to decay to a lepton a W boson: $\nu_{R}\rightarrow l^{\pm}W^{\mp}$; the W boson would then decay predominantly into a quark-antiquark pair, and the signature of this decay chain was searched for, with either two electrons or muons in the final state.999 There is an ongoing dispute in the literature on whether the rate of LNV processes at collider experiments are always suppressed by the small parameters $\epsilon_{i}$ and $\mu$ in Eq. (2.5) and therefore unobservably small (roughly of the order of the ”naive seesaw estimate”) Kersten:2007vk ; Moffat:2017feq or whether coherent flavour oscillations can lead to LNV signatures in spite of the smallness of these parameters Anamiati:2016uxp ; Antusch:2017ebe ; Antusch:2017pkq . In the range of $M_{I}$ below the electroweak scale under consideration here, the strongest direct search constraints do not come from experimental signatures that rely on LNV, and our results are therefore only mildly affected by the outcome of this discussion. Hence, in our analysis, $\mu\propto U_{\alpha I}^{4}$, $\alpha=e,\mu$. The original analysis was carried out under the assumption of Majorana RHNs, so no additional correction is necessary.

The analysis placed $95\%$ CL limits on the two mixing angles in the mass range of $100-500$ GeV. Details on the number of observed/expected events and background is available and could be cast into a likelihood function combining Poissonian and Gaussian errors; however, we find that implementing the limits in GAMBIT as a half-Gaussian with zero mean and error set at 1.64$\sigma$ reproduces the experimental limits well enough for the purpose of a global fit.

3.3.8 CMS

With the LHC having run with a center-of-mass energy of 13 TeV, the CMS detector searched for different event signatures of the same process as ATLAS. $95\%$ CL limits were calculated for $U_{eI}$ and $U_{\mu I}$ for RHNs with mass between 1 GeV and 1.2 TeV Sirunyan:2018mtv .

As before, Majorana RHNs were assumed in the analysis, and our implementation of the limits mirrors that of ATLAS.

Note that updated bounds from ATLAS Aad:2019kiz and CMS CMS:2018szz ; Sirunyan:2018xiv have been released, but are not included, since these papers came out after our scans were completed. However, the new bounds from ATLAS are comparable to those from DELPHI, and the newer dilepton search from CMS only produces stronger bounds for RHN masses above $\sim 500$ GeV, which is beyond our range of study.

3.3.9 LHCb

LHCb has performed direct searches for heavy neutrinos. The most recent results Aaij:2014aba were derived with an inconsistent model and have been corrected in Ref. Shuve:2016muy . They are subdominant in the mass range considered here. In Ref. Antusch:2017hhu the results of a generic long lived particle search Aaij:2016xmb has been re-interpreted in the context of heavy neutrinos. We do not include these results here because the conservative interpretation does not yield stronger bounds than the ones we include.

3.3.10 Other experiments

Further measurements at Borexino Back:2003ae , Bugey Hagner:1995bn , SIN Abela:1981nf , BEBC CooperSarkar:1985nh , JINR Baranov:1992vq , TRIUMF Britton:1992pg ; Britton:1992xv , OKA Sadovsky:2017qsr ; Aguilar-Arevalo:2019owf , ISTRA Duk:2011yv , NOMAD Astier:2001ck , NA62 CortinaGil:2017mqf , Belle Liventsev:2013zz , KEK Hayano:1982wu ; Yamazaki:1984sj and T2K Abe:2019kgx have both published constraints on RHNs. We do not indculde them here because, with the present data, they are subdominant or cover a different mass range.

4 Scanning strategy and parameter ranges

In this work, we focus on the exploration of the RHN parameter space using frequentist statistics. Our main goal is to establish the ranges of RHN parameters that are not yet explored by experiments, and a frequentist approach delivers a suitable and prior-independent method. We are dealing with a high dimensional parameter space, which we have to project into two-dimensional plots. To this end, the central quantity of interest is the profile likelihood,

[TABLE]

which is, for fixed parameters of interest $\theta_{1}$ and $\theta_{2}$, the maximum value of the (log-)likelihood function that can be obtained when maximizing over the remaining parameters $\vec{\eta}$.

We emphasize that the main goal of this work is to establish conservative constraints on RHN mixings and masses by profiling over all relevant parameters. We do not perform a proper goodness-of-fit analysis to experimental data, which would require sampling of experimental results; given the large range of included experimental results and the sometimes limited knowledge about individual experiments this is beyond the scope of the current work. Instead, likelihoods are included in a approximate fashion that allows to reproduce published experimental results, and we use Wilks’ theorem Wilks:1938dza to approximate the sampling statistics of log likelihood ratios and estimate confidence contours when necessary.

Our scanning strategy is designed in order to explore the complex parameter space of the RHN model such that we obtain reliable results for the projections shown in this work. To this end, we perform a large set of scans with different settings which we then merge into a single dataset. We study the normal (NH) and inverted (IH) hierarchy independently, in order to avoid artificially favouring one over the other due to the different normalisation of the active neutrino likelihoods (c.f. Section 3.1). Hence, we make independent scans for each of the neutrino mass hierarchies, normal and inverted, for the full set of scans described below.

4.1 Parameters and priors

The parameter ranges and priors for the original scans can be seen in Table 5. We emphasize that ‘priors’ do here not correspond to priors in the Bayesian sense, but rather determine the efficiency with which different regions of the parameter space are explored. For convergent scans, the results are prior-independent. We have chosen to split the complex angles $\omega_{ij}$ into their real and imaginary parts. The active-sterile mixings depend strongly on the imaginary parts of $\omega_{ij}$ $\left(\Theta^{2}\sim\frac{\exp(2\text{Im}(\omega))}{M}\right)$ and large values of Im$\omega$ produce mixings that are too large to pass any constraints, so we take a conservative range Im$\omega\in[-15,15]$, and also pre-emptively disallow choices that lead to $|\Theta|_{ij}^{2}>1$. As discussed in 2.5, a condition for an approximate $B-\bar{L}$ symmetry to be realized is for two RHNs to have almost degenerate masses, which extends the range of the mixings so that they can be probed by experiments. This provides motivation for using a logarithmic prior on the RHN masses, also allowing the scanner to sample better the region close to the limits of the most constraining experiments/observables.

The C-I parametrisation, as defined in Section 2.4, together with the particular parametrisation choice of $R$ in eq. 30, was found to not fully cover the entire parameter space. To circumvent this and ensure that all possible couplings are covered by the scans, we introduce an additional parameter to the scan $R_{\rm{order}}$ with discrete values $[1,6]$ corresponding to each of the possible permutations of the definition of $R$ in terms of $R^{ij}$. This allows full coverage of the coupling space and, since the likelihood is conceptually independent of the order in $R$ (and confirmed by the data), it ensures an uniform distribution of values in the parameter $R_{\rm{order}}$.

Out of the active neutrino parameters, only $\alpha_{1}$ and $\alpha_{2}$ are unconstrained by oscillation data, hence they are allowed to vary freely from 0 to $2\pi$ with flat priors. The ranges for the other neutrino phases and angles are taken as the widest of the $3\sigma$ ranges, for normal or inverted hierarchy, from the NuFit collaboration NuFit , also with flat priors. The mass of the lightest active neutrino, $m_{\nu_{0}}$, has a definite impact on the lower bound of $U_{I}^{2}$ (23) Drewes:2015iva , so we choose a logarithmic prior, which enables us to examine this impact in greater detail than a flat prior would allow and keeps the BBN limits relevant Hernandez:2014fha . The upper limit on $m_{\nu_{0}}$ is chosen as the broad cosmological bound given by Planck Ade:2015xua , $\sum m_{\nu}<0.23$ eV101010This upper limit is not very conservative in light of Planck data, a more conservative bound would be $\sum m_{\mu}<0.6$ eV Abe:2018emu . However, there is no effect of this constraint on our data as most high likelihood data points lie in the limit $m_{\nu_{0}}\to 0$. We have, nevertheless, studied a subset of cases with the conservative bound and indeed found them to not be relevant.. In order to better fit the active neutrino data, the mass splittings $\Delta m^{2}_{21}$ and $\Delta m^{2}_{3l}$ are chosen as scan parameters, where $l=1$ and $\Delta m^{2}_{3l}>0$ for normal hierarchy and $l=2$ and $\Delta m^{2}_{3l}<0$ for inverted hierarchy.

Since the construction of the mixing matrix in the C-I parametrisation depends on $m_{H}$ (1-loop correction), as seen in 2.4, we take $m_{H}$ as a nuisance parameter with a Gaussian distribution around its averaged measured value PDG17 and a flat prior. Other SM parameters are fixed to their PDG values PDG17 .

4.2 Targeted scans

We encountered a number of challenges while sampling the full RHN parameter space. One reason is connected to the behaviour of the likelihood function over the whole parameter range. The adopted scanning algorithm (Diver, see below for details) is designed to find regions of maximum likelihood across the parameter space. However, as we will discuss later when we study the effect of each individual observable, most constraints have flat contributions to the likelihood in a large portion of the parameter space. Hence, the scanner often does not fully explore large regions with equal or worse likelihood. This happens especially near the experimental bounds. Furthermore, although high couplings are possible between active and sterile neutrino sector, they often lie in the symmetry protected regime, as described in Section 2.5 and/or require severe fine-tuning of the parameters. Again, exploring these regions turned out to be challenging.

Therefore, we designed and performed a large set of targeted scans to fully saturate the experimental bounds, the list of which can be found in Table 6. The design strategies we adopted for these targeted scans can be summarised as follows.

First, all targeted scans were performed using a differential RHN model, where the parameter $M_{2}$ is replaced by $\Delta M_{21}$, with a logarithmic prior. This allows the exploration of the symmetry protected region, with near degenerate masses for two right-handed neutrinos.

Most of the experimental bounds occur at high couplings, thus in order to encourage the scanner to explore the high coupling regions, we added an artificial likelihood to the scan to drive the scan to the unexplored boundaries. To saturate the experimental bounds for each coupling $U_{\alpha I}^{2}$, $\alpha=e,\mu,\tau$, different targeted scans were performed using this coupling slide likelihood on each of the couplings, of the form $s\log U_{\alpha I}^{2}+m\log M_{I}$. Table 6 shows the parameter that is optimised in each scan, $\alpha$, and the coefficients, $(s,m)$. This contribution was later removed from the data in the postprocessing stage.

The targeted scans were further split along the $M_{I}$ axis following the limits of the various experimental constraints (mostly from direct searches). This ensures that each coupling (with the selection above) saturates the most relevant experimental upper bound in each mass range. Additionally, some scans used different values of $\Delta M_{21}$ and/or $m_{\nu_{0}}$ to further force the scan into fine-tuned regions of parameter space. The ranges used for $M_{I}$, $\Delta M_{21}$ and $m_{\nu_{0}}$ for each scan are specified in Table 6.

A similar strategy was used to saturate the BBN bound at low couplings. Three scans were performed for each hierarchy, with slide coefficients $(s,m)=(-0.5,-0.5)$ on each coupling $U_{\alpha I}^{2}$, $\alpha=e,\mu,\tau$ . To further optimise on low couplings, these scans were performed for fixed $m_{\nu_{0}}=10^{-10}$ and a narrow range on Im$\omega$ $\in[-0.5,0.5]$. With these settings the BBN bound was fully saturated in the explored mass range.

We found that some of the experimental likelihoods provide positive contributions to the total likelihood in specific regions of the parameter space. This forced the scan towards those regions, leaving others unexplored. Although this is a rather interesting feature, and will be discussed in detail later, it prevented a thorough exploration of the full parameter space. We thus chose to remove the likelihood contribution of $R^{K}_{e\mu}$ from the total likelihood that drives the scan, adding it later in postprocessing. Other likelihoods with positive contributions, $\Gamma_{\rm{inv}}$, CKM and $R^{\tau}_{e\mu}$, tended to force the scan towards large $U_{\tau I}^{2}$ couplings. Although desirable to saturate the limits, this also left regions with low $\tau$ coupling undersampled. Thus, a cut on the coupling $U_{\tau I}^{2}$ was enforced in some scans to fully sample all regions.

The adopted strategy for scanning was driven by the need to fully sample the parameter space. The results from all the diverse scans were combined into a single dataset after some postprocessing (see below). This does not pose a problem for the statistical interpretation, since we are interested in the profile likelihood, which only becomes more accurately estimated when adding additional chains.

4.3 Scanning framework

To perform the detailed scans, we make use of the GAMBIT framework, as described in Appendix A, and the differential evolution scanner Diver, version 1.0.4 ScannerBit , which is a self-adaptive sampler, capable of sampling the profile likelihood more efficiently than other scanners. We choose a population size of NP = 19200 and a convergence threshold of convthresh = $10^{-10}$. After some tests, we have concluded that the aggressive $\lambda$jDE setting in Diver provides an improvement on the sampling of the parameter space, since it is more suited for sampling fine-tuned regions.

These scanner settings, including the very low convergence threshold, together with the scanning strategy described above, ensure a thorough exploration of the parameter space, albeit at the price of CPU time. Despite the fact that none of the observables used required heavy computation or simulations, most scans took between 2 and 10 hours of running time on a large number of supercomputer cores varying between 250 and 780. All tests and scans were carried out across several supercomputer facilities, including the MareNostrum supercluster in Barcelona, Marconi in Bologna, LISA/Surfsara through the University of Amsterdam and Prometheus in Krakow.

4.4 Data postprocessing

Upon completion of the scans, a number of postprocessing tasks were performed on the data to prepare it for plotting. As previously mentioned, the first of these tasks was to remove the artificial coupling slide likelihood used to drive the scans to high couplings.

Due to the large amount of scans performed and the low convergence threshold used, the size of the samples surpassed 1TB for each hierarchy, rendering them unmanageable for most plotting routines. We hence performed a few operations on the scan results prior to combining them. With the target of showing profile likelihood plots in the $M_{I}$ vs $U_{\alpha I}^{2}$ planes, we hence extracted a subset of the data points optimised in these planes, with a resolution of $10^{-5}$. Since most scans were targeted to saturate the limits for a particular coupling (see Table 6) we perform this reduction of the data in the respective mass vs coupling two-dimensional planes. The combined set will hence be optimised for all couplings. Additionally, and independent reduction of the data is performed on the planes $m_{\nu_{0}}$ vs $U_{\alpha I}^{2}$, since we intent to study the effect of $m_{\nu_{0}}$ cuts on the coupling limits.

The flavour label of the heavy neutrinos is arbitrary, and the experimental constraints on a heavy neutrino with a given mass cannot depend on the labelling. However, for reasons explained in more detail in appendix C, the scanning strategy outlined in Sec. 4.2 introduces a bias that suggests that the constraints differ for $N_{1}$, $N_{2}$ and $N_{3}$. Hence, to remove this bias in the labels, after combining the reduced datasets for all the scans, we conduct a symmetrization procedure over the combined datasets. We therefore symmetrize over $M_{I}$ as well as $U_{\alpha I}$, which will increase the size of the datasets six fold.

Lastly, in order to compare with the $n=2$ case, two further datasets were obtained, for normal and inverted ordering, where the data points are required to lie in the symmetry protected region.

Out of the incalculable amount of data points we collected through our scanning procedures, a total of 40.7 million valid data samples were used for plotting. Of which 11M correspond to normal hierarchy and 10M for inverted hierarchy, optimised on $M_{I}$ vs $U_{\alpha I}^{2}$ planes, and 9.9M for normal and 9.7M for inverted hierarchy, optimised on $m_{\nu_{0}}$ vs $U_{\alpha I}^{2}$ planes. The datasets with points in the symmetry protected region have over 71k and 20k valid data samples for normal and inverted hierarchy, respectively. These samples can be found in Zenodo Zenodo_RHN .

4.5 Capped likelihood

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

5 Results and discussion

5.1 General constraints on the RHN mass and mixing

The constraints are shown in Figs. 1-4 for the couplings $U_{\alpha I}^{2}$ to the active neutrino flavours $\alpha=(e,\mu,\tau)$, as well as their combination $U_{I}^{2}=\sum_{\alpha}U_{\alpha I}^{2}$, as functions of the heavy neutrino masses $M_{I}$. Here, the second index can refer to any of the heavy neutrino flavours $I=(1,2,3)$, because their labelling is not physical. Figs. 5-7 show the combinations of couplings $U_{\alpha I}U_{\beta I}$ with $\alpha\neq\beta$. The allowed profile likelihood regions are flat for most of the parameter space, in particular for small couplings $U_{\alpha I}^{2}$, and drop smoothly at high couplings following the relevant upper limits. The white lines around the experimental limits mark the 1$\sigma$ and 2$\sigma$ contours, which are estimated assuming Wilks’ theorem with 2 degrees of freedom111111All profile likelihood plots were created using pippi pippi ..

The largest values of mixings $U_{\alpha I}^{2}$ and $U_{\alpha I}U_{\beta I}$ for all flavours are allowed for $M_{I}$ above the masses of the weak gauge bosons. In this regime the direct searches at colliders are sub-dominant, and the heavy neutrino properties are primarily constrained from above due to electroweak precision observables, lepton flavour violation and CKM constraints. The upper limits on the couplings $U_{\alpha I}^{2}$ and $U_{\alpha I}U_{\beta I}$ within $2\sigma$ of the highest likelihood for each hiearchy and flavour in the high mass region can be found in Table 7. It can be readily noticed that the upper limits for the $\tau$ couplings is much larger than for the other two flavours, which can be understood because the limits from EWPO and LFV are stronger for $e$ and $\mu$ (see also Sec. 5.4). In particular, the combination $U_{eI}U_{\mu I}$ has the smallest of upper limits, as shown as well in Fig. 5, due to strong constraints from LFV observables, specifically $\mu\to e\gamma$ and $\mu-e$ conversion (see Figs. 34 and 35 in Appendix  D).

For $M_{I}$ between the masses of the $D$ mesons and the $W$ boson the limits from direct searches dominate because the heavy neutrinos can be produced efficiently via the $s$-channel exchange of on-shell $W$ bosons. In the range between the $D$ meson masses and the $W$ boson mass, the limits from the DELPHI Abreu:1997uq and CMS Sirunyan:2018mtv experiments compete to impose the strongest bound.

Below the $D$ meson mass the constraints on $U_{eI}^{2}$ and $U_{\mu I}^{2}$ are dominated by direct search constraints from fixed target experiments, in particular CHARM Bergsma:1985is and NuTeV Vaitaitis:1999wq above the kaon mass, PS-191 Bernardi:1987ek and E949 Artamonov:2014urb between the pion and kaon mass and pion decay experiments at even lower masses. In this regime the global constraints on $U_{eI}^{2}$ and $U_{\mu I}^{2}$ are in good approximation given by the direct search constraints, as discussed in Sec. 5.2 and Figs. 8-10. This is in contrast to the model with $n=2$, where the global fits rule out a significant mass range below the kaon mass that appears to be allowed if one simply superimposes the direct constraints in the mass-mixing planes Drewes:2016jae . For $U_{\tau I}^{2}$, the direct search constraints are much weaker, the limit from long-lived particle searches by DELPHI remains the most significant one in our scans. Figure 3 shows that direct searches become subdominant for the $\tau$ coupling and the EWPO limit is saturated for a considerable range of masses below the kaon mass.

For masses below roughly $0.3$ GeV the global constraints are stronger than the sum of their ingredients due to an interplay of the lower bound from BBN on the mixings, the upper bounds on $U_{eI}^{2}$ and $U_{\mu I}^{2}$ from direct searches and the constraints on the heavy neutrino flavour mixing pattern from neutrino oscillation data (discussed further below in Sec. 5.3). The latter disfavours large hierarchies amongst the couplings to individual SM flavours, though these constraints are weaker than in the model with $n=2$ Drewes:2016jae ; Drewes:2018gkc . This implies that upper bounds on combinations of $U_{eI}^{2}$ and $U_{\mu I}^{2}$ indirectly constrain $U_{\tau I}^{2}$. The BBN constraint on the lifetime does not impose a constraint on any individual coupling $U_{\alpha I}^{2}$, but requires at least some of them to be sizeable and practically translates into a lower bound on $U_{I}^{2}$ that is visible in Figure 4. Both, the BBN constraint and the constraint on the flavour mixing pattern (that will be discussed in more detail in Sec. 5.3 and is visible in Fig. 11) leads to the lower and upper bounds on $U_{\tau I}^{2}$ that are visible in Figure 3.

The upper bound on the total mixing $U_{I}^{2}$ from the global constraints can roughly be identified with the bound on $U_{\tau I}$ across the entire mass range as it is constrained the weakest. The lower bound is again given by the lifetime constraint from BBN. In addition, there is a lower bound from the requirement to explain the light neutrino oscillation data that depends on $m_{\nu_{0}}$ and is therefore only visible if one imposes a cut on this unknown quantity. Our results agree with the analytic estimates made in Ref. Drewes:2019mhg , as will be discussed in Sec. 5.3, and are illustrated in Figure 4.

5.2 Discussion of individual bounds

Figures 8 and 9 show explicitly the effect of direct searches on the upper limits of the $e$ and $\mu$ couplings in the mass range $M_{I}\in[0.1,10]$ GeV. Most of the limits shown are at 90%CL, with the exception of DELPHI at 95%. As expected, they lie between the 1$\sigma$ and $2\sigma$ contours. Some of the experimental limits, PS-191 and CHARM, do not directly constrain an individual coupling, but rather the combination $U_{e/\mu I}^{2}\times\sum_{\alpha}c_{\alpha}U_{\alpha I}^{2}$ (as mentioned in Sections 3.3.2 and 3.3.3), with the coefficients $c_{\alpha}$ from eq. (73). As we profile over the other two couplings, the strongest limit for the $\alpha$ flavour for these experiments would correspond to $(U_{I}^{\rm{exp}})^{2}/c_{\alpha}$, with $(U_{I}^{\rm{exp}})^{2}$ being the reported limit by the experiment. Hence the former ratio is what is shown in the figures as the PS-191 and CHARM limits. As observed, $U_{eI}^{2}$ is constrained by PS-191 and CHARM in the lowest and next-to-lowest mass regions, whereas they are superseded by the limits from E949 and NuTeV for $U_{\mu I}^{2}$. In the lowest mass region for the $\mu$ coupling it would appear that the E949 bound is in fact not saturated as the experimental limit falls below the data. This is however just an artifact of binning and interpolation in that region and the fact that the E949 limit is quite jagged. To illustrate this, we show in Figure 10 a zoom into the lowest mass region from Figure 9, where it can be seen clearly that the profile likelihood follows the limits of E949.

Neutrino oscillation data imposes very strong constraints on the parameter space and disfavours vast volumina in the 18 dimensional model parameter space. In the scenario with $n=2$ this has a visible effect on the projections of the global constraints on the $M_{I}$-$U_{\alpha I}^{2}$ planes Drewes:2016jae , in particular for heavy neutrinos lighter than the kaon, where the interplay between neutrino oscillation data, BBN and direct searches rules our most values of $M_{I}$. This effect strongly depends on the light neutrino mass ordering, and varying the light neutrino oscillation parameters within their experimentally allowed limits leads to visible differences Drewes:2016jae . In the present analysis with $n=3$ the impact of neutrino oscillation data on the likelihoods in the $M_{I}$-$U_{\alpha I}^{2}$ planes is much smaller. This is primarily visible in the third generation and the total mixing, cf. Figs. 3 and 4, where the dependence on the light neutrino mass ordering is weak. The reason is that the larger dimensionality of the parameter space with $n=3$ makes it easier to avoid conflicts with direct or other indirect bounds. With $n=2$ neutrino oscillation data also imposes strong constraints on the flavour mixing pattern Hernandez:2016kel ; Drewes:2016jae ; Drewes:2018gkc . These are also visible in the present analysis, cf. Figs. 11-13, but can be avoided by choosing a sufficiently large value for $m_{\nu_{0}}$. The constraints on the flavour mixing pattern strongly depend on the light neutrino mass ordering, and varying the light neutrino oscillation parameters within their experimentally allowed range has a considerable impact on the predictions.

As mentioned above, EWPO (including $\Gamma_{\rm{inv}}$, $m_{W}$, $W$-decays and $s_{w}$), LFV and CKM constraints become relevant for very large couplings and are thus the dominant limit in the high mass region, as well as a small region at small masses for the $\tau$ coupling (Figure 3). Besides providing constraints, in particular $\Gamma_{\rm{inv}}$ and CKM observables are also responsible for the slight excesses in the total likelihood, which we will discuss in Sec. 5.4. Other constraints included in the analysis have little to no effect on the profile likelihoods as shown above.

Among the leptonic decays, only $R^{K}_{e\mu}$ has some impact on the likelihood, with a negative contribution at masses below $0.45$ GeV. Both $R^{K}_{e\mu}$ and $R^{\tau}_{e\mu}$ show minor excesses in total likelihood, which again will be discussed later. Other lepton universality constraints have only little effect on the likelihood.

Neutrinoless double beta decay sets strong upper bounds on $U_{eI}^{2}$ for generic parameter choices, which strongly disfavours considerable regions of parameter space. However, in the limit where lepton number is approximately preserved the expected signal from $0\nu\beta\beta$ is suppressed. Since many of our parameter points are in this symmetry protected scenario, particularly at high couplings, the impact of this constraint on the likelihoods in the projection on the mass-mixing plane is minimal. This is consistent with what was found in Refs. Drewes:2015iva ; Drewes:2016jae ; Bolton:2019pcu .

The effect of BBN can be seen in the lower limits of Figure 3 and 4. The lower limit on $U_{I}^{2}$ is a direct consequence of BBN, as lower couplings would mean that right-handed neutrinos would not decay before BBN and thus affect the abundance of primordial elements. Although no individual limits are imposed by the BBN constraint on the couplings, the strong upper limit on the $e$ and $\mu$ flavours at low masses has the side effect of setting a lower limits on $U_{\tau I}^{2}$, as seen in Figure 3.

For a more detailed explanation of the effect of each partial likelihood, and associated figures, we refer to Appendix D.

5.3 Lightest neutrino mass and flavour mixing

Oscillation data strongly constrains most of the active neutrino parameters, in particular the mass splittings $\Delta m^{2}_{21}$ and $\Delta m^{2}_{3l}$, the mixing angles $\theta_{ij}$ and CP phase $\delta_{CP}$, whereas the lightest neutrino mass $m_{\nu_{0}}$ remains unknown. There are upper bounds from cosmology on the sum $\sum_{i}m_{i}$ that depend on the active neutrino mass hierarchy, the data set used and the underlying cosmological model. The value quoted by the Planck collaboration for a standard cosmological model is $\sum_{i}m_{i}<0.12$ eV Aghanim:2018eyx , a discussion of how this value changes with different assumptions can e.g. be found in the Particle Data Group Report Tanabashi:2018oca . In fact, using the best fit values for the mass splittings from the NuFit data NuFit and the conservative value $\sum_{i}m_{i}<0.23$ eV, we can infer the upper limits of $m_{\nu_{0}}<7.12\times 10^{-2}$ eV for normal and $m_{\nu_{0}}<6.55\times 10^{-2}$ eV for inverted hierarchy.

The value of $m_{\nu_{0}}$ strongly impacts on the lower limit on $U_{I}^{2}$. One can obtain a reliable estimate of the lower bounds on $U_{I}^{2}$ by setting $\mathcal{R}=1$ Drewes:2019mhg . This makes the PMNS matrix unitary, and the lower limit one the smallest mixing can be estimated as $U_{I}^{2}\gtrsim m_{\nu_{0}}/M_{I}$. Using this approximation, we show in Figure 4 the lower limits on $U_{I}^{2}$ that we obtain in our scans for different values of $m_{\nu_{0}}=(0.05,10^{-2},10^{-3},10^{-4})$ eV. In the case of $m_{\nu_{0}}=0$ there is no absolute lower limit on the coupling from the seesaw mechanism, and the residual lower limit on $U_{I}^{2}$ is due to the BBN constraint.

The lightest neutrino mass has an important effect on the pattern of flavour mixing. In the limit of large $m_{\nu_{0}}$, there is almost no constraint on the allowed flavour ratios $U_{\alpha I}^{2}/U_{I}^{2}$. This is shown in Figure 11 by the black solid (dashed) contours, which indicate the allowed region within $1\sigma$ ($2\sigma$) where the lightest neutrino mass is $m_{\nu_{0}}$ < 10 meV (close to the cosmological bound stated above). In this case, there is no visible upper limit on $U_{\mu I}^{2}/U_{I}^{2}$ or $U_{\tau I}^{2}/U_{I}^{2}$ for normal hierarchy, whereas $U_{eI}^{2}/U_{I}^{2}$ is constrained $\lesssim 0.95$. Conversely, for inverted hierarchy there is an upper limit for the $\mu$ and $\tau$ flavours, but none for the $e$ flavour. However, for smaller values of $m_{\nu_{0}}$, the allowed range for the flavour mixing pattern becomes significantly constrained.121212 When constraining $m_{\nu_{0}}$ to very small values, we almost decouple one right handed neutrinos. The contribution of this feebly coupled state to the generation of light neutrino masses is negligible, which in return implies that its properties are almost unconstrained by neutrino oscillation data, and such is its flavour mixing pattern. Thus extreme ratios $U_{\alpha I}^{2}/U_{I}^{2}$ can in principle occur for this particular heavy neutrino, although the absolute values of $U_{I}^{2}$ remains negligible, and it has no effect on any near future experiment. Since our focus is primarily on heavy neutrinos that make a measurable contribution to the generation of light neutrino masses and/or may be discovered in experiments, we applied a cut on $M_{I}U_{I}^{2}>10^{-10}$ GeV in Figure 11 and Figure 12 to remove artefacts arising from states that are practically decoupled. The value of the cut is motivated by experimental sensitivities as demonstrated. The NA62 experiment, for instance, will only be able to probe up to sensitivities of $M_{I}U_{I}^{2}\approx\mathcal{O}(10^{-8})$ Drewes:2018gkc ; under optinistic assumptions the LHC may test $M_{I}U_{I}^{2}\approx\mathcal{O}(10^{-8})$ Drewes:2019fou and the FCC $M_{I}U_{I}^{2}\approx\mathcal{O}(10^{-11})$ Antusch:2017pkq . This is shown by the lines for $m_{\nu_{0}}<1$ meV (blue), $m_{\nu_{0}}<0.1$ meV (green) and $m_{\nu_{0}}<0.01$ meV (red). For masses lower than 0.01 meV the constraints saturate and the size of the ellipse remains almost constant. This can be also seen in Figure 12, where the largest coupling ratio is plotted for each flavour as function of neutrino mass.

It is instructive to compare our results to the constraints on the flavour mixing pattern in the scenario with $n=2$ that were found in Refs. Hernandez:2016kel ; Drewes:2016jae ; Drewes:2018gkc . For this purpose it is not sufficient to simply insert very small values for $m_{\nu_{0}}$ in the parameterisation (29) because such values can also be achieved due to accidental cancellations in the light neutrino mass matrix (without decoupling of any of the heavy neutrinos), cf. Section 2.6. To remove such fine tuned points we impose the following cuts

[TABLE]

Here $\epsilon$ is an arbitrarily small number, which we choose as $\epsilon=0.01$ for convenience. In addition, we work in the limits as defined by $|\mathrm{Im}\omega_{23}|\gg 1$ and $\mathrm{Re}\omega_{13}\sim\pi/2$ for normal hierarchy, and $|\mathrm{Im}\omega_{12}|\gg 1$ for inverted hierarchy (c.f. Appendix C). Note that we randomised the order of the matrices $\mathcal{R}^{ij}$, and hence for normal hierarchy we can only reproduce the true symmetry protected regime for the permutation $\mathcal{R}=\mathcal{R}^{23}\mathcal{R}^{13}\mathcal{R}^{12}$. The inverted hierarchy limit is independent of permutations as two of the $\omega_{ij}$ are zero. In Figure 13, we show the triangle plots with 1$\sigma$ and 2$\sigma$ contours for NH and IH in the symmetry protected region after applying the aforementioned cuts to remove fine-tuned points. The results are consistent with what was found in Ref. Drewes:2018gkc for $n=2$ RHNs. It is worth noting that there is a sharp upper limit on $U_{eI}^{2}/U_{I}^{2}$ where the contours do not show. This is due to the hard upper limit imposed on $m_{\nu_{0}}$ in order to reach the $n=2$ case and it is, as before, consistent with the results in Drewes:2018gkc .

5.4 Discussion of excesses likelihoods

In the previous subsections, we have made use of an exclusion-only ‘capped’ profile likelihood to study the constraining effect of the various observables on the parameter space (for a justification see Sec. 4.5). The total likelihood, however shows a pattern of excesses in some small regions of the parameter space. As discussed in Sec. 4.5, experimental results with a preference for specific heavy neutrino masses and mixings would in general not show up as localized excesses in the total profile likelihood. This is due to the fact that for each value of $M_{I}$ it would be in general possible to find a value of $M_{J}$ (with $J\neq I$) and associated couplings that would maximize the excess likelihood, irrespective of the values of $M_{I}$. In order to extract the specific masses and couplings preferred by an excess likelihood, we adopt throught this subsection the following strategy. We only allow one of the three RHNs (which we take to be $I=1$) to acquire the required masses and couplings, while disallowing the other two RHNs to enter the preferred region. This is emphasized in the plots by specifying $M_{1}$ and $|U_{\alpha 1}|^{2}$ instead of $M_{I}$ and $U_{\alpha I}^{2}$. Mind that these results would be identical for $M_{2}$ and $M_{3}$.

The invisible width of the $Z$ boson is modified by the presence of the right-handed neutrinos through their mixing, as described in Section 3.2.1. For very high $\tau$ couplings, $U_{\tau I}^{2}>10^{-3}$, the prediction from the RHN model is actually a better fit to the experimental measurement than the SM, and thus there is a slight ($<2\sigma$) excess. A similar effect occurs for the CKM and $R^{\tau}_{e\mu}$ constraints, where the modified contribution on the neutrino mixing in the decay products of $K$-mesons and $\tau$, enhances the prediction with respect to that of the SM. Figure 14 shows the excesses on the total profile likelihood in the $M_{1}$ vs $U_{\tau I}^{2}$ plane, zoomed in at high couplings (as discussed above, we excluded $M_{2}$ and $M_{3}$ from entering the excess regions). Since there are no constraints from direct searches at masses above $M_{1}>80$ GeV or in the range $0.3<M_{1}<0.5$ GeV, there is a combined excess shown of about $2\sigma$.

In order to study the impact of the different partial likelihoods on the total likelihood excess, we show in Figures 15 and 16 the partial one-dimensional likelihoods for $\Gamma_{\rm{inv}}$ (blue), CKM (green) and $R^{\tau}_{e\mu}$ (pink) with respect to the total likelihood (red) for $M_{1}<1$ GeV and $M_{1}>60$ GeV, respectively. All likelihoods are normalised so that they show up as a bump over the combination of all other likelihoods $\mathcal{L}_{0}$ (grey). These plots show that the combination of excesses from all three sources amounts to a deviation of around (high mass) or above (low mass) $2\sigma$ with respect to the background. As observed in the figures, the effect of $R_{e\mu}^{\tau}$ is rather negligible compared to the other two relevant likelihoods. Even larger couplings are severely penalised by the steep drop in the $\Gamma_{\rm{inv}}$ likelihood.

The excesses shown in Figures 14–16 in $|U_{\tau 1}|^{2}$, for both low and high masses, are the most significant excesses arising in our three RHN scenario, but not the only ones. At masses around the $K$-meson resonance, there is an even dimmer excess in $|U_{e1}|^{2}$, arising from the constraint on fully leptonic decays of $K$-mesons, $R^{K}_{e\mu}$. As seen in Figure 17, for both normal and inverted hierarchy, there is a $\sim 1\sigma$ excess at $M_{1}\sim 0.45$ GeV. As before, we show in Figure 18 the one-dimensional likelihoods for $R_{e\mu}^{K}$ (purple) with respect to the total likelihood (red), over the background of the combination of the rest of constraints (grey). Although the $R^{K}_{e\mu}$ likelihood keeps increasing for larger values of $|U_{e1}|^{2}$, the total likelihood drops at the limit shown in the figures due to the constraints from the CHARM experiment (orange).

Although the identified excesses provide interesting hints towards specific regions of the RHN parameter space, they should not be over-interpreted, since their significance remains rather small and probably consistent with statistical fluctuations. The presence of such excesses was already observed before, identified in EWPOs Antusch:2014woa (cf. also Akhmedov:2013hec ) and CKM constraints, and particularly in $\tau\to s$ transitions Drewes:2015iva .

6 Conclusions & Outlook

We presented here the first frequentist global analysis of the extension of the Standard Model by three heavy right-handed Majorana neutrinos for a large range of their masses, from 60 MeV to 500 GeV, and for normal and inverted hierarchy of the active neutrino masses. As detailed in Section 1.2, our analysis improves on previous studies in numerous ways. Most notable is the inclusion of a larger number of experimental constraints than in previous studies, such as EWPOs, all LFV decay channels, active neutrino mixing and masses, as well as many direct searches. Furthermore, we have performed a proper statistical combination of all constraints using a composite likelihood approach, and studied the overall constraints on the parameter space using robust profile likelihood methods. To this end, we have used the advanced BSM inference tool GAMBIT gambit , which we appropriately extended with the relevant model specifications and experimental constraints.

The results shown in Section 5 cover the full studied mass range for all couplings down to $U_{I}^{2}\sim 10^{-16}$. The profile likelihood contours are consistent with the results found in previous studies. The upper limits on the heavy neutrino mixing with electron and muon flavour mostly follow the confidence levels provided by direct search experiments. In the projection of the likelihoods on the $M_{I}$-$U_{\alpha I}^{2}$ planes the interplay becomes visible only in the constraints on the third generation and for masses below a GeV. This is qualitatively different from the model with only two heavy neutrinos ($n=2$), where combination of direct, indirect and cosmological bounds imposes stronger constraints than each of them individually on the mixing with all three SM flavours, and this interplay can rule out a considerable mass region below the kaon mass Drewes:2016jae . We for the first time studied the global constraints on the heavy neutrino flavour mixing pattern, which strongly depends on the mass of the lightest SM neutrino $m_{\nu_{0}}$. We explicitly studied the limit of vanishing lightest neutrino masses, where we have shown that the flavour mixing pattern becomes significantly constrained for small values of $m_{\nu_{0}}$. For $m_{\nu_{0}}<$ 0.01 meV these constraints become independent of the precise value of $m_{\nu_{0}}$ in both mass hierarchies, which suggests that one heavy neutrino has effectively decoupled. In this regime we demonstrated that one can recover the results that have previously been found in the model with only two RHN in earlier works.

Furthermore, we identified a few excesses in the profile likelihood, which are due to the invisible decay width of the $Z$-boson, the CKM unitarity constraint and $R_{e\mu}^{K}$. Our best fit has a significance (w.r.t. SM) slightly above $2\sigma$. Although these excesses are not significant enough to favour the $n=3$ right-handed neutrino model in favour of the SM at the moment, an improvement on the measurements of the relevant observables will increase/decrease their significance in the future. Future $e^{+}e^{-}$ colliders, such as the ILC, FCC-ee or CEPC, might measure EW observables, including the $Z$ decay width, with higher precision Fan:2014vta than the current value from LEP Abbiendi:2000hu . The NA62 experiment, which targets kaon decays, might be able to improve the measurments of the CKM matrix elements $V_{us}$ and $V_{ud}$, as well as the lepton universality ratio $R_{e\mu}^{K}$ through more precise measurement of the fully leptonic decays of kaons Goudzovski:2010uk .

Since the strongest constraints on the absolute value of the couplings come from direct searches, it is expected that the results obtained in this analysis will change significantly with the next generation of direct search experiments. An overview of projected sensitivities can e.g. be found in Refs. Antusch:2016ejd ; Beacham:2019nyx ; Alimena:2019zri . Many of these searches can be performed at existing facilities, including the LHC, NA62, T2K or the DUNE near detector. The sensitivity of the LHC will soon be upgraded with the recently approved FASER experiment Feng:2017uoz and other proposed dedicated detectors Chou:2016lxi ; Kling:2018wct ; Gligorov:2017nwh ; Curtin:2018mvb ; Dercks:2018wum ; Alpigiani:2018fgd ; Helo:2018qej . In the more distant future the SHiP experiment Alekhin:2015byh ; Anelli:2015pba can search for heavy neutrinos in the GeV mass range SHiP:2018xqw , while future folliders such as FCC Abada:2019zxq or CEPC CEPCStudyGroup:2018ghi can explore larger masses. These experimental perspectives make the study of right handed neutrinos an exciting topic for the years to come. Additional motivation for such searches comes from cosmology because the baryon asymmetry of the universe can be explained by low scale leptogenesis for all experimentally allowed values of the mixing angles in the model considered here if the heavy neutrino masses lie below the electroweak scale Akhmedov:1998qx ; Abada:2018oly . If any heavy neutral leptons are found in experiments then our results for their properties, such as the flavour mixing pattern as a function of light neutrino parameters, provide a powerful test to assess whether these particles are responsible for the generation of light neutrino masses and/or the baryon asymmetry of the universe Chun:2017spz , and to distinguish the model with three heavy neutrinos considered here from the model with two heavy neutrinos or other extensions of the SM.

Acknowledgements

The authors would like to thank S. Antusch, C. Weiland and R. Ruiz as well as the rest of the GAMBIT Community for helpful discussions and comments. We acknowledge PRACE for awarding us access to Marconi at CINECA, Italy, and the Red Española de Supercomputación (Spain; FI2016-1-0021), for access to MareNostrum, Spain. Part of this work was carried out on the Dutch national e-infrastructure (Lisa cluster) with the support of SURF Cooperative. The authors are also grateful for the computing grant on the Prometheus computer from PLGRID Infrastructure (Poland). The work of M.C. is funded by the Polish National Agency for Academic Exchange under the Bekker program. M.C. is also grateful to Foundation for Polish Science (FNP) for its support. The work of S.K. and C.W. is funded by NWO through the Vidi research grant 680-47-532. T.E.G is funded by the Australian Research Council (Discovery Project DP180102209). The work of J.H. is funded by the DFG Emmy Noether Grant No. HA 8555/1-1, and was earlier supported by the Labex ILP (reference ANR-10-LABX-63) part of the Idex SUPER, and received financial state aid managed by the Agence Nationale de la Recherche (ANR), as part of the programme Investissements d’avenir under the reference ANR-11-IDEX-0004-02.

Appendix A GAMBIT Implementation

GAMBIT 131313gambit.hepforge.org. (the Global and Modular BSM Inference Tool) gambit is a global fitting software framework that allows for extensive calculations of observables and likelihoods in particle and astroparticle physics. It provides, out-of-the-box, a suite of statistical methods and parameter scanning algorithms, together with a hierarchical model database, a strong interface to external tools and a host of other utilities that make it one of the most powerful global fitting tools on the market.

In a nutshell, the fundamental building blocks of GAMBIT are its module functions, which calculate all physical and mathematical quantities revelant to an analysis. Each module function provides a capability which, together with the return type of the function, unequivocally specifies the quantity calculated.

These module functions are sorted in the different physics modules, according to their purpose, e.g. functions calculating dark matter relic density lie in DarkBit DarkBit . Since most observables and likelihoods computed for this analysis do not belong naturally to any of the existing GAMBIT modules, we introduce a new GAMBIT physics module, NeutrinoBit, which contains all calculations relating to (active and sterile) neutrino physics. All of the module functions and capabilites described below are implemented in the new module NeutrinoBit, unless otherwise stated.

A.1 Neutrino models

In GAMBIT, models are defined by a set of parameters and relations to other models gambit . All the SM and active neutrino parameters are defined in a model called StandardModel_mNudiff, a daughter model of the GAMBIT model StandardModel_SLHA2, which includes SM parameters written in the SLHA2 convention Allanach:2008qq . The former contains the parameters mNu_light, dmNu21 and dmNu3l, which give the lightest neutrino mass and mass splittings, respectively. Other parameters in this model that are relevant for this study and are scanned over, as described in Section 4, are alpha1, alpha2, delta13, theta12, theta23 and theta13.

The right-handed neutrino sector is defined in another model, RightHandedNeutrinos, which contains the RHN masses $M_{I}$ and the real and imaginary parts of the $\omega_{ij}$ parameters in the C-I parametrisation (c.f. Section 2.4). These are M_1, M_2, M_3, ReOm12, ReOm13, ReOm23, ImOm12, ImOm13, ImOm23. To better explore the symmetry preserved region (Section 2.5), a differential model, inheriting from RightHandedNeutrinos, has been defined, RightHandedNeutrinos_diff. This model swaps the parameters M_2 for delta_M21, and defines a translation function from the parameters of the daughter model to the parent model as

[TABLE]

Lastly, the parameter Rorder encodes the ordering of the matrices $\mathcal{R}^{ij}$ in Eq. (30), which allows us to fully cover all the parameter space with the C-I parametrisation.

There is a number of useful quantities and observables that can be constructed from the neutrino parameters, and these are all implemented in NeutrinoBit.cpp. In the active neutrino sector we calculate the neutrino mass matrix, m_nu, and their mixing matrix, UPMNS, as well as useful quantities such as the type of hierarchy, ordering, the squared mass splittings, md21, md31 and md32, and the minimal neutrino mass min_mass. It is worth noting that it is possible to fix the hierarchy of a scan by providing the option ordering to the capability m_nu in the configuration file. For example, in order to scan the normal hierarchy we would define in the YAML file

⬇

Rules:

• capability: m_nu

  options:

  ordering: 1

The right-handed neutrino sector also contains a couple of useful capabilities, SeesawI_Vnu, which is the active neutrino mixing matrix in type-I seesaw, effectively UPMNS corrected by the presence of the right-handed neutrinos, and SeesawI_Theta, the active-sterile mixing matrix in type-I seesaw, currently implemented using the C-I parametrisation.

Another useful capability defined here is Unitarity, which is fulfilled by two module functions according to whether the model scanned is the SM or a RHN model, and checks whether the full mixing matrix is unitary. All these capabilities relating to neutrino masses and mixings and the module functions that fulfill them, along with their dependencies and options can be seen in Table 8.

Lastly, NeutrinoBit.cpp also contains likelihoods for the active neutrinos, which are implemented following the results from the NuFit collaboration (c.f. Section 3.1). The capabilities associated with these are md21_lnL, md3l_lnL for the mass splittings, deltaCP_lnL, theta12_lnL, theta23_lnL and theta13_lnL for the phases and mixing angles, and sum_mnu_lnL for the cosmological limit on the sum of neutrino masses. All these capabilities, with their module functions and dependencies are listed in Table 9.

A.2 Right-handed neutrino likelihood functions

Every observable and likelihood described in Section 3 has an assigned capability within GAMBIT. Most of these have been implemented in the new GAMBIT module NeutrinoBit, since they concern mostly neutrino physics. Their module functions are coded in the file RightHandedNeutrinos.cpp, to keep them separated from the likelihoods and observables concerning only active neutrinos in NeutrinoBit.cpp. The exception to this is the LFV observables and semileptonic lepton universality tests, which can be found in FlavBit FlavBit , implemented in FlavBit.cpp and the electroweak precision observables, which are coded up in PrecisionBit.cpp in PrecisionBit SDPBit .

The implementation details for each specific observable are as follows:

*Electroweak precision observables

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

*Lepton flavour violation

*The capabilities related to lepton flavour violation can be found in FlavBit and are muegamma, tauegamma, taumugamma, mueee, taueee, taumumumu, taumuee, taueemu, tauemumu, taumumue, mueTi, mueAu and muePb. Table 11 shows these capabilities, the module functions that provide them, implemented in FlavBit.cpp, and their dependencies. The likelihoods, shown in Table 12, are collated into three capabilites, according to the type of process, l2lgamma_lnL for $l\to l\gamma$, l2lll_lnL for $l^{-}\to l^{-}l^{-}l^{+}$ and mu2e_lnL for $\mu-e$ conversion in nuclei.

*Lepton universality

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

*CKM unitarity

*The NeutrinoBit capability calc_Vus, implemented in RightHandedNeutrinos.cpp, calculates the value of $V_{us}$ that maximizes the likelihood for a given mixing matrix $\Theta$. The capabilities lnLckm_Vus and lnLckm_Vusmin compute the log-likelihood using $V_{us}$ as a scan parameter and as calculated by the profiling of calc_Vus, respectively. The capabilities, module functions and dependencies defined in GAMBIT for the calculation of the observable and the likelihood connected to CKM unitarity are listed in Tab. 15.

*Neutrinoless double beta decay

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

*Big Bang nucleosynthesis

*There are a number of processes that contribute to the decay width of the right-handed neutrinos, relevant for Big Bang nucleosynthesis, and each of them is computed by a capability. These are Gamma_RHN2piplusl, Gamma_RHN2Kplusl, Gamma_RHN2Dplusl, Gamma_RHN2Dsl, Gamma_RHN2Bplusl, Gamma_RHN2Bcl, Gamma_RHN2pi0nu, Gamma_RHN2etanu,  Gamma_RHN2etaprimenu, Gamma_RHN2etacnu, Gamma_RHN2rhoplusl, Gamma_RHN2Dstarplusl, Gamma_RHN2Dstarsl, Gamma_RHN2rho0nu, Gamma_RHN2omeganu, Gamma_RHN2phinu, Gamma_RHN2Jpsinu, Gamma_RHN23nu, Gamma_RHN2llnu, Gamma_RHN2null, Gamma_RHN2nuuubar, Gamma_RHN2nuddbar and Gamma_RHN2ludbar. The total decay width of each of the right-handed neutrinos is given by Gamma_BBN and the log-likehood for BBN by lnL_bbn. Tab. 17 shows the capabilities, as defined in GAMBIT, that pertain to Big Bang nucleosynthesis, and the module functions that satisfy them, along with dependencies that other module functions fulfill. The decay process considered in each function is mentioned below its name.

*Direct searches

*As detailed in Section 3.3, the likelihoods for direct searchs depend on the active-sterile matrix elements $U_{\alpha I}$. Hence, for simplicity the capabilities Ue1, Ue2, Ue3, Um1, Um2, Um3, Ut1, Ut2 and Ut3 are implemented in NeutrinoBit, as well as the phases of each of the matrix elements Ue1_phase, Ue2_phase, Ue3_phase, Um1_phase, Um2_phase, Um3_phase, Ut1_phase,Ut2_phase and Ut3_phase. These can be seen in Table 18 where I=1,2,3. All the capabilities UaI can take a pair of options upper_limit and lower_limit to force the values within the given range. Using these quantities, the likelihoods for the different direct search experiments are calculated, and their capabilities are lnLpienu , lnLps191e, lnLps191mu, lnLcharme, lnLcharmmu, lnLcharmtau, lnLdelphi_shortlived, lnLdelphi_longlived, lnLatlase, lnLatlasmu, lnLe949 and lnLnutev. The capabilities, module functions and their dependencies for all relevant direct search experiments are tabulated in Tab. 19.

*Other capabilities

*The theoretical constraint for perturbativity of the Yukawa couplings has been implemented in NeutrinoBit as well. The capability perturbativity_lnL calculates a step function likelihood for this constraint. Tab. 20 shows the module function that provides this capability and its dependencies.

Lastly, the artificial coupling slide likelihood that was introduced to drive the scan towards high couplings, as described in Section 4, was also implemented in NeutrinoBit with capability coupling_slide. The module function and dependencies of this capability can also be seen in Table 20.

Appendix B Full expressions for the relevant observables

B.1 Decay widths and form factors for LFV observables

The decay widths of LFV processes, as described in Section 3.2.2, are given by Kuno:1999jp ; Abada:2014kba

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where we used $\hat{A}^{V}_{XY}\equiv A^{V}_{XY}+e^{2}K_{1}^{X}$. The couplings $e$, $g_{1}$, $g_{2}$ correspond to the electromagnetic, hypercharge and weak couplings of the SM.

The form factors $K_{1}^{X}$, $K_{2}^{X}$, $A^{S}_{XY}$ and $A^{V}_{XY}$ are taken in the flavour basis where the charged lepton mass matrix is diagonal.

The dipole form factors $K_{1}^{X}$ and $K_{2}^{X}$ are given as Abada:2014kba

[TABLE]

The four lepton form factors $A^{V}_{XY}$ and $A^{S}_{XY}$ corresponding to the process $l^{-}_{\alpha}\to l^{-}_{\beta}l^{-}_{\gamma}l^{+}_{\delta}$, with a vector or scalar mediator respectively, are Abada:2014kba

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where sum over $a$ and $c$ is assumed, $x_{a}=m_{\nu_{a}}^{2}/m_{W}^{2}$, $s_{w}=\sin\theta_{w}$ and $g_{\pm}=g_{1}\sin\theta_{w}\pm g_{2}\cos\theta_{w}$.

The $\mu-e$ conversion ratio described in section 3.2.2, for a general nucleus $N_{Z}^{Z+N}$, can be written as

[TABLE]

where $B_{XY}^{K}\leftrightarrow C_{XY}^{K}$ for up-type quarks and the numerical factors $G_{K}$ are given in Kosmas:2001mv . The nuclear form factor $F_{p}$, the effective atomic number $Z_{\text{eff}}$ and the capture rate $\Gamma_{\text{capt}}$ of the nucleus Alonso:2012ji can be seen, for the cases we are studying, $\text{Ti}_{22}^{48}$, $\text{Au}_{79}^{197}$ and $\text{Pb}_{82}^{208}$, in Table 21.

The form factors $K_{1}^{X}$, $K_{2}^{X}$ are defined above and $B_{XY}^{K}$ and $C_{XY}^{K}$ are given by

[TABLE]

with $g_{d}=3g_{2}\cos\theta_{w}+g_{1}\sin\theta_{w}$, $g_{u}=-3g_{2}\cos\theta_{2}+g_{1}\sin\theta_{w}$, $x^{(u,d)}_{c}=m_{(u,d)_{c}}^{2}/m_{W}^{2}$ and $V_{ij}$ is the CKM matrix of quark mixing.

Lastly, the loop functions used in (82)-(99) are defined as Abada:2014kba

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

B.2 Decay widths relevant for Big Bang Nucleosynthesis

The various decay widths of RHNs, relevant for the BBN as described in Sec. 3.2.6, are listed here. These expressions are taken from Gorbunov:2007ak ; Atre:2009rg ; Canetti:2012kh ; Bondarenko:2018ptm ; Ballett:2019bgd , among which there are slight differences that will be commented upon when relevant. We list here the decays for Majorana fermions, which differ by a factor of 2 with respect to the rates for Dirac fermions as shown, for instance, in  Bondarenko:2018ptm .

The decay width of a RHN, $N_{I}$, to a lepton, $l_{\alpha}$ and a charged pseudoscalar meson, $P^{+}=\pi^{+},K^{+},D^{+},D_{s},B^{+},B_{c}$, is Gorbunov:2007ak ; Atre:2009rg ; Bondarenko:2018ptm ; Ballett:2019bgd

[TABLE]

with $f_{P}$ the decay constant of the meson $P^{+}$, which can be seen in Table 22, $V_{P}$ the CKM matrix element corresponding to $P^{+}$ and

[TABLE]

Similarly, the decay width of a RHN to a neutrino and a neutral pseudoscalar meson, $P^{0}=\pi^{0},\eta,\eta^{\prime},\eta_{c}$, is Gorbunov:2007ak ; Bondarenko:2018ptm ; Ballett:2019bgd

[TABLE]

with $f_{P}$ the meson decay constant, Table 22, and $x_{P}$ as in (110). The expression from Atre:2009rg for this decay missed a factor of 2, which was corrected by the later work Ballett:2019bgd .

The decay width of a RHN to a lepton and a charged vector meson, $V=\rho^{+},D^{*+},D_{s}^{*}$, is Atre:2009rg ; Bondarenko:2018ptm ; Ballett:2019bgd

[TABLE]

where $f_{V}$ is the decay constant of $V^{+}$, in Table 22, $V_{V}$ is the CKM matrix element associated with $V^{+}$, $\lambda$ is defined in (111) and

[TABLE]

The results from Bondarenko:2018ptm use a different definition of decay constant $g_{V}=m_{V}f_{V}$, but the final values agree nevertheless. The value of $g_{\rho}$ in Gorbunov:2007ak differs from that of Bondarenko:2018ptm and there is a factor of $\tfrac{1}{2}$ missing as well with respect to the other works.

The last of the semileptonic decays of RHNs is the decay to a neutrino and a neutral vector meson, $V^{0}=\rho^{0},\omega,\phi,J/\psi$, with decay width Bondarenko:2018ptm

[TABLE]

with $f_{V}$ the meson decay constant, Table 22, $x_{V}$ as in (114) and $\kappa_{V}$ a neutral current correction factor Bondarenko:2018ptm

[TABLE]

The expressions from Bondarenko:2018ptm , Atre:2009rg and Ballett:2019bgd agree but for a different definition of $\kappa_{V}$, whereas Gorbunov:2007ak misses the $\kappa_{V}$ factor altogether.

The fully-leptonic three body decays of RHNs can be to three neutrinos, with decay width given by Gorbunov:2007ak ; Bondarenko:2018ptm

[TABLE]

and to two charged leptons and a neutrino. If the charged leptons have the same flavour, the decay width is Gorbunov:2007ak ; Bondarenko:2018ptm

[TABLE]

with $x_{l}$ as in (109) with $\alpha\leftrightarrow\beta$, the coefficients $C_{i}$ are

[TABLE]

and the functions $S(x,y)$ and $L(x)$

[TABLE]

If, on the other hand, the charged leptons are of different flavour, the decay width is given by Canetti:2012kh

[TABLE]

where $x_{\alpha}$ and $x_{\beta}$ are as in (109), $S(x,y)$ is as in (123) and

[TABLE]

The expressions for leptonic decays agree across all works, except minor differences in the functional forms of the decays to different flavour leptons between Canetti:2012kh and Bondarenko:2018ptm . These differences are, however, numerically negligible.

For large masses, above the hadronisation scale, the full hadronic decay width of the right-handed neutrinos is better approximated by computing their decay to free quarks, instead of the individual meson channels. Hence the decays of RHNs to quarks through the neutral current is Bondarenko:2018ptm

[TABLE]

where the 3 upfront accounts for the number of colours, $x=M_{q}/M_{I}$, $L(x)$ as in Eq.(124), and

[TABLE]

for up-type quarks and

[TABLE]

for down-type quarks.

And lastly, the decay of RHNs to free quarks through the charged current can be written as Bondarenko:2018ptm

[TABLE]

where $x_{i}=M_{i}/M_{I}$ and the function $I(x_{u},x_{d},x_{l})$ is

[TABLE]

The differences between the expressions for the decay to free quarks between Canetti:2012kh and Bondarenko:2018ptm mirror those of the leptonic decays. In this case, however, the differences can be substantial for decays to third generation quarks, in which case we opt to use the expressions from Bondarenko:2018ptm .

Appendix C Distinguishing symmetry protected from tuned parameter choices

One goal of the present work is to fully understand the experimentally allowed range of parameters for heavy neutrinos, with minimal theoretical bias. To achieve this we employ the Casas-Ibarra parametrisation (29) and adapt agnostic priors for the parameters in table 5. On the other hand it is also instrictive to understand what fraction of the parameter space can only be realised at the cost of fine tuning in the parameters. This requires to distinguish fine-tuned parameter choices from symmetry protected ones.

The Casas-Ibarra parametrisation (29) is inherently motivated from “bottom up”, and it is not easy to see directly from the values of its fundamental parameters whether they exhibit a symmetry protection. A full analytic exploration of all the possible solutions and their classification between symmetric and fine-tuned would be a useful exercise, but lies outside the scope of this work. In our numerical scan we take a more pragmatic approach. We first generate a huge amount of parameter choices by randomising the parameter values and the order of the matrices $\mathcal{R}^{ij}$ to ensure a maximal coverage of the parameter space. We then use the cuts (5.3) to distinguish the symmetry protected points a posteriori. This cut practically enforces the structure (2.5) on the masses and couplings.

Using the cut (5.3) requires some care for two reasons. First, the form (2.5) does not capture all symmetry protected points, cf. footnote 6. We may therefore misidentify some symmetry protected points as tuned. We find, however, that the number of such points is small. Second, when using the Casas-Ibarra parametrisation, it is possible to generate points that mimic the form (2.5) and hence pass the cut (5.3), but in fact exhibit a significant amount of tuning.

To illustrate the second point we work at tree level and approximate $\tilde{M}^{\text{diag}}\simeq M_{M}$, which yields

[TABLE]

For the inverted hierarchy, it is straightforward to show that one qualitatively gets a pseudo-Dirac pair of heavy neutrinos for

[TABLE]

with $|{\rm Im}\omega|\gg 1$. In addition we have to set $m_{\nu_{0}}=0$ to find the symmetry protected region, as a non-zero lightest neutrino mass is not consistent with $\epsilon_{\alpha}^{\prime}\to 0$. In this case the upper left block of the matrix $\mathcal{R}\sqrt{M_{M}}$ in (29) is large and the third row that multiplies $m_{3}=m_{\nu_{0}}$ is small, thereby mimicking the structure in (2.5). One can therefore interpret the decoupling $F_{\alpha 3}=0$ and the vanishing mass of the lightest neutrino physically as results of the symmetry.

If we choose normal ordering, then choosing (135) still yields a structure that passes the cut (5.3), but it is in fact a tuned solution that just mimics this structure. In that case $m_{\nu_{0}}=m_{1}$ multiplies the large components of $\mathcal{R}\sqrt{M_{M}}^{-1}$ in (29). Hence, the approximate symmetry makes the wrong light neutrino mass small ($m_{3}$ instead of $m_{\nu_{0}}=m_{1}$). Of course one can set $m_{1}=0$ by hand in (29), but this choice cannot be justified by the symmetry. Though the limit (135) leads to a pseudo-Dirac structure amongst the $\nu_{Ri}$ as predicted by the $B-\bar{L}$ symmetry, the vanishing mass of the lightest neutrino is not a result of that symmetry.

The problem is that the Casas-Ibarra parametrisation allows on to set $m_{\nu_{0}}=0$ by hand and gives no warning if this leads to accidental cancellations.

Realising the symmetry requires to choose the eigenvalues of $M_{M}$ and the non-zero $\omega_{ij}$ consistently in a way that the large block in the matrix $\mathcal{R}\sqrt{M_{M}}^{-1}$ in (29) multiplies the two non-zero light neutrino masses. For normal ordering this is achieved with

[TABLE]

again with $m_{\nu_{0}}=0$. In particular, one cannot choose $\nu_{R3}$ to be the particle that decouples. This is clearly no fundamental problem because the labels of the $\nu_{RI}$ have no physical meaning, but it means that the labelling and the order of the matrices $\mathcal{R}^{ij}$ have to be taken into consideration when applying a cut to identify symmetry protected points in the numerical data.

The situation is yet more tricky if one considers small perturbations around the choices (136) or (135). In our numerical scan we randomise the order of the three matrices $\mathcal{R}^{ij}$ in (30) to generate more points. If one exactly takes the choice (136) or (135) for normal or inverted neutrino mass ordering, respectively, then the approximate $B-\bar{L}$ conserving limit is reproduced irrespectively of the ordering of the $\mathcal{R}^{ij}$. However, the effect that small perturbations around this limit have strongly depends on this ordering. The effect of perturbing $\mathcal{R}$ is the smallest if the matrices $\mathcal{R}^{ij}$ are ordered in a way that the one with large entries (controlled by $\omega$) directly multiplies $\sqrt{M_{M}^{-1}}$ in (134). For normal ordering this is the case with $\mathcal{R}=\mathcal{R}^{23}\mathcal{R}^{13}\mathcal{R}^{12}$, and for inverted ordering with $\mathcal{R}=\mathcal{R}^{12}\mathcal{R}^{13}\mathcal{R}^{23}$. This procedure was crucial to reproduce the constraints on the heavy neutrino flavour mixing pattern in the $n=2$ model found in ref. Drewes:2018gkc , cf. fig. 13.

It is worth noting that (135) and (136) are not the only combinations of parameters that yield the symmetry protected scenario, but rather the simplest. The non-trivial structure of the complex rotation matrix $\mathcal{R}$ yields many solutions to the required block layout described before. In fact, we will take advantage of this fact further below to recover the $n=2$ case from the $n=3$ Lagrangian for normal ordering by taking

[TABLE]

since in this work we will focus mainly on the case where $M_{1}$ and $M_{2}$ are almost degenerate.

Appendix D Partial likelihoods

The final result of a global fit shows the combined effect of all likelihoods on the parameter space of the model. It is, however, often useful to understand the effect on the individual partial likelihoods. Therefore, we show here a comprehensive set of scatter plots that show the contribution of each relevant partial likelihood in the $M_{I}$ vs $|U_{\alpha I}|^{2}$. In all figures throughout this section the colourbar measures the relative partial log likelihood for each observable with respect to the global best fit value. As we have seen before, away from the massless neutrino limit there is little difference between NH and IH, and thus we will only show the partial likelihoods for normal ordering.

Figures 19 - 24 show the most constraining likelihoods on the $|U_{eI}|^{2}$ coupling. The likelihood values are normalised to the best fit value for each partial likelihood. Consistently with the results above, various direct searches constrain large values of the coupling, with PS191 dominating for $M_{I}\lesssim 0.45$ GeV, CHARM for $M_{I}\sim(0.45,2)$ GeV, the long-lived particle search from DELPHI for $M_{I}\sim(2,4)$ GeV and DELPHI prompt search for $M_{I}\sim(4,80)$ GeV.

As seen in Figure 23, direct searches from CMS compete in a small mass range with DELPHI prompt searches, the statistical combination of the two setting stronger limits than each of them individually. Recent and future results from CMS and ATLAS not included in this study are expected to dominate in this range.

Figure 24 shows that the larger mass range is unconstrained by direct searches, hence electroweak precision observables, in particular $\sin\theta_{W}$, are responsible for the upper limits in this range.

Similar to the case above, the coupling $|U_{\mu I}|^{2}$ is constrained from above by several direct and precision searches. Figures 25-31 show the effect of the individual likelihoods on the upper limit of $|U_{\mu I}|^{2}$. As opposed to the electron case, where for most mass ranges only one constraint dominated, in this case several mass ranges show competing effects from various constraints. For $M_{I}<0.45$ GeV both PS191 and E949 are relevant; in the range $M_{I}\sim(0.45,2)$ GeV searches at NuTeV are the most constraining, with a small contribution from the results from CHARM; the long-lived particle search from DELPHI remains unchallenged for $M_{I}\sim(2,4)$ GeV whereas, as before, the DELPHI prompt search competes in the range $M_{I}\sim(4,80)$ GeV, with searches at CMS.

Larger masses are not constrained by direct searches, but rather by a combination of precision limits. Contrary to $|U_{eI}|^{2}$, where only $\sin\theta_{W}$ dominated at large masses, upper values of $|U_{\mu I}|^{2}$ are also mildly constrained by lepton flavour violating decays, particularly $\mu\to e\gamma$. Hence, for $M_{I}\gtrsim 80$ GeV, the combination of EWPO, $\sin\theta_{W}$, and LFV decays, are the most constraining, as seen in Figures 32 and 33. The effect of the LFV constraints can be better appreciated in Figs. 34 and 35 on the combination $|U_{eI}U_{\mu I}|$ where $\mu\to e\gamma$ is the dominant constraint, supplemented slightly by $\mu-e$ conversion.

The couplings of heavy neutrinos to the $\tau$ flavour, $|U_{\tau I}|^{2}$, are not as strongly constrained from above by direct searches. In Figures 36-37, one can see that for low masses, $M_{I}\lesssim 0.3$ GeV, only the direct searches from CHARM in the $\tau$ channel set an upper limit on the couplings. In the mass range $M_{I}\sim(0.5,80)$ GeV, long-lived and prompt searches by DELPHI dominate. At low masses, the $|U_{\tau I}|^{2}$ coupling is constrained from below, as seen in Figure 39. This lower bound is a consequence of BBN, which sets a lower limit on the sum of couplings $|U_{I}|^{2}$, and PS191, which forces the $e$ and $\mu$ couplings to be small at low masses.

In the mass range $M_{I}\sim(0.3,0.5)$ GeV, as well as for large masses $M_{I}\gtrsim 80$ GeV, direct searches do not constrain $|U_{\tau I}|^{2}$. Hence in these ranges, the strongest constraints come from the invisible decay of the $Z$ boson, as seen in Figure 40. This figure uses the "capped" likelihood defined previously, so the excesses in $\Gamma_{\rm inv}$ discussed in Section 5 will not be visible.