Neutrinoless double beta decay versus other probes of heavy sterile neutrinos

We make a comparative study of the neutrinoless double beta decay constraints on heavy sterile neutrinos versus other direct and indirect constraints from both lepton number conserving and violating processes, as a sensitive probe of the extent of lepton number violation and possible interference effects in the sterile sector. We introduce a phenomenological parametrisation of the simplified one-generation seesaw model with one active and two sterile neutrino states in terms of experimentally measurable quantities, such as active-sterile neutrino mixing angles, CP phases, masses and mass splittings. This simple parametrisation enables us to analytically derive a spectrum of possible scenarios between the canonical seesaw with purely Majorana heavy neutrinos and inverse seesaw with pseudo-Dirac ones. We then go on to constrain the simplified parameters of this model from various experiments at the energy, intensity and cosmic frontiers. We emphasise that the constraints from lepton number violating processes strongly depend on the mass splitting between the two sterile states and the relative CP phase between them. This is particularly relevant for neutrinoless double beta decay, which is weakened for small mass splitting and opposite CP parities between the sterile states. On the other hand, neutrinoless double beta decay is especially sensitive for Majorana sterile neutrinos with masses around $0.1-10$ GeV.


Introduction
The observation of neutrino oscillations in solar, atmospheric, reactor and accelerator neutrino data [1] implies that at least two of the three active neutrinos have a small but non-zero mass and that individual lepton flavour is violated. In the Standard Model (SM), neutrinos have only one helicity state ν L , and therefore, cannot acquire a Dirac mass through the Higgs mechanism, unlike the charged fermions. A Majorana mass term of the formν C L ν L (where ν C L ≡ ν T L C −1 , C being the charge conjugation matrix) is also forbidden in the SM due to its gauge structure and particle content. Specifically, the SM does not contain an SU(2) L triplet Higgs which could give rise to theν C L ν L term. By adding a SM-singlet, right-handed (RH) neutrino field ν R per generation to the SM, one could in JHEP03(2020)170 principle generate a Dirac mass term; however, to get sub-eV left-handed (LH) neutrino masses as required by the neutrino oscillation data, one needs the Dirac Yukawa couplings to be 10 −12 . While theoretically allowed, such a scenario would be rather uninteresting from an experimental point of view. A more appealing choice is to break the accidental (B − L) symmetry of the SM to generate a Majorana neutrino mass at tree or loop-level.
The simplest tree-level realisation of the (B − L)-breaking is through the effective dimension-5 Weinberg operator (L T H)(L T H)/Λ, where L and H are the SU(2) L lepton and Higgs doublets of the SM and Λ is the scale of new physics that induces (B − L)breaking [2]. Here the intermediate heavy particles integrated out in the low-energy theory are SM-singlet fermions, identified as the RH Majorana neutrinos ν R,i (with i = 1, 2, · · · ) with mass m N i . This is widely known as the type-I seesaw mechanism [3][4][5][6][7]. In the minimal type-I seesaw extension of the SM, the RH Majorana neutrinos, also known as the sterile neutrino states (or heavy neutral leptons in some literature), being SM gauge-singlets, can only interact with the SM particles through their mixing with the active neutrino states. In the traditional 'vanilla' seesaw mechanism, this active-sterile neutrino mixing is given by due to the smallness of the light neutrino mass m ν 0.1 eV, as inferred from neutrino oscillation data [1], as well as the cosmological limit on the sum of active neutrino masses [8]. Thus for a low seesaw scale in the sub-TeV range, the experimental effects of the sterile neutrino are expected to be small, unless they have additional interactions, e.g. when they are charged under an extra gauge U(1) B−L . However, there exists a class of minimal SM plus low-scale type-I seesaw scenarios [9][10][11][12][13][14][15][16][17][18][19][20][21], where V N can be sizeable while still satisfying the light neutrino data. This is made possible by assigning specific textures to the Dirac and Majorana mass matrices. The stability of these textures can in principle be guaranteed by enforcing extra symmetries in the lepton sector [12,13,19,[22][23][24]. We will generically assume this to be the case for our subsequent discussion and freely vary the active-sterile mixing up to O(1), without referring to any particular texture or model-building aspects.
For the conventional seesaw scenarios mentioned above, the active neutrino masses are inversely proportional to the lepton number violating (LNV) Majorana mass scale m N (hence the name 'seesaw'). There exist important variations, such as inverse seesaw [25][26][27], linear seesaw [28][29][30][31] and generalised inverse seesaw [32][33][34], where the active neutrino masses are directly proportional to the lepton-number breaking scale. In such scenarios, large active-sterile neutrino mixing can be achieved rather naturally, irrespective of the sterile neutrino mass spectrum. Experimentally, the main distinguishing feature of these variants from pure type-I seesaw is the pseudo-Dirac nature of the sterile neutrinos, which suppresses the LNV signals, in contrast with the purely Majorana nature of the sterile neutrinos in the type-I seesaw scenario. However, this distinction may not be as clear in the presence of additional CP phases in the sterile sector, depending on the sterile neutrino mass spectrum [35][36][37], with important implications for collider physics and leptogenesis [38].
The aim of this paper is to show how both purely Majorana and pseudo-Dirac limits (and the spectrum of possible cases in between) can be understood from a simple phe-JHEP03(2020)170 nomenological parametrisation in terms of the a priori measurable mixing angles, CP phases and mass eigenvalues involving the sterile neutrinos. In order to show our results analytically as much as possible, we will work in a simplified single-generational picture (involving only the electron flavour) and introduce the generalised inverse seesaw with two SM-singlet Weyl fermions ν R,1 and ν R,2 . In this case, the unitary matrix V that diagonalises the full neutrino mass matrix M ν is a 3 × 3 matrix containing three mixing angles (ϑ e1 , ϑ e2 , ϑ 12 , with the first two being the mixing angles of the active neutrino with the two sterile states, while the last one being the mixing angle in the sterile sector) and three CP phases (δ, φ 1 , φ 2 , with δ being the Dirac CP phase and φ 1 , φ 2 being the Majorana phases). We will apply this general parametrisation to identify the regions of parameter space allowed by consistency relations among the neutrino mass matrix elements. We will look at how the Majorana and pseudo-Dirac limits of the sterile state pair depend on the phases φ 1 and φ 2 = φ 2 − 2δ and how these in turn are completely determined by the activesterile squared mixing strengths s 2 e1 ≡ sin 2 ϑ e1 and s 2 e2 ≡ sin 2 ϑ e2 as a result of the (1, 1) element of M ν being zero. We will also use the (1, 3) element, which can be set to zero by a particular rotation and therefore parametrisation of M ν , to constrain the sterile-sterile squared mixing strength s 2 12 ≡ sin 2 ϑ 12 and linear combination of phases δ = 2φ 2 + δ. The angle ϑ 12 is an unobservable parameter in the SM because it is not contained in the mixing strengths V eN i . Nonetheless, the solution for s 2 12 in the chosen parametrisation gives a parametrisation-independent value for the light active neutrino mass at one-loop, which we enforce to be considerably smaller than the tree-level mass.
In order to put the neutrinoless double beta (0νββ) decay sensitivity into context, we review the experimental constraints on the sterile neutrino sector from both lepton number conserving (LNC) and violating channels from high-energy collider searches, high-intensity beam dump and meson decay experiments, beta decays and other nuclear processes, activesterile neutrino oscillation experiments, electroweak precision data and other indirect laboratory searches, as well as cosmological and astrophysical observations. We pay particular attention to possible interference effects from two sterile states on the LNV constraints as a function of their mass splitting. We give special emphasis on the 0νββ decay constraint, which has been argued to be the most stringent one for active-sterile neutrino mixing in the electron flavour; see e.g. summary plots in refs. [39][40][41][42]. We re-evaluate the 0νββ decay constraints using the general parametrisation discussed above and show how these are affected by the sterile neutrino mass splitting and CP phases, in comparison to other laboratory constraints. In fact, under certain conditions, we find the 0νββ decay constraints to be weaker than the direct search limits from colliders, thus reinforcing the importance of independent direct searches for sterile neutrinos in all flavours. On the other hand, for the theoretically interesting mass regime m N ≈ 0.1 − 10 GeV, 0νββ decay is comparable to current and future direct searches. As we work in a simplified one-generation framework containing a single active neutrino state which we identify as the electron neutrino, we cannot model the coherent contribution of the other two light states to 0νββ decay. Our main focus is on the constraints on sterile neutrinos in a simplified yet consistent seesaw picture but we will comment on the omission of the other light states (see the discussion at the end of section 5.1).

JHEP03(2020)170
The rest of this paper is organised as follows. In section 2 we introduce the generalised inverse seesaw for the neutrino mass matrix M ν , which reduces to the type-I, inverse and linear seesaw scenarios under different limits. We also investigate the masses at tree-level and at one-loop for the light mostly-active neutrinos as a function of the model parameters. In section 3 we introduce a phenomenological parametrisation of the unitary matrix V that diagonalises M ν in terms of three mixing angles (ϑ e1 , ϑ e2 , ϑ 12 ) and three CP phases (δ, φ 1 , φ 2 ), and identify the regions of parameter space allowed by the consistency relations implied by (M ν ) 11 = 0. In section 4 we review the current upper limits and future sensitivities on the active-sterile mixing strength |V eN | 2 as a function of the sterile mass m N from 1 eV to 10 TeV. In section 5, we re-evaluate the 0νββ decay constraints in the generalised inverse seesaw, particularly for different values of the splitting ∆m N , and make a comparison with other constraints discussed earlier. We conclude in section 6. For completeness, the summary plots for constraints on |V µN | 2 and |V τ N | 2 are given in appendix A.
2 Generalised seesaw and neutrino mass spectrum

Model setup
We consider the addition of two SM-singlet Weyl fermions ν R,1 and ν R,2 to the SM particle content. We restrict ourselves to the first generation of SM fermions, which is the most relevant for 0νββ decay, and also allows us to present the gist of our results analytically. The SM Lagrangian is then extended to Here, L e = (ν L,e , e L ) T is the first-generation SM lepton doublet, H = iσ 2 H * with H = (H 0 , H − ) T being the SM Higgs doublet and σ 2 being the second Pauli matrix, M S is the Majorana mass term for the sterile states, and a summation over the sterile states is assumed (with i = 1, 2). After electroweak symmetry breaking by the vacuum expectation value H 0 = v 174 GeV, we obtain the neutrino Dirac mass terms (M D ) i = y ei v, and the Lagrangian (2.1) gives rise to the following neutrino mass matrix in the basis ν C L,e , ν R,1 , ν R,2 : The above mass matrix can be diagonalised by a 3 × 3 unitary matrix V such that giving rise to three mass eigenvalues m ν , m N 1 , m N 2 which can be chosen to be real and non-negative. We have denoted the mass eigenvalues suggestively for the case we will focus on, with one light, dominantly active, state (m ν 1 eV) and two much heavier,

JHEP03(2020)170
dominantly sterile, states (m N 1 , m N 2 m ν ). Accordingly, we conventionally order the mass eigenstates ν C L,e , ν R,1 , ν R,2 defined by The above minimal first-generation extension of the SM incorporates simplified versions of various seesaw scenarios: The type-I seesaw [3][4][5][6][7] is realised for ||M D || ||M S || (where ||M|| ≡ Tr(M † M) is the norm of matrix M). In fact, only one sterile state is minimally required to give mass to the one active neutrino considered here, i.e.
The light neutrino mass in this case is given by The minimal inverse seesaw [25][26][27] incorporates (M D ) 2 = 0 and (M S ) 11 = 0, so the neutrino mass matrix (2.2) becomes with µ S , m D m S . The light neutrino mass in this limit is given by The generalised inverse seesaw [33,34] incorporates (M D ) 2 = 0, but (M S ) 11 = µ R = 0, so the neutrino mass matrix (2.2) becomes with µ S , m D m S . This does not affect the mass of the light neutrino given by eq. (2.9) at tree-level, but will generate a one-loop correction [9,33] as discussed in section 2.2.

JHEP03(2020)170
The minimal linear seesaw [28][29][30][31] has (M D ) 2 = µ F = 0, but (M S ) 11 = (M S ) 22 = 0: with µ F , m D m S . The light neutrino mass in this case is given by Note that the mass matrix (2.11) can always be rotated to the form given by eq. (2.10) with appropriately defined µ R and µ S [32]. We will take advantage of this fact later to simplify our analysis, without loss of generality.
In the above scenarios we have not specified the source of LNV. Whether any one of the terms in eq. (2.2) violates lepton number will depend on the L assignment for the two sterile neutrinos ν R,1 , ν R,2 . For example, making the choice L(ν R,1 ) = L(ν R,2 ) = L(ν L,e ) = +1, suggested by treating the sterile neutrinos as RH counterparts to the LH active neutrinos, will mean that both terms in M D conserve L whereas all terms in M S violate L by two units. On the other hand, if L(ν R,1 ) = L(ν L,e ) = +1, L(ν R,2 ) = −1, the LNV terms are (M D ) 2 and (M S ) 12 = (M S ) 21 . While the choice of the origin of LNV is crucial to describe the underlying model, from a phenomenological point of view, the lepton number assignment does not need to be fixed. Also, any observable LNV effect crucially depends on the relative CP phase between the two sterile eigenstates, as we will see below. In any case, the smallness of the parameters µ R,S,F in the three seesaw variants discussed above is technically natural in the 't Hooft sense [43], i.e. in the limit of µ R,S,F → 0, lepton number symmetry is restored and the light neutrino ν L,e is exactly massless to all orders in perturbation theory, as in the SM.

Radiative corrections to the neutrino mass
The light neutrino mass acquires a one-loop radiative correction from the self-energy diagrams involving the SM gauge and Higgs bosons [9,44,45], induced by the Lagrangian (2.1). In terms of the 1 × 2 matrix M D and the 2 × 2 matrix M S as defined through eq. (2.2), the finite loop contribution in our single-generation case can be written as [33] δm 1-loop JHEP03(2020)170 In our analysis, we will require that the one-loop corrections are subdominant to the tree-level mass, using a 10% contribution as the limit, (2.15) Using different loop-to-tree contribution ratios will not change our results qualitatively.
3 Phenomenological parametrisation of the mixing matrix As noted before, we will neglect the flavour structure of the lepton sector and work in a single-generation picture with only an electron flavour active neutrino field and two sterile fields; ν C L,e , ν R,1 and ν R,2 . In this case the general neutrino mass matrix M ν in eq. (2.2) can be diagonalised by a 3 × 3 unitary matrix V as described in eq. (2.3). It is simple to reverse this diagonalisation in order to express the mass matrix in terms of the a priori measurable mixing angles, CP phases and mass eigenvalues, We can first consider a parametrisation of V analogous to that of the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) mixing matrix accompanying charged currents in the SM, in terms of the cosine c ij ≡ cos ϑ ij and sine s ij ≡ sin ϑ ij of the three mixing angles ϑ e1 , ϑ e2 and ϑ 12 . They describe, respectively, the mixing between the mostly-active light neutrino mass eigenstate ν e and the first mostly-sterile mass eigenstate N 1 , ν e and the second mostlysterile mass eigenstate N 2 and finally between N 1 and N 2 . The angles can in principle lie in the range ϑ ij ∈ [0, π/2] and the equivalent of the Dirac CP phase in the range δ ∈ [0, 2π]. D is a diagonal matrix containing the remaining two Majorana phases φ 1,2 ∈ [0, 2π], As for the light active neutrino PMNS mixing matrix, only two physical Majorana phases survive because an overall phase can be rotated away.

JHEP03(2020)170
Rather than this phenomenological approach we can instead write V in a form explicitly imposing existing constraints from neutrino oscillations. A convenient way to do this is the so-called Casas-Ibarra parametrisation [46], which has been generalised in ref. [47] to include the complete parameter space of sterile neutrino masses and mixings. Here, in the three-generation picture and for two sterile states the active-sterile mixings are related to the light active neutrino masses m i (assuming m 1 = 0), heavy neutrino masses m N i and PMNS mixing matrix elements by where R is an arbitrary 2 × 2 orthogonal matrix parametrised by a complex mixing angle ϑ 45 + iγ 45 and H is a hermitian matrix encoding deviations from unitarity in the light neutrino sector. For fixed values of m i and m N i the size of mixings V eN i depend on ϑ 45 and γ 45 . In the phenomenological single-generation picture this translates to choices of the CP phases φ 1 , φ 2 and δ. We will proceed with our phenomenological approach because it is not our immediate goal to reproduce the observed light neutrino data, which is an implicit input to the Casas-Ibarra parametrisation. Our goal is to investigate in the most direct way the phenomenology of active-sterile mixing in the generalised inverse seesaw.

Consistency relations
We will now apply this general parametrisation to the seesaw scenarios discussed in section 2.1. Without a triplet Higgs extending the SM field content, the active neutrinos cannot acquire a mass of the formν C L ν L and thus the (1, 1) entry of M ν in eq. (2.2) is strictly zero at tree-level. This requirement must be satisfied irrespective of the remaining mass matrix structure (i.e. type-I, inverse or linear seesaw). Written in terms of the phenomenological parameters this condition may be written as where we have divided the sum by the heavy neutrino mass m N 1 . We note first that this constraint has no dependence on the sterile-sterile mixing angle ϑ 12 . It can also be seen that such a constraint is equivalent to the vanishing of the effective 0νββ decay mass m ββ = i (U PMNS ) 2 ei m i , where in that case the summation is over the three light neutrino mass eigenstates. While this would be an accidental cancellation -possible for a normally ordered light neutrino spectrum with specific values of the Majorana phases in U PMNS (as opposed to V ) -the condition in eq. (3.5) must always be satisfied at tree-level, putting requirements on the values of the three masses, three mixing angles and three CP phases. Instead of the parameter m N 2 it is equally valid to use the mass splitting ∆m N = m N 2 − m N 1 , which will be of importance later.
As illustrated in figure 1, the condition in eq. (3.5) can be visualised as a triangle in the complex plane, formed by three sides with lengths The angles between these sides are determined by the phase φ 1 and the linear combination φ 2 − 2δ, which we label φ 2 for JHEP03(2020)170 convenience. Not all combinations of the masses and mixings allow a triangle to be formed with side lengths L ν , L 1 and L 2 . Specifically, the triangle can only be closed (for some values of φ 1 and φ 2 ) if the longest length is smaller (or equal) to the sum of the shorter lengths, max(L ν , L 1 , L 2 ) ≤ min(L ν , L 1 , L 2 ) + med(L ν , L 1 , L 2 ) . The allowed regions for the squared active-sterile mixing strengths s 2 e1 and s 2 e2 are shown in figure 2 (left) for different choices of the light and heavy neutrino masses. The centre shape (light blue) corresponds to the choice m ν /m N 1 = 10 −10 and ∆m N /m N 1 = 10 −2 . This for example could correspond to a light neutrino mass m ν = 10 −3 eV and heavy neutrino masses m N 1 = 10 MeV and m N 2 = 10.01 MeV. The allowed mixing strengths form a region centred around s 2 e1 ≈ s 2 e2 ≈ m ν /m N 1 . Thin, virtually line-like extensions to large s 2 e1 = s 2 e2 , small s 2 e1 and small s 2 e2 are also possible. As can be seen from the dark blue and green regions, increasing (decreasing) r ν = m ν /m N 1 will move the bulk of the region along the diagonal to higher (smaller) mixing. As can be seen from the yellow region, increasing the splitting ∆m N shifts the allowed region to smaller values of s 2 e2 but not s 2 e1 . The red region, on the other hand, shows the scenario in which ∆m N becomes negative (when m N 2 < m N 1 ). The allowed region instead moves up to larger s 2 e2 for the same s 2 e1 . We will investigate this behaviour more quantitatively below. Figure 2 (right) shows the same regions but with the axes given by the ratio and sum of the mixing strengths, s 2 e2 /s 2 e1 and s 2 e1 + s 2 e2 , respectively. It especially illustrates that there exists a lower limit on the total active-sterile mixing strength s 2 e1 + s

CP phases
The limiting behaviours for small and large mixing strengths can be related to the CPconserving cases when the phases adopt values such that e iφ 1 = ±1, e iφ 2 = ±1, corresponding to the relative CP parity of the sterile fields. The CP parity of the m ν state is defined by convention as +1. Three possibilities emerge for the CP parities of the other states: Because s 2 e2 ≥ 0 this can only be satisfied if s 2 e1 ≤ m ν /(m N 1 + m ν ) m ν /m N 1 , i.e. for s 2 e1 up to the ordinary single heavy-state seesaw mixing s 2 e1 = m ν /m N 1 . Consequently, s 2 e2 can range from s 2 e2 = 0 (when s 2 e1 = m ν /m N 1 ) to s 2 e2 ≈ m ν /m N 2 (when s 2 e1 = 0). This scenario corresponds to the canonical seesaw with two heavy Majorana states; the active state can mix with either of them with adjustable strength. In figure 2 (left), this particular limit corresponds to the line-like extensions towards vanishing s 2 e2 at the bottom (N 2 decouples, s 2 e1 → m ν /m N 1 ) and vanishing s 2 e1 to the JHEP03(2020)170 . Intermediate solutions lie on the lower left edge of the allowed region in figure 2 (left). Rearranging eq. (3.7) for small ∆m N gives s 2 e1 + s 2 e2 ≈ m ν /m N 1 , a behaviour that can clearly be seen in figure 2 (right).
(C) e iφ 1 = +1, e iφ 2 = −1: The contributions of the heavy states can (partially) cancel among each other, L ν + (L 1 − L 2 ) = 0. Again, we can solve for the mixing angle s 2 e2 , Here, no upper bound on s 2 e1 exists and it can in principle take values between 0 ≤ s 2 e1 ≤ 1. For a small mass splitting ∆m N m N 1 this case corresponds to the inverse seesaw scenario where the two heavy Majorana states form a pseudo-Dirac neutrino pair. In figure 2 (left), this limit corresponds to the thin extension of the allowed region to large mixing strengths. It should be noted that this phenomenological parametrisation does not enforce a small mass splitting and ∆m N can be arbitrarily large for a given light neutrino mass m ν . As we will discuss below, however, this will induce large loop corrections to m ν .
For arbitrary values of the phases φ 1 and φ 2 the interior of the shaded regions in figure 2 is covered. In order to simplify the following discussion we make use of the dimensionless ratios r ν = m ν /m N 1 and r ∆ = ∆m N /m N 1 as already introduced in figure 2. For arbitrary phases, eq. (3.5) in fact represents two conditions; Re{(M ν ) 11 } = 0 and Im{(M ν ) 11 } = 0. These relations can be rearranged to find two equivalent expressions for s 2 e2 , where the first and second equalities are derived from the real and imaginary conditions, respectively. We can also rearrange eq. (3.9) to solve for the tangent of φ 2 , where we also indicate approximate solutions for the different limits of s 2 e1 . In effect, the condition in eq. (3.5) has allowed us to eliminate two parameters, s 2 e2 and φ 2 , by expressing them in terms of a subset of the remaining free parameters, r ν , r ∆ , s 2 e1 and φ 1 . The freedom to divide eq. (3.5) by m N 1 and using instead the ratios r ν and r ∆ also effectively removes a mass degree of freedom. This can be seen from the behaviour of the allowed regions in figure 2; a shift in the s 2 e1 -s 2 e2 plane only occurs when r ν and r ∆ are changed. It must however be remembered that the other elements of M ν (e.g. m D , m S , µ R , µ S ) have been divided by m N 1 , so this factor must be taken into account when calculating these flavour-basis parameters as functions of the phenomenological mass-basis parameters.

JHEP03(2020)170
Alternatively a more physical choice would be to solve for cos φ 1 and cos φ 2 using the cosine rule for the (M ν ) 11 = 0 constraint triangle in figure 1, where the approximate expressions hold for small mixing s 2 e1 , s 2 e2 1. In this way the phases φ 1 and φ 2 are determined (up to a pair of solutions in the range [0, 2π], modulo π) by the neutrino masses through the ratios r ν and r ∆ and the mixing strengths s 2 e1 and s 2 e2 , all of which are in principle experimentally measurable. If the solution for φ 1 lies in the first or second quadrant (i.e. φ 1 ∈ [0, π]), in order to close the triangle in figure 1 it is necessary for φ 2 to be in the third or fourth quadrants (φ 2 ∈ [π, 2π]) and vice versa.
An important parameter in determining the nature of the two heavy states is the phase difference ∆φ = φ 1 − φ 2 = φ 1 − φ 2 + 2δ between N 1 and N 2 . If ∆φ ≈ 0 we expect the heavy states to behave like Majorana fermions, whereas for ∆φ ≈ ±π they should form a pseudo-Dirac pair with an associated suppression of LNV effects. Using the solutions eqs. (3.11) and (3.12), or alternatively using the cosine rule for the third angle of the triangle in figure 1, ∆φ is given in terms of the other parameters by This phase difference is plotted in figure 3 (left) as a function of the mixing strengths s 2 e1 and s 2 e2 within the region allowed by the (1, 1) element constraint eq. (3.5). Note that the active-sterile mixing strengths s 2 e1 and s 2 e2 are normalised by r ν and r ν /(1 + r ν + r ∆ ) respectively, making the plot generically applicable for an arbitrary choice of the light and heavy neutrino masses. The edges of the allowed region correspond to the CP-conserving combinations of phases: (i) φ 1 = φ 2 = π to the lower left corresponding to the canonical seesaw with two Majorana heavy states and (ii) φ 1 = 0 (π), φ 2 = π (2π) on the top (lower right) edge, corresponding to an inverse seesaw -like scenario. Intermediate scenarios interpolating between these limiting cases are characterised by the phase difference |∆φ| increasing from 0 to π as shown.
We have so far seen that is it possible to eliminate the two phases φ 1 and φ 2 from the nine initial phenomenological parameters. The next question is whether additional relationships can be found between these parameters. While eq. (3.5) is a parametrisationindependent condition, we can make convenient choices for the remaining parameters in M ν which assist in this effort. For example, as discussed briefly in section 2.1, without lack of generality we can assume that the (1, 3) element of the neutrino mass matrix M ν in eq. (2.2) vanishes, (M ν ) 13 = 0. This can always be achieved by rotating the heavy states appropriately, even in the linear seesaw scenario [32]. Using our phenomenological JHEP03(2020)170  parametrisation this corresponds to Note that in this condition the linear combination φ 2 = φ 2 − 2δ does not appear explicitly. As we would like to continue using the relations for cos φ 1 and cos φ 2 in eqs. (3.11) and (3.12) we introduce the linear combination δ = 2φ 2 + δ orthogonal to φ 2 . The phases φ 1 , φ 2 and δ can consequently be written as linear combinations of φ 1 , φ 2 and δ . As the (1, 1) element constraint we can take both the real and imaginary part of eq. (3.14), rearranging for s 2 12 as a function of s 2 e1 , s 2 e2 and the phases, (3.16) The sterile-sterile mixing strength s 2 12 is shown in figure 3 (right) as a function of s 2 e1 and s 2 e2 for δ = 0. Further, proceeding as before, we can equate the real and imaginary solutions of JHEP03(2020)170 s 2 12 in eq. (3.15), i.e. C R = C I . Rewriting in terms of the phases φ 1 , φ 2 and δ and making use of the solutions for cos φ 1 and cos φ 2 in theory allows to solve for the final phase δ in terms of r ν , r ∆ , s 2 e1 and s 2 e2 . In practice it is difficult to do this analytically, but numerically δ can be found by finding the intersecting points of the curves C R (δ ) and C I (δ ).
We have therefore seen that, given values of the parameters r ν , r ∆ , s 2 e1 and s 2 e2 and assuming a particular parametrisation of the neutrino mass matrix M ν , the remaining parameters s 2 12 , φ 1 , φ 2 and δ are uniquely determined. Thus, if the absolute neutrino mass scale m ν were known and an experiment were to observe two sterile states with mass splitting ∆m N and mixing strengths s 2 e1 and s 2 e2 , in the generalised inverse seesaw the sterile-sterile mixing strength s 2 12 and CP phases φ 1 , φ 2 and δ are predicted quantities. As we will see in section 4, direct searches for the production and decay of heavy states can probe, if not sensitive to the lepton numbers of the final states, the mixing matrix for particular values of m N 1 or m N 2 . If the splitting ∆m N is large enough for the two states to be resolved, |V eN 1 | 2 and |V eN 2 | 2 can be measured independently, constraining the values of the other parameters. If the splitting is below the energy resolution of an experiment it will instead be sensitive to the sum |V eN 1 | 2 + |V eN 2 | 2 . As seen in figure 3 (left), this can only put a lower bound on ∆φ while s 2 12 and δ are left unconstrained. Most current and future direct searches are still probing the regime |V eN 1 | 2 ≈ |V eN 2 | 2 r ν , where the generalised inverse seesaw predicts the phase difference ∆φ = ±π. Some experiments like the KATRIN upgrade TRISTAN [48] and the future long-baseline neutrino oscillation experiment DUNE [49] may however reach mixing strengths |V eN 1 | 2 r ν , thus being able to pin down any phase difference in the range |∆φ| ∈ [0, π], cf. figure 7.
The next question to ask is whether the parameters s 2 12 , φ 1 , φ 2 and δ can be measured in order to confirm the predictions of the generalised inverse seesaw. The Majorana and pseudo-Dirac limits (governed by φ 1 and φ 2 ) are primarily distinguished by the magnitude of LNV. In the case where the sterile mass splitting is not too small, LNV searches are currently probing mixing strengths in the pseudo-Dirac limit. It is unlikely for future LNV searches to be able to reach the mixing strengths |V eN 1 | 2 r ν required for the Majorana limit. Put differently, if an experiment sees two sterile states with mixings |V eN 1 | 2 ≈ |V eN 2 | 2 r ν , but also a large LNV signal (e.g. from a large asymmetry in the pseudorapidity distribution at the ILC [50]), it would strongly imply some other source of LNV [20]. For example, the states N 1 and N 2 could possess additional strong couplings to SM particles from a TeV-scale type-III seesaw mechanism, or the light neutrino masses are not generated by the seesaw (e.g. instead, radiatively) [18].
We next consider s 2 12 . In the small mixing limit s 2 e1 , s 2 e2 1, the matrix in the basis that the charged lepton Yukawa coupling is diagonal. In ref. [38] it was noted that the Dirac sub-matrix M D can always be redefined as measure the angle ϑ 12 , making it unphysical (see also ref. [51]). If right-handed currents are introduced, for example in left-right symmetric models, s 2 12 in theory becomes observable because the lower two sub-matrices of V in eq. (3.2) rotate the W R gauge boson interaction. In the pseudo-Dirac case it becomes possible to observe heavy neutrino mixing via the ratio of same-sign to opposite-sign charged lepton production rates in colliders [37,52,53], where Γ N is the average decay width of the sterile neutrinos. The distinguishing signal here is that R can take an intermediate value between 0 (Dirac limit) and 1 (Majorana limit). While the elements of the mixing matrix V N containing s 2 12 appear in the same-sign and opposite-sign rates, they cancel in the numerator and denomator for ∆φ = ±π. This is generally not true if |∆φ| < π.

Including loop corrections
The mixing strength s 2 12 is nonetheless important for evaluating the radiatively generated neutrino mass at one-loop in eq. (2.13) (exact expression) and eq. (2.14) (in the limit µ R,S m S ). When written in terms of the masses, mixing angles and CP phases (in the particular parametrisation setting the (1, 3) element of M ν to zero), the flavour-space parameters m D , m S , µ S and µ R are functions of s 2 12 . In evaluating these parameters for the purpose of evaluating δm 1-loop ν , we will for simplicity assume δ = 0 from the start instead of numerically solving C R = C I for s 2 12 and δ for given values of m ν , m N 1 , r ∆ , s 2 e1 and s 2 e2 . We reiterate that m ν and m N 1 must be chosen independently (instead of just the ratio r ν ) because an overall factor m N 1 cannot be eliminated from m D , m S , µ S and µ R as for the (M ν ) 11 = 0 and (M ν ) 13 = 0 constraints.
In this scenario we can investigate the value of ϑ 12 for the limiting cases of φ 1 and φ 2 = φ 2 along the edges of the allowed region in figure 3. Applying the limits s 2 e1 , s 2 e2 1 and m ν m N 1 to the expression for s 2 12 in eq. (3.15), the cases resolve to: (A) e iφ 1 = e iφ 2 = +1: No solution.
(B) e iφ 1 = e iφ 2 = −1: In this case we have where s 2 e1 ≤ r ν as discussed before in this case, making the root well defined.
The general behaviour for s 2 12 is shown in figure 3 (right) as a function of the active-sterile mixing strengths s 2 e1 and s 2 e2 . At each point in the allowed region the phases φ 1 and φ 2 JHEP03(2020)170 are calculated according to eqs. (3.11) and (3.12) as shown in figure 3 (left), while δ is set to zero. We see that the sterile-sterile mixing is ϑ 12 = π/2 when s 2 e1 r ν . As s 2 e1 approaches r ν along the canonical seesaw side of the allowed region the mixing angle falls to ϑ 12 = 0. These two values are physically equivalent, signifying an exchange in the role of the two heavy states as one state becomes decoupled while the other state's mixing strength increases to r ν or r ν /(1 + r ν + r ∆ ). In the inverse seesaw limit the sterile-sterile mixing angle approaches ϑ 12 = π/4, i.e. maximal mixing.
With the sterile-sterile mixing strength s 2 12 taken care of, we now return to the neutrino mass generated at one-loop. So far in this section we have worked at tree-level. From gauge invariance of the SM Lagrangian under SU(2) L , it is not possible to write a Majorana mass termν C L ν L for the left-handed neutrino field, and thus the (1, 1) element of the neutrino mass matrix is zero. The inclusion of loop corrections will to first loop-order however lead to the appearance of a finite value for the (1, 1) element in eq. (2.10), where δm 1-loop ν is given by eq. (2.13). This will contribute to the mass eigenvalue of the lightest state as is the tree-level mass from the diagonalisation of the mass matrix eq. (2.10) as discussed in section 2.1. When using m ν from now on we assume that it is the physical mass as measured by an experiment, including both the tree-level and one-loop contributions.
In figure 4 (left), we plot the exact formula for δm 1-loop ν in eq. (2.13) as a function of the heavy neutrino mass m N 1 and the mixing strength s 2 e1 . The parameters m ν , r ∆ , φ 1 and φ 2 (for δ = 0) are fixed as indicated in the figure, while s 2 e2 and s 2 12 are calculated according to eqs. (3.9) and (3.15), respectively. Specifically the tree-level mass and the relative heavy neutrino splitting are given for the benchmark values m ν = 10 −3 eV and r ∆ = 10 −2 , while the Majorana phases are chosen such that the scenario is located on the right edge of the allowed parameter space in figure 3 (left). We also plot the 'seesaw' line s 2 e1 = r ν = m ν /m N 1 in grey. Below this line s 2 e2 will tend to the constant value r ν /(1 + r ν + r ∆ ) ≈ r ν , while s 2 12 tends towards π/2. Above this line is the inverse seesaw limit with s 2 e2 = s 2 e1 /(1 + r ∆ ) ≈ s 2 e1 and s 2 12 = π/4. This plot demonstrates the strong dependence of |δm 1-loop ν | on the model parameters. For large m N 1 , we can already see that the one-loop corrections are dangerously large as a consequence of the comparatively large splitting between the heavy states ∆m N = r ∆ m N 1 . Looking at the approximate loop formula in eq. (2.14) and recalling that m D , m S , µ R contain terms proportional to m N 1 (when written in terms of the mass-basis parameters and mixing angles), the strong dependence on m N 1 is not surprising because δm 1-loop ν naively scales as m 3  . Solid lines are found by using the exact formula eq. (2.13), while the dashed lines use this same formula but in the limit µ R,S m S , given by ref. [33].
As stated before, in this work we will maintain the assumption that the loop corrections to the light neutrino mass are sub-dominant, i.e. we assume that the neutrinos largely acquire their masses via the tree-level seesaw mechanism. A reasonable requirement that |δm 1-loop ν | < 0.1m ν [cf. eq. (2.15)] can subsequently be used to set an upper limit on the active-sterile mixing strengths. This is shown in figure 4 (right) as a function of the heavy neutrino mass m N 1 for different values of r ∆ = ∆m N /m N 1 . It can be seen that as the relative splitting r ∆ becomes smaller, the associated upper limit on the mixing strength becomes weaker. The solid and dashed lines correspond to the upper limit derived from the exact formula eq. (2.13) and the approximation eq. (2.14), respectively. It can be seen that the exact and approximate upper limits diverge for small m N 1 and s 2 e1 -this is because µ R,S m S no longer holds in this particular region of the parameter space.
In figure 5 we again plot the region satisfying the tree-level constraint (M ν ) 11 = 0 in eq. (3.5), but now also exclude the region not satisfying the |δm 1-loop ν | < 0.1m ν loop requirement for different values of the relative splitting r ∆ . It can be seen that as r ∆ increases the allowed region is reduced, excluding much of the inverse seesaw region. It is worth mentioning that in order to see this effect around s 2 e1 ∼ r ν requires large relative splittings, otherwise the loop requirement only excludes much larger mixings strengths s 2 e1 ≈ s 2 e2 in the inverse seesaw limit. While combining the constraints (M ν ) 11 = 0 and where we take the neutrino mass on the l.h.s. to be the tree-level mass to first approximation. Writing m tree ν = m ν − δm 1-loop ν via eq. (3.22), eq. (3.23) can be rearranged as before to solve for s 2 e2 and cos φ 2 , but now as a function of the loop mass. Paradoxically, s 2 e2 and φ 2 are themselves required to evaluate the loop mass in eq. (2.13) as a function of m N 1 and s 2 e1 . Inserting the new expressions for s 2 e2 and φ 2 , the loop mass can be evaluated iteratively by first setting (M ν ) 11 = 0 and then re-inserting each new value back into the one-loop formula. We find that the difference between the initial (setting (M ν ) 11 = 0) and iterated loop mass is negligibly small when the initial loop mass satisfies |δm 1-loop ν | < 0.1m ν . When the initial loop mass is larger this iterative approach is strictly no longer valid, but we assume that is viable up to |δm 1-loop ν | ∼ 0.1m ν . This should then not significantly affect the upper bounds on s 2 e1 derived from the loop condition. In other words, we keep the constraints derived using (M ν ) 11 ≈ 0 and |δm 1-loop ν | < 0.1m ν .

Constraints on heavy sterile neutrinos
In this section we will summarise the results of experimental searches for sterile neutrinos and hence constraints on the active-sterile mixing |V N | 2 over the sterile neutrino mass range m N ∈ [1 eV, 10 TeV]. For lighter masses it becomes possible for one of the sterile states to form a quasi-Dirac state with the active state. A large portion of this parameter space is constrained by solar neutrino oscillations [54,55]. For heavier masses m N 10 TeV,

JHEP03(2020)170
sterile neutrinos can generate the light active neutrino masses via the conventional seesaw mechanism. These neutrinos, however, are not kinematically accessible to direct searches. The constraints from existing searches and observations in the m N −|V eN | 2 parameter space are shown in figure 6 by various shaded regions, whereas figure 7 illustrates the sensitivity of expected future experiments and observations. As our ultimate focus is on a comparison with constraints from 0νββ decay in section 5, we focus on the first generation mixing element |V eN | 2 . However, for the sake of completion and future reference, we also compile and update the constraints on |V µN | 2 and |V τ N | 2 in appendix A. For earlier summary plots showing a partial list of these constraints, see e.g. refs. [39][40][41][42]. Most limits shown in the plots were derived assuming a single heavy neutrino. For small splitting, the limit is applicable on |V eN 1 | 2 +|V eN 2 | 2 ; for large splitting the limits apply separately for each species.

High-energy collider searches
Heavy states are produced in charged-current (CC) and neutral-current (NC) processes through their admixture with the active states, and thus their decay products can be searched for at high-energy colliders which copiously produce W and Z bosons. For sufficiently small mixing angles, the macroscopic decay-ength of the heavy neutrinos can result in displaced vertices with distinct detector signatures. We consider the following searches (keywords in bold match the corresponding regions in figures 6 and 7): • The LHC collaborations ATLAS and CMS have searched for N production and decay through various channels. Both have recently searched for decays of Wproduced N to three charged leptons, W ± → ± N, N → ± ∓ ν ( = e, µ), either in the LNC or LNV mode. ATLAS used the prompt final state of three isolated leptons and no opposite-charge same-flavour lepton pairs (LNV channel) to reject Drell-Yan, W + jets and tt backgrounds. CMS broadened the search to the LNC channel with a sensitivity to displaced decays. The analyses impose the limits |V eN | 2 , |V µN | 2 < 10 −5 − 10 −4 over the mass range 5 GeV < m N < 50 GeV [56,57]. ATLAS and CMS have also conducted searches for the LNV same-sign dilepton + jets channel, W ± → ± N, N → ± jj [58,59]. Above the Z boson mass limits can be improved in future by ATLAS and CMS during the high luminosity (L = 3 ab −1 ) LHC phase (HL-LHC) and by a future √ s = 27 or 100 TeV Future Circular Collider (FCC-hh) [60,61]. Around the Higgs mass, limits can also be set from the SM Higgs decay to sterile neutrinos [62].
• In the future, ATLAS, CMS and LHCb are expected to probe smaller |V N | 2 through displaced vertex searches. For a given mixing, m N must lie in a specific range in order to avoid N decaying promptly or outside the detector. The best projected limit is |V eN | 2 , |V µN | 2 10 −9 for m N ≈ 30 GeV [56].
• At the LEP collider, the collaborations L3 [63,64] and DELPHI [65] searched for N produced through on-shell Z production, e + e − → Z → N ν , followed by the decays N → ∓ W ± , N → ν Z and N → ν H. Using N → e ∓ W ± and W ± → jj, L3 enforced a limit of |V eN | 2 < 10 −4 in the range 5 GeV < m N < 80 GeV. This was reduced to JHEP03(2020)170 |V eN | 2 < 10 −5 by an improved DELPHI analysis. At a future linear electron-electron collider such as the ILC [66], for a benchmark √ s = 500 GeV and L = 100 fb −1 limits may be improved to |V eN | 2 < 10 −4 above the Z mass. At a proposed Compact Linear Collider (CLIC), for √ s = 3 TeV and L = 1 ab −1 limits are |V eN | 2 10 −5 − 10 −4 for 600 GeV < m N < 2.3 TeV [67,68]. Furthermore, a future FCC-ee collider, acting as a powerful e + e − → Z factory and exploiting low backgrounds in displaced vertex searches, can improve the sensitivity drastically; down to |V eN | 2 10 −11 for m N ≈ 50 GeV [69]. At the ILC it may also be possible to distinguish LNC and LNV W ± exchange channels between the e + e − pair by measuring the asymmetry of the outgoing lepton pseudorapidity distribution [50]. Finally, the proposed Large Hadron-Electron Collider (LHeC) LHC upgrade may also provide competitive constraints for m N > m Z [68,70,71]. An overview of proposed collider sensitivities is given in ref. [72].

On the LNV signal at colliders
As for the LNV signature at colliders, in a natural seesaw scenario with approximate lepton number conservation, the LNV amplitude for the on-shell production of heavy neutrinos at average four-momentum squareds = (m 2 N 1 + m 2 N 2 )/2 can be written as [35,80] A for ∆m N Γ N , i.e. for a small mass splitting |∆m N | = |m N 2 − m N 1 | between the heavy neutrinos compared to their average decay width Γ N ≡ (Γ N 1 + Γ N 2 )/2. Thus, the LNV amplitude in eq. (4.1) will be suppressed by the small mass splitting, except for the case ∆m N Γ N when it can be resonantly enhanced [35,81].
For the 5 − 50 GeV range of sterile neutrino masses probed by the ATLAS and CMS same-sign trilepton and dilepton + jets analyses, the total sterile neutrino decay width, if decays only takes to place to SM leptonic and hadronic degrees of freedom, is given by

JHEP03(2020)170
where the complete expressions for the factors a (m N ) are given in refs. [39,82]. The factors a (m N ) include the contributions from two-body semi-leptonic and three-body leptonic decays, and are approximately given by where N 2−body and N 3−body are the number of decay channels open for each decay topology. Γ 2−body and Γ 3−body are given roughly by where f M is the order of magnitude of the meson decay constants [39]. For m N ≈ 50 GeV all three-body leptonic decays and two-body semi-leptonic decays to pseudoscalar mesons (π 0 , η, η , η c , η b , π ± , K ± , D ± , D ± s , B ± , B ± c ) and vector mesons (ρ 0 , ω, φ, J/ψ, Υ(4S), are open, and so the total decay width (for |V µN | 2 = |V τ N | 2 = 0) is approximately For small splittings, e.g. r ∆ = 10 −4 and hence ∆m N ≈ 5 MeV for m N ≈ 50 GeV, and the |V eN | 2 ∼ 10 −5 mixing probed by the LNV analyses, eq. (4.5) implies that Γ N /∆m N ∼ 10 −6 . Collider searches specifically looking for an LNV signal in figure 7 are therefore still valid for this splitting and splittings down to r ∆ ∼ 10 −10 . As will be discussed later, this is important for the comparison with 0νββ decay in this mass range. We finally note that the analysis of ref. [83] gives an estimate for the regions of the m N − |V N | 2 parameter space where the ratio R in eq. (3.18) is less than or greater than a third. Comparing with figure 1 of that work, we again confirm that LNV signals searched for by colliders below the electroweak scale remain unsuppressed, particularly for ∆m N of order the light neutrino mass splittings (motivated by naturalness).

Meson decays and beam-dump experiments
At the intensity frontier N can be produced abundantly in beam-dump experiments and through various meson decays. We consider the following limits: • The TRIUMF PIENU experiment [84] conducted a search for N produced in pion decays at rest. Utilising the helicity suppression of the π → eν decay channel in comparison to π → µν channel, the presence of N induces extra peaks in the lower positron energy region. Improving on previous results limited by the background µ + → e + ν eνµ , the collaboration set limits at the level of |V N | 2 10 −8 in the range 60 MeV < m N < 129 MeV [85][86][87].
• The NA62 experiment [88] used a secondary 75 GeV hadron beam containing a fraction of kaons, and has been able to probe the decays K + → + N ( = e, µ). For small |V N | 2 the N decay length is much longer than the 156 m detector volume and the process is characterised by a single detected track -a positive signal is JHEP03(2020)170 a peak in the missing mass distribution. Limits |V eN | 2 , |V µN | 2 < 10 −8 − 10 −9 in the range 170 MeV < m N < 450 MeV (up to the kaon mass) have been made. In future NA62 will be converted to a beam-dump configuration and will be able to probe hadronic decays to N , followed by N decays, up to the D meson mass. The projected sensitivity is |V eN | 2 , |V µN | 2 < 10 −8 for 1 GeV < m N < 2 GeV [89]. A recent recalculation of the impact of sterile neutrinos on kaon decays was conducted in ref. [90].
• The Belle experiment [91] was a B factory that extended the peak search method to higher energies -using BB pairs collected at the Υ(4S) resonance, the decay mode B → (X) N , with X a charmed meson D ( * ) or light meson, could be followed by N → π ( = e, µ). Constraints were made between the K and B meson masses and at best were |V N | 2 3 × 10 −5 for m N ≈ 2 GeV [92].
• The NA3 experiment [93] collided a secondary 300 GeV π − beam with an iron absorber, producing hadronic states which subsequently decayed to leptonic, semileptonic or fully hadronic final states. N decays producing leptonic or semi-leptonic final states could be produced from the decays of π, K, D and B mesons. NA3 was most sensitive up to the D meson mass, setting limits of |V eN | 2 , |V µN | 2 < 10 −4 for 1 GeV < m N < 2 GeV.
• Accelerated neutrino beam experiments have conducted a variety of parallel searches. The CHARM [94,95] and PS191 [96] experiments and the IHEP-JINR neutrino detector [85,97] searched for a small fraction of N in a predominantly ν µ beam. The beams were produced by colliding a primary beam of protons with an iron or copper fixed target, with the hadronic products decaying as π/K/D → ν(N ) ( = e, µ). If sufficiently massive, N may decay before reaching the detector via the channel N → + − ν . CHARM also used a wide-band neutrino beam to constrain the NC process ν µ n(p) → N X followed by N → µX within the detector. IHEP-JINR and PS191 provide constraints (down to |V eN | 2 10 −7 and |V eN | 2 10 −9 , respectively) up to the kaon mass. CHARM provides constraints up to the D meson mass, at best |V eN | 2 , |V µN | 2 10 −7 for 1 GeV < m N < 2 GeV.
• The long-baseline neutrino oscillation experiment T2K [98] searched for an admixture of N in its initial neutrino beam flux, produced by colliding 30 GeV protons with a graphite target at J-PARC. Daughter K ± of a given charge are focused and decay via K → ν(N ). The off-axis near-detector at a baseline of 280 m searched for N decays via the channel N → π, improving on the constraints made by PS191. In future, the near detector of the oscillation experiment DUNE will be highly sensitive for m N up to the D s meson mass [99,100].
• The future beam-dump experiment SHiP [101] is purposely designed to look for exotic long-lived particles. Utilising a 400 GeV proton beam from the CERN Super Proton Synchrotron, it is expected to be sensitive to sterile neutrinos with m N up to the B c meson mass (∼ 6 GeV). In a benchmark scenario where the electron-sterile JHEP03(2020)170 coupling dominates, SHiP is expected to be sensitive down to |V eN | 2 10 −10 for m N ≈ 1.6 GeV [102].
• In parallel with collider searches it is possible to look for LNV Decays of tau leptons and pseudoscalar mesons as discussed in refs. [39,82,103,104]. One issue is that if the LNV process is mediated by the light neutrinos the amplitude is proportional to and suppressed by the small m 2 ν , while if mediated by heavy neutrinos it is suppressed by 1/m N and |V N | 2 . LNV decay widths however can be strongly enhanced if a sterile state is produced on-shell. The sensitivity of NA62 to three-body LNV light mesons decays (K + → + + π − ), BESIII to charmed meson decays (D + /D + s → + + π − /K − ) and BaBar, Belle and LHCb for B meson decays (B + → + + π − /K − /D − /ρ − /K * − ) for , = e, µ were estimated most recently by ref. [104]. The BESIII has also conducted its own analysis on the (D + → + + π − /K − ) decay channel [105]. Finally, the Future LNV decay sensitivities of NA62, LHCb, Belle-II, MATHUSLA, SHiP and FCC-ee have been explored in ref. [106]

Beta decays and nuclear processes
Active neutrinos are produced in the β-decays of unstable isotopes and in nuclear fission processes. Heavy sterile neutrinos can also be produced via the active-sterile mixing if the sterile mass is smaller than the energy release (Q-value) of the relevant nuclear process. The production of a sterile state results in a distortion or 'kink' in the β-decay spectrum and associated Kurie plot. It is also possible for the sterile state to decay before detection. We include the following searches: • Heavy neutrinos produced in β-decays significantly alter the energy spectrum of the emitted β electron. In order to be kinematically accessible the sterile neutrino mass must be smaller than the Q-value of the process, m N < Q β . If this is satisfied and the sterile states are sufficiently more massive than the active states, the β-decay spectrum becomes the incoherent sum where m 2 β = k |U ek | 2 m 2 k is the usual scale probed by β-decay [107]. This expression can give rise to multiple kinks in the spectrum of relative size |V eN i | 2 and at energies E kink = Q β − m N i . Such an effect for a single sterile neutrino has been probed for a variety of isotopes with a range of different Q-values, and therefore sensitive to different m N . Isotopes include 3 [118]. In the future, strongly improved limits by the operating tritium β-decay experiment KATRIN and the proposed TRISTAN upgrade are expected [48]. The capability of the PROJECT 8 experiment, which uses the alternative method of cyclotron radiation emission spectroscopy, has also been briefly explored [119].

JHEP03(2020)170
• Reactor neutrino experiments are sensitive to sterile neutrinos with masses in the range 1 MeV < m N < 10 MeV. At these masses it is possible for N to decay within the detector via the channel N → e + e − ν. Limits have been set by searches at the Rovno [120] and Bugey [121] reactors. This effect was also searched for by the Borexino experiment [122], which detected neutrinos produced by the fission processes in the Sun -heavy neutrinos with masses up to 14 MeV can be produced in the decay of 8 B. Borexino has set the best limits; |V eN | 2 10 −6 − 10 −5 for m N ∼ 10 MeV.

Active-sterile neutrino oscillations
Persistent anomalies in neutrino oscillation experiments are still providing intriguing hints for the existence of an additional mass squared splitting ∆m 2 ∼ 1 eV 2 to the wellestablished solar and atmospheric mass squared splittings ∆m 2 sol = 7.55 × 10 −5 eV 2 and |∆m 2 atm | = 2.5 × 10 −3 eV 2 , respectively [123,124]. This apparent splitting has been established in the measurement of multiple oscillation processes, including ν µ → ν e accelerator neutrino appearance (LSND anomaly),ν e →ν e reactor neutrino disappearance (reactor anomaly) and the ν e → ν e disappearance of 37 Ar and 51 Cr electron capture decay neutrinos (gallium anomaly). Attempts have been made to fit the data to models with additional eV-scale neutrinos, e.g. (3+1) and (3+2) phenomenological models. While recent reactor experiments such as DANSS [125] and NEOS [126] have improved the statistical significance of an additional eV-scale sterile state, when combined with the ν e appearance data of MiniBooNE they are in strong tension with the observed ν µ → ν µ accelerator neutrino disappearance of the MINOS, NOνA and IceCube experiments.
In the context of the single-generation simplification of this work we interpret the mass squared splitting to be ∆m 2 41 = m 2 N −m 2 ν . As we are focused on the electron-sterile coupling it is thus only the ν e → ν e andν e →ν e experiments sensitive to sin 2 2θ ee ≈ 4|V eN | 2 that are relevant. For sub-eV sterile neutrino masses the Daya Bay [127], KamLAND [54] and upcoming JUNO [128] experiments can probe the mixing down to |V eN | 2 10 −3 . However it should be noted that if one wants to fit the solar and atmospheric mass splittings in a minimal (3+1) or (3+2) extension, solar data excludes the region 10 −9 eV < m N < 0.6 eV [47,55]. Below this region is the pseudo-Dirac scenario and above the mini-seesaw extending to the conventional high-scale seesaw. Light sterile neutrinos can be implemented in the context of an inverse seesaw as considered in refs. [34,129,130].
In figures 6 and 7 we therefore start m N at the eV-scale. Above this the DANSS and NEOS experiments provide limits down to |V eN | 2 10 −2 (as both exclusions are similar, figure 6 shows NEOS only) while the operating PROSPECT [131] experiment provides constraints up to m N = ∆m 2 41 + m 2 ν ∼ 5 eV. Over the same mass range Super-Kamiokande, IceCube and DeepCore (SK+IC+DC) provide complementary limits [132]. We note that the above limits are from oscillations conserving total lepton number. While it is in principle possible to observe LNV in oscillations, this requires new physics beyond sterile neutrinos such as right-handed currents [133].

Electroweak precision data and other indirect laboratory constraints
Any mixing between active and sterile neutrinos necessarily induces non-unitarity effects among the active neutrinos visible in CC and NC processes [134][135][136]. This is most easily parametrised by a non-unitary light neutrino mixing matrix where the matrix η measures deviations from unitarity. The elements of η are given in a generic seesaw model by 2|η | = i V N i V * N i and alter electroweak precision data (EWPD) observables. These include leptonic and hadronic measurements of the weak mixing angle s 2 W , the W boson mass m W , ratios of fermionic Z boson decay rates R l , R c , R b and σ 0 had , the Z invisible decay width Γ inv Z and ratios of leptonic weak decays testing EW universality R π , R W , R K and R l . Furthermore, by modifying G F , the non-unitarity of U ν impacts the values of CKM mixing matrix elements extracted from experiments. Numerous weak decays have been used to pin down the CKM elements V ud , V us , V ub and the unitarity condition |V ud | 2 + |V us | 2 + |V ub | 2 = 1. Assuming a single sterile state coupling to just the first generation, all of these measurements enforce a constant bound of 2|η ee | = |V eN | < 0.050 for m N 1 GeV [136][137][138][139][140][141].
Another indirect measurement of η and hence different combinations of the activesterile mixings comes from the non-observation of lepton flavour violating (LFV) processes → γ and µ − e conversion in nuclei [142]. Due to the different flavours of charged leptons involved in these processes, active-sterile mixings to at least two active generations are required. For the purpose of our single active generation picture we may convert the constraint on |V eN V * µN | obtained from the limits Br(µ → eγ) < 4.2 × 10 −13 [1] and R Ti µ→e < 4.3 × 10 −12 [143] to a constraint in the m N − |V eN | 2 parameter space by assuming |V µN | = |V eN |. We find |V eN | 2 10 −3 for m N ≈ 10 GeV, improving to |V eN | 2 10 −5 for 100 GeV m N 10 TeV. In making the assumption |V µN | = |V eN | however, the constraints in the m N − |V µN | 2 parameter space equally apply for |V eN | 2 . For clarity and consistency we therefore do not show the LFV constraints in figure 6.

Cosmological and astrophysical constraints
The presence of sterile states with masses m N and mixings |V N | 2 (and therefore predicted production rates, decay lengths and active-sterile oscillations) can have drastic consequences on early-universe observables, and have been explored extensively in the literature [144]. These include the abundances of light nuclei formed during Big Bang Nucleosynthesis (BBN), temperature anisotropies in the Cosmic Microwave Background (CMB) radiation and the large-scale clustering of galaxies. Deviations from the standard smooth, isotropic background evolution and perturbations around this background impose severe constraints, especially for sterile states with masses m N 100 MeV. The limits are however highly sensitive to the production and decay mechanism of the sterile state and can be relaxed in certain models. For the purpose of comparison we consider the following scenarios: • Sterile neutrinos with masses m N 1 GeV can be sufficiently long-lived to disrupt the standard formation of light nuclei 4 He, D, 3 He and 7 Li during BBN [145,146].

JHEP03(2020)170
For larger masses the decay products from the accessible two-body and three-body decays have enough time to thermalise with the plasma. For decay times τ 1 s occuring below T 1 MeV, i.e. roughly after ν/N -e ± decoupling and the onset of BBN, both the modified background expansion due to the presence of non-relativistic N and the altered weak processes n + ν ↔ p + e − and p +ν ↔ n + e + involving nonthermal decay product neutrinos lead to modified nuclei abundances. The condition τ = Γ −1 N 1 s naively translates to a lower limit of |V eN | 2 10 −11 (GeV/m N ) 5 for N → 3ν, N → νe + e − and the sub-dominant radiative decay N → νγ. Above the pion mass threshold the already considerably less stringent constraints are made even weaker by including the decays N → νπ 0 and N → e ± π ∓ .
• Sterile neutrinos decaying at later times (with τ t rec ≈ 1.2 × 10 13 s) to nonthermally distributed active neutrinos can modify the amount of dark radiation measured (beyond the usual value including active neutrino oscillations, N eff 3.046) at recombination, ∆N eff . Decays after recombination but before the current epoch (t rec τ t 0 ≈ 4.3 × 10 17 s) can also be important. Useful probes of these effects on the smooth, isotropic expansion history include the CMB shift parameter R CMB (related to the position of the first acoustic peak in the CMB temperature power spectrum), the first peak of Baryon Acoustic Oscillation (BAO) sound waves imprinted on the large-scale distribution of galaxies and finally the value of the Hubble parameter H(z) inferred from Type Ia supernova, BAO and Lyman-α survey data. These exclude values of m N and |V eN | 2 corresponding to lifetimes up to t 0 , where the condition that N does not make up more than the observed matter density Ω sterile < Ω DM ≈ 0.12 h −2 and thus overcloses the Universe also applies. This constraint can naturally be evaded in exotic models [147][148][149][150], for example those that inject additional entropy and dilute the dark matter (DM) energy density. We indicate the combined constraints from ref. [151] in figure 6 as CMB+BAO+H 0 .
• Sterile neutrinos with masses 1 keV m N 100 keV can avoid the global constraints above if the active-sterile mixing is sufficiently small, i.e |V eN | 2 10 −10 − 10 −8 . With lifetimes longer than the current age of the Universe these sterile states are viable DM candidates if efficiently produced [119,152,153]. Depending on the size of the lepton-antilepton asymmetry η L ≡ n L /n γ , population can occur either through resonant (η L > 10 6 η b ) or non-resonant (η L ≈ 0) active-sterile oscillations. The former (Shi-Fuller mechanism [154]) is independent of |V N | 2 while the latter (Dodelson-Widrow mechanism [155]) requires values of |V N | 2 now excluded by the global constraints. If DM is composed entirely of keV sterile neutrinos their fermionic nature limits the phase space density of DM-rich dwarf galaxies and imposes the Tremaine-Gunn bound, m N 0.4 keV. It is also possible to search for anomalous X-ray lines from the radiative decays N → νγ in the diffuse X-ray background and from DMrich astrophysical objects. An intriguing signal at E 3.55 keV implying a sterile neutrino with a mass of 7.1 keV has continued to persist in observations of stacked galaxy clusters [156], the Perseus galaxy cluster and Andromeda M31 galaxy [157] and the centre bulge of the Milky Way [158]. In figure 6 we include the most recent JHEP03(2020)170 observations of M31 and the Milky Way by NuSTAR [159,160]. In figure 7 we show the slightly improved future sensitivity of ATHENA [161]. These constraints assume Ω DM = Ω sterile , but can be multiplied by Ω sterile /Ω DM to account for other DM species [151].
• Active-sterile mixings can be excluded for sterile neutrinos in the mass range 10 eV m N 10 keV by examining their impact on Type II Supernovae. Active-sterile neutrino oscillations hinder the standard neutrino reheating of the reflected shock wave which becomes stalled in the first fraction of a second after the core bounce. For the explosion to proceed and additionally produce the observed SN1987Aν e signal of terrestrial detectors such as Kamioka [162] and IMB [163], a certain region of the m N − |V N | 2 parameter space must be excluded. Refs. [164][165][166][167][168][169][170] have studied in detail the resonant conversion ν e → N in the dense medium of collapsing stars and the necessary conditions to prevent impeding the supernova explosion. Refs. [171][172][173] have similarly investigated ν µ,τ → N conversions for which the Mikheyev-Smirnov-Wolfenstein resonance conditions are different. An open question is whether the conditions for r-process nucleosynthesis to produce heavy elements in the supernova outflows are met in these cases [166,169]. Lastly, sterile neutrinos that escape supernovae can subsequently decay radiatively via N → ν e γ and N → ν e e + e − γ, producing an excess of gamma rays arriving soon after the detection of the ν e . The non-observation of such an excess for SN1987A provides a stringent limit in the mass range 1 MeV m N 30 MeV [174]. Given the various assumptions and calculational differences of the constraints discussed we show for illustration in figure 6 the excluded region from ref. [165].
• Sufficiently stable and light sterile neutrinos with masses m N 50 eV can be produced with quasi-thermal temperatures before the decoupling of active neutrinos via active-sterile oscillations [119,175,176]. While relativistic they contribute themselves towards the extra effective number of light fermionic degrees of freedom ∆N eff . Once becoming non-relativistic they contribute towards the matter density as Ω sterile h 2 = (m sterile eff /94.1 eV) while also damping density perturbations below a mass-dependent free-streaming scale. The most simple case of a single sterile neutrino thermalising through oscillations at the active neutrino temperature has ∆N eff = 1 and m sterile eff m N [54,177,178] which is now likely excluded [179]. The Planck collaboration has made fits of CMB (TT+lowP+lensing+BAO) data to the parameters ( m ν , N eff ) and (m sterile eff , ∆N eff ) [8]. In refs. [151] and [180] these constraints are mapped to the (∆m 2 41 , sin 2 2θ ee ) parameter space which we use to plot the grey dot-dashed CMB constraints in figure 6.
Of particular importance is the dependence of the 0νββ decay rate on the sterile neutrino mass m N and the average momentum exchange squared of the process p 2 . We will see that if m 2 N > p 2 the contribution from a 'heavy' sterile neutrino is suppressed by 1/m N and |V eN | 2 . In the limit m 2 N p 2 the heavy states are integrated out and 0νββ decay becomes a probe of generic short and long-range exchange mechanisms with dimension-7 and above effective operators (depending on the model of interest) at the interaction vertices [193]. If m 2 N p 2 the 'light' sterile neutrino contributes much like a light active neutrino. In this case the condition (M ν ) 11 = i V 2 ei m i = 0 suppresses the total 0νββ decay rate [T 0ν 11 . Multiple sterile states, some with masses above and some below p 2 is an intriguing intermediate scenario. It was observed in ref. [181] that here the 'light' sterile neutrino contribution may even dominate over the light active contribution; the necessary and contradictory prerequisites are a large source of LNV and a small loop contribution to the light neutrino masses. This was found to be possible either in an extended seesaw or by having fine-tuned cancellations between generations.
We will also give a broad comparison between the discussed 0νββ decay constraints and those from the numerous searches discussed in section 4, particularly where the 0νββ decay constraints become relevant (m N 100 keV). One of the most interesting aspects of this comparison is the change of the 0νββ decay constraints as a function of the mass splitting between the heavy states ∆m N . Because 0νββ decay is an LNV process we know specifically in the inverse seesaw that it must vanish in the LNC limit µ R,S → 0. The LNV matrices µ R,S also control the splitting between the heavy states, so in the limit ∆m N → 0 (the heavy states form a pseudo-Dirac fermion) the 0νββ decay limits vanish. Following section 4.2, we will compare this with the suppression of LNV collider and meson decay constraints. No such suppression occurs for the LNC search constraints discussed generally in section 4.
It is also crucial to consider how the sterile neutrino mass splitting ∆m N affects the interpretation of the direct searches. For example, the analyses of β-decay kink searches and meson decay peak searches assume a single sterile state and constrain the associated mixing |V eN | 2 and mass m N . On the other hand, it could be the case that there are two sterile neutrinos with a splitting ∆m N below the energy resolution of the experimentthe searches are then sensitive to the sum of mixings |V N 1 | 2 +|V N 2 | 2 . It is easy to see that, again in the single-generation case, there is a lower limit on this sum from the (M ν ) 11 = 0 condition (or the requirement to produce the observed light neutrino mass m ν ), where we assume r ∆ 1. This is qualitatively identical to the discussion of ref. [194], where it is made clear that for any individual mixing V N i it is not possible to impose a lower limit from the seesaw relation because we are free to set |V N 1 | 2 = 0 and |V N 2 The equivalent freedom in the three-generation picture can be for example the choice of orthogonal matrix R entering the generalised Casas-Ibarra parametrisation [46,47]. If JHEP03(2020)170 ∆m N is instead larger than the energy resolution of direct searches, the non-observation of a sterile state excludes regions in both the m N 1 − |V N 1 | 2 and m N 2 − |V N 2 | 2 parameter spaces. As direct searches have so far only probed mixing strengths viable in the inverse seesaw region of the parameter space, |V N 1 | 2 ≈ |V N 2 | 2 (1 + r ∆ ), the excluded region in m N 2 − |V N 2 | 2 excludes additional portions of m N 1 − |V N 1 | 2 . In our subsequent figure 12 this is simply represented in the excluded region shifted to smaller m N 1 and larger |V N 1 | 2 by the factor (1 + r ∆ ).

Coherent contribution of light and heavy neutrinos
The 0νββ decay rate or inverse half-life, taking into account the exchange of both three active and n S sterile neutrinos, can be written as where G 0ν is a kinematic phase space factor for the outgoing electron pair, g A the axial coupling strength, m p the proton mass and M 0ν (m i ) the nuclear matrix element (NME) of the process for an exchanged Majorana neutrino of mass m i [195].
The most recent calculations of G 0ν for relevant 0νββ decay isotopes have included effects such as the Coulomb distortion of the electron wave functions due to the finite size of the daughter nucleus and electron screening [196][197][198]. The NMEs are in principle far more difficult to compute as they encode the non-trivial transition between the initial and final state nuclei in the process. The NMEs entering eq. (5.2) take the form where J µ is the hadronic current, R the nuclear radius and ω i = p 2 + m 2 i the energy of the exchanged neutrino. It is necessary to sum over all possible intermediate nuclear states n between the initial and final states I and F respectively, and µ = E n − 1 2 (E I + E F ) is the relative energy of these virtual states with respect to the average energy of the process. This sum, along with the non-perturbative nature of the hadronic currents, has made the calculation of eq. (5.3) extremely difficult, and at present there are still large theoretical uncertainties in computed values. Four common simplifying assumptions are (i) the closure approximation, (ii) the impulse approximation, (iii) J P = 0 + final nuclear states and (iv) electrons emitted in s-wave. (i) assumes that only exchanged neutrino momenta |p| of similar size to the nucleon-nucleon spacing contribute to the amplitude -this allows the denominator in eq. (5.3) to be pulled out of the sum and removes the contribution of intermediate odd-odd nuclei. (ii) allows the expression of the hadronic current matrix elements in terms of the nucleon-level current form factors associated with the vector (g V ), axial-vector (g A ), induced weak-magnetic (g M ) and induced pseudo-scalar (g P ) couplings. As 0νββ decay parent and daughter isotopes have even numbers of protons and neutrons, their ground state is always J P = 0 + , while decays to excited states are suppressed, thus justifying the assumption (iii). Finally, p-wave emitted electrons are also suppressed and the computation of G 0ν is greatly simplified in the s-wave case, as assumed in (iv).
A useful interpolating formula for the NMEs can be derived examining the limits of eq. (5.3) for the neutrino mass much smaller and much larger than the average momentum exchange, where M 0ν ν and M 0ν N are dimensionless 'light' and 'heavy' NMEs respectively. It is possible to write an approximate interpolating formula that includes both of these scaling behaviours, so that the half-life formula (5.2) including sterile states becomes [103,199] 1 T 0ν (5.6) Using the above-discussed approximations in eq. (5.3) the values of |M 0ν ν | and |M 0ν N | have been calculated in a variety of different frameworks. These include the quasiparticle random phase approximation (QRPA) [199,200], interacting boson model (IBM-2) [201,203,204] and interacting shell model (ISM) [202]. A review of these methods as well as their respective strengths and weaknesses is given in ref. [205]. In table 1 we show the light and heavy NMEs and their associated fractional uncertainties for the 0νββ decay isotopes 76 Ge and 136 Xe. The QRPA calculations of the Tübingen and Jyväskylä groups and the IBM-2 calculations of the Yale group give NME values for quenched (g A = 1) and non-quenched (g A = 1.269) values of the axial coupling and also for phenomenological Argonne [206] and CD-Bonn [207] forms of the Jastrow potential describing two-nucleon short-range correlations. We use the average of these NME values and take the uncertainty JHEP03(2020)170 to be half the maximum spread. It was noted in ref. [200] that the QRPA Jyväskylä and IBM-2 Yale heavy NMEs change by a common factor when changing potentials, while for an unknown reason the changes for the QRPA Tübingen heavy NMEs are significantly different. There are now numerous other computational tools being used for ab initio calculations of light NMEs for both light and heavy nuclei, including improved chiral effective field theory [208], renormalisation group [209,210] and lattice QCD techniques [211].
In figure 8 we plot the 76 Ge and 136 Xe NMEs as a function of the exchanged neutrino mass m N i using the interpolating formula of eq. (5.5) and the different light and heavy NMEs given in table 1. It can clearly be seen that the NMEs are constant below p 2 ∼ 100 MeV 2 and suppressed by 1/m 2 N i above. If all masses are below p 2 we will see that it is instead the seesaw relation suppressing the 0νββ decay rate. To plot the uncertainty bands in figure 8 we propagate the uncertainties of |M 0ν ν | and |M 0ν N | through eq. (5.3) as It can be seen that the largest uncertainties are in the IBM-2 NMEs -for illustrative purposes and to give conservative estimates we use these NMEs in the following discussion. In our single-generation simplification the summation appearing in the interpolating formula is approximately where we have used the approximate seesaw relation m ν + e iφ 1 m N 1 s 2 e1 = −e iφ 2 m N 1 (1 +

JHEP03(2020)170
r ∆ )s 2 e2 to eliminate s 2 e2 and rewrite the summation using the factors Alternatively, one could eliminate s 2 e2 and φ 2 using the exact seesaw relations eqs. (3.9) and (3.10). However, taking the small mixing approximation c 2 e1 ≈ c 2 e2 ≈ 1 as done above makes a very small difference to the following results. It is easy to see that these two substitutions are equivalent -if we set s 2 e1 = 0 in eq. (5.8) we would be left with the contributions from the light and second heavy state. There is a relative minus sign between terms because in this limit φ 2 = π in both the canonical seesaw φ 1 = π and inverse seesaw φ 1 = 0 cases (and any intermediate φ 1 value), as can be seen in figure 3. Taking the square of the summation in eq. (5.8) and inserting into eq. (5.6) now gives Experimental lower bounds on the 0νββ decay half-life T 0ν 1/2 > (T 0ν 1/2 ) exp (or χ 2 < χ 2 exp ) can therefore be used to put an upper bound on s 2 e1 as a function of and (through dependence on p 2 and χ exp ) the light and heavy NMEs |M 0ν ν | and |M 0ν N |, Of course there is another limit derived from the quadratic inequality χ 2 < χ 2 exp , that is a lower bound on s 2 We will see that for most choices of parameters this is negative and unphysical. It will be important when cos φ 1 < 0 and α > χ exp . As detailed earlier, we work in a one-generation framework with one light state ν L,e which we identify with the electron neutrino. As such, the effective 0νββ mass is not a coherent sum as usually defined, but is simply given by the electron neutrino mass, m ββ = m ν . (5.14) In our parametrization, m ν is always real and positive. We calculate 0νββ consistently in this framework; specifically, we include a coherent summation in the contributions of the light neutrino state and the two heavy neutrino states including the relative phases as JHEP03(2020)170  detailed above. In this sense, m ν is a surrogate for the general effective 0νββ mass m ββ in eq. (5.13), but we cannot include the potential destructive interference therein due to the Majorana phases within the light PMNS matrix. This effect has been discussed extensively in the literature (for a review, see e.g. [212]) whereas our focus is on the constraints on the heavy neutrino parameters. The precise value of the light neutrino contribution may only be important if it saturates the limit from 0νββ decay searches. This is discussed in figure 11 (right) and the accompanying text where the choice m ν = 6 × 10 −2 eV is near the excluded m ββ limit and thus the constraints on the extra contributions of the heavy neutrinos become overly restrictive. These may instead be relaxed if there is a sizeable cancellation among the light neutrino contributions reducing m ββ . A full analytic discussion of 0νββ in presence of three active neutrinos mixing with sterile neutrinos is beyond the scope of the present work and will be pursued in a follow-up work.

Sensitivity to sterile neutrino parameters
In figure 9 we display the upper bounds on the sum of squared active-sterile mixings |V eN 1 | 2 + |V eN 2 | 2 ≈ s 2 e1 + s 2 e2 as a function of the first sterile neutrino mass m N 1 for three

JHEP03(2020)170
small values of the sterile neutrino mass splitting ratio r ∆ 1 and for benchmark values of the light neutrino mass m ν = 10 −3 eV and Majorana phase φ 1 = 0. The sum is used in the assumption that for small splitting the energy resolutions of direct searches are larger than ∆m N and consequently constrain |V eN 1 | 2 + |V eN 2 | 2 as a function of the mass m N 1 ≈ m N 2 . Making use of the s 2 e1 inequality in eq. (5.11) we take the most recent lower limits on T 0ν 1/2 from the 136 Xe KamLAND-Zen [213] and 76 Ge GERDA-II [214] experiments and the IBM-2 light and heavy NMEs in table 1 to plot the solid (and dashed) curves in the upper right portion of figure 9. The bands illustrate the uncertainty on |V eN 1 | 2 + |V eN 2 | 2 as a function of m N 1 found by propagating the conservative IBM-2 uncertainties through eq. (5.11). The red curves in figure 9 depict the upper limits on |V eN 1 | 2 + |V eN 2 | 2 when including only the contribution of a single sterile state (neglecting light active exchange) towards 0νββ decay. Finally, we show for these choices of r ∆ the upper limits on |V eN 1 | 2 + |V eN 2 | 2 from the requirement that |δm 1−loop ν | < 0.1m ν , taken directly from figure 4 (right). We compare these 0νββ decay bounds to the direct search limits discussed in section 4. These include the current (blue-shaded) and future (blue dot-dashed line) sensitivities of LNC probes including β-decay kink searches, meson decay peak searches, beam dump experiments and collider constraints. We also display separately the current (red-shaded) and future (red dot-dashed line) sensitivities of LNV meson decay and collider probes. Faint grey regions correspond to the cosmological excluded regions. Finally, the dark grey shaded region below the seesaw line |V eN 1 | 2 + |V eN 2 | 2 = mν m N 1 is excluded as explored at the start of section 5.
In figure 10 we similarly show the upper bounds from 0νββ decay and loop considerations for the same (small) values of the sterile mass splitting ratio but instead use the predicted sensitivity of future experiments, T 0ν 1/2 10 28 y. This reach may be achievable at the proposed PandaX-III [215]  We first observe that the upper bounds are most stringent for m N 1 ∼ p 2 ∼ 200 MeV, reaching |V e | 2 10 −7 . Towards lower m N 1 both sterile states are 'light' and the 0νββ decay rate is suppressed by the seesaw relation, eventually erasing the upper bounds for m N 1 1 MeV. For higher m N 1 both sterile states are 'heavy' and the limits become weaker as m N 1 increases due to the growing NME suppression by 1/m 2 N 1 . We also see a strong dependence on the sterile mass splitting ratio; decreasing r ∆ by a factor of ∼ 10 2 weakens the upper bound by a similar factor both above and below m N 1 ∼ p 2 . This is to be expected as r ∆ → 0 corresponds the pseudo-Dirac limit in which lepton number is approximately conserved and the 0νββ decay process is forbidden. Comparing the 76 Ge and 136 Xe bounds it is interesting to note that those for the former are slightly more stringent despite the smaller experimental half-life lower bound. As can be seen in figure 8 this is counteracted by 76 Ge possessing larger NMEs on average compared to 136 Xe. Comparing with direct searches we see that for these small choices of r ∆ the current upper bounds are at best comparable with non-resonant meson decay limits for 1 MeV < m N 1 < 1 GeV and more stringent than collider constraints for m N 1 > 5 GeV.
We saw in section 4.2 that when the sterile mass splitting ratio r ∆ is decreased the LNV collider constraints shaded in red do not weaken significantly -this is because the JHEP03(2020)170  amplitide of LNV is controlled by the ratio Γ N /∆m N . By considering the open sterile neutrino decays to SM particles we found Γ N /∆m N = Γ N /(r ∆ m N 1 ) 1 in the mass range 5 GeV m N 1 50 GeV for r ∆ 10 −10 . Thus, when r ∆ 10 −2 the 0νββ decay constraints become less stringent than the same-sign dilepton and LNV trilepton collider constraints.
The behaviour of the 0νββ decay upper bound in the 'light' and 'heavy' regimes can be quantified by taking the Taylor expansion of eq. (5.11) in the opposing limits m N 1 / p 2 1 and m N 1 / p 2 1. In the light regime we derive , (5.15) while in the heavy regime The m N 1 dependence of these upper bounds agrees qualitatively with figure 9 -in the light regime the upper bounds scale as 1/m 3 N 1 and in the heavy regime as m N 1 . The dependence on r ∆ is also in agreement -for r ∆ 1 both eqs. (5.15) and (5.16) are inversely proportional to r ∆ . Thus decreasing or increasing r ∆ shifts the entire upper bound to higher and lower mixings for the whole range of m N 1 . In figure 11 we study more closely the |V eN 1 | 2 + |V eN 2 | 2 = s 2 e1 c 2 e2 + s 2 e2 upper bound in the r ∆ = 10 −2 case for different values of the Majorana phase φ 1 (left) and the light neutrino mass m ν (right). To the left it is clear that changing φ 1 has little effect on the 0νββ decay constraints for these choices of parameters. As shown in eq. (5.15), in the light regime the s 2 e1 upper bound is independent of φ 1 because the suppression of the 0νββ decay rate through the seesaw relation is also independent of φ 1 . From eq. (5.16) we see that in the heavy regime changing φ 1 has little effect for these parameter choices because m ν / p 2 χ exp , i.e. the light neutrino contribution is negligible. 0νββ decay is therefore driven by the two heavy states. This is the limit α 1 and s 2 e1 χexp β in eq. (5.11). For m N 1 p 2 1/2 we have
To the right we see that the effect of increasing m ν for φ 1 = 0 is to strengthen the upper bound in the heavy regime. This again is described by eq. (5.16) -there is a cancellation between the two terms in the brackets as m ν / p 2 approaches χ exp . In this limit the light active contribution becomes non-negligible compared to the difference between the heavy sterile contributions. For the inverse seesaw region of the parameter space If for example (φ 1 , φ 2 ) = (0, π), the light contribution adds constructively with the difference and the upper bound on s 2 e1 (multiplying the heavy contributions) must be smaller to JHEP03(2020)170 account for the observed half-life lower bound. If on the other hand (φ 1 , φ 2 ) = (π, 0), the light and heavy contributions add destructively and the s 2 e1 upper bound can be relaxed. If m ν / p 2 > χ exp (which may be the case for a large lower limit on T 0ν 1/2 ) no value of s 2 e1 in the heavy regime is permitted for φ 1 = 0. In eq. (5.11) this corresponds more generally to the case α > χ exp in which the upper bound on s 2 e1 becomes negative and unphysical. Constructive interference between the light active contribution and the difference between the heavy sterile contributions, e.g. as for (φ 1 , φ 2 ) = (0, π), now gives a T 0ν 1/2 less than the experimental lower limit, or χ > χ exp . Conversely, if the light and heavy contributions interfere destructively, e.g. for (φ 1 , φ 2 ) = (π, 0) above the seesaw line and (φ 1 , φ 2 ) = (π, π) below, then s 2 e1 multiplying the heavy contributions can be made large enough to meet the condition χ < χ exp (but not so large as to dominate over the light contribution). As well as an upper bound, this sets a lower bound on s 2 e1 in the heavy regime. This is the lower bound in eq. (5.12) becoming non-negative.
It is worth reminding the reader that we are considering a value of r ∆ in the range [0, ∞] and so the introduced quantities in eq. (5.9) satisfy α > 0 and β > 0. As explained in section 3, a value of r ∆ in the range [−1, 0] is equivalent to swapping the roles of the sterile states, now having m N 2 < m N 1 . In this equally valid range the introduced quantities satisfy α > 0 and β < 0. Because of this we see by examining eqs. (5.11) and (5.12) that the behaviours of the active-sterile mixings (s 2 e1 , s 2 e2 ) and Majorana phases (φ 1 , φ 2 ) are also swapped, with cancellation between light active and heavy sterile contributions taking place for (φ 1 , φ 2 ) = (0, π).
In figure 11 we also see how the loop constraints change when varying φ 1 and m ν . For the extreme values φ 1 = 0, π and intermediate value φ 1 = π 4 the loop constraints are broadly the same. However for φ 1 = π 2 the upper bound becomes nearly two orders of magnitude more stringent. As m ν is increased by an order of magnitude (we do not go to m ν > p 2 χ exp ≈ 0.083 eV for the reasons discussed previously) we can also see that the loop constraints are correspondingly weakened by an order of magnitude.
In figure 12 we display the active-sterile mixing |V eN 1 | 2 ≈ s 2 e1 as a function of m N 1 for three large values of the sterile neutrino mass splitting ratio r ∆ ≥ 1 and for benchmark values of the light neutrino mass m ν = 10 −3 eV and Majorana phase φ 1 = 0. We do not show the sum in this case because it is assumed that the splittings are large enough for the two states to be resolved individually in direct search experiments. We compare these bounds to the direct search limits discussed in section 4. Due to the large splitting, shifted versions of the excluded region depending on the value of r ∆ now apply -this is a shift to smaller m N 1 and to larger |V eN 1 | 2 by a factor (1 + r ∆ ). For example, if the T2K experiment excludes a second state of mass m N 2 and mixing |V eN 2 | 2 , it also implies the non-existence of the first state at m N 1 ≈ m N 2 /(1 + r ∆ ) and |V eN 1 | 2 ≈ |V eN 2 | 2 (1 + r ∆ ). These particular relations apply because the T2K bounds are in the inverse seesaw region of the parameter space. For these large splittings we immediately see that the 0νββ decay constraints converge towards the upper bound in the limit of single heavy neutrino exchange (commonly used in the literature), shown by the thin red curve in figures 9, 10 and 12.
We have so far neglected the one-loop contribution to the neutrino mass δm which would be expected to alter the suppression and 1/m 3 N 1 scaling due to the tree-level seesaw relation. However, when we look at figure 4 we see that |δm 1−loop ν | ∼ 10 −12 eV in the light regime, even after the iterative procedure on δm 1−loop ν is applied. Thus we safely expect this effect on the 0νββ decay constraint curves to be negligible.

Conclusions
Heavy sterile neutrinos represent one of the most interesting candidates for particles beyond the Standard Model. They are conspicuously absent from the Standard Model particle content which means SM neutrinos are the only fermions that do not have an electroweak singlet partner field. It is not far-fetched to assume that this has something to do with the fact that sterile neutrinos, as their name implies, are singlets under all the SM gauge groups and a Majorana mass term breaking total lepton number is therefore not protected.
There is a strong ongoing and planned effort to search for sterile neutrinos over a wide range of masses and active-sterile mixing strengths. The main focus of this work is to JHEP03(2020)170 compare direct searches such as at the LHC and in meson decays with constraints from 0νββ decay. The latter is the most important probe of lepton number violation and light Majorana neutrino masses. Heavy neutrinos will generically contribute to 0νββ decay as well. They are thus constrained by current searches and can be probed in future 0νββ decay experiments.
In this work, we have introduced a phenomenological parametrisation of a onegeneration seesaw model in terms of experimentally measurable quantities, such as activesterile neutrino mixing angles, CP phases, masses and mass splittings. We have identified the regions of parameter space allowed by consistency conditions in the neutrino mass matrix in the single-generation case, and have showed how the type-I and inverse seesaw limits can be recovered (cf. figure 2). Imposing the additional consideration that the loop contribution to the active neutrino mass must be less than 10% of the tree-level mass further reduces this allowed parameter space, as shown in figure 5.
We summarise current and future experimental constraints on the sterile neutrino mass-mixing parameter space over a wide range of interest, including both lepton number conserving and violating processes (cf. figures 6 and 7), emphasising that the LNV constraints could change depending on the mass splitting between the two sterile states and the relative CP phase between them. This is particularly relevant for 0νββ decay searches, which are significantly weakened for quasi-Dirac sterile neutrinos, as shown in figure 9, while for large mass splitting, the 0νββ decay constraint remains strong in the electron sector; cf. figure 12, and it is especially relevant for heavy neutrino masses in the region m N ≈ 100 MeV to a few GeV, where the future 0νββ decay sensitivities can reach a level close to the small active-sterile mixing strengths expected in a vanilla seesaw scenario.  Figure 13. Constraints on the mass m N of the sterile neutrino and its squared mixing |V µN | 2 with the muon neutrino (above) and |V τ N | 2 with the tau neutrino (below). The shaded regions are excluded by the searches listed in the appendix.

JHEP03(2020)170
NOνA [236], CDHS [237] and CCFR [238]. It should be noted that limits on the disappearance channel, i.e. active to sterile oscillation, tend towards a constant upper bound as a function of ∆m 2 and hence m N . This bound can in principle be extended to arbitrarily large m N , covering the region between 100 eV and 1 MeV. Finally, limits have been set from considerations of supernovae [173] and the non-observation of X-rays [159,160]. Alternatively, upper limits on the tau-sterile mixing strength have been placed above m N ∼ 1 MeV by NOMAD [239], CHARM [240], Super-Kamiokande, T2K, lepton universality and B decays [241], L3, DELPHI and electroweak precision data. From oscillations upper bounds have been placed by IceCube [232,233], Super-Kamiokande [235] and NOνA [236]. The supernovae and X-ray constraints also apply for |V τ N | 2 .
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.