Looking for the minimal inverse seesaw realisation

In this work we consider a simple extension of the Standard Model involving additional fermionic singlets and assume an underlying inverse seesaw mechanism (with one or more right-handed neutrinos and one or more sterile fermions) for neutrino mass generation. Under the assumption that both sterile states and right-handed neutrinos are present, our goal is to determine which is the minimal inverse seesaw realisation that accounts for neutrino data while at the same time complying with all experimental requirements (electroweak precision tests and laboratory constraints). This study aims at identifying the minimal inverse seesaw realisation for the 3-flavour and for the 3 +~more-mixing schemes, the latter giving an explanation for the reactor anomalies and/or providing a possible candidate for the dark matter of the Universe. Based on a perturbative approach, our generic study shows that in the class of inverse seesaw models giving rise to a 3-flavour flavour mixing scheme, only two mass scales are relevant (the light neutrino mass scale, $m_\nu$ and the mass of the right-handed neutrinos, $M_R$) while in the case of a 3 + 1-mixing scheme, an additional mass scale ($\mu$ $\in[m_\nu,M_R]$) is required. For each of the two obtained inverse seesaw frameworks, we conduct a thorough numerical analysis, providing predictions for the hierarchy of the light neutrino spectrum and for the effective mass in neutrinoless double beta decay.


Introduction
Oscillation experiments have established a clear evidence for two oscillation frequencies (∆m 2 ij )implying that at least two neutrino states are massive -as well as the basic structure of a 3-flavour leptonic mixing matrix (for a recent review, see [1]). However, current reactor [2], accelerator [3,4] and Gallium anomalies [5] suggest that there might be some extra fermionic gauge singlets with mass(es) in the eV range. This would imply that instead of the three-neutrino mixing scheme, one would have a 3 + 1-neutrino (or 3 + more) mixing schemes (for a global overview, see [6]).
Sterile fermionic states (not necessarily light) are present in several neutrino mass models and their masses can range from well below the electroweak scale up to the Planck scale. Other than the reactor and accelerator anomalies, the existence of sterile states is also motivated by certain indications from large scale structure formation [7,8]. Nevertheless, due to the mixings of the sterile fermionic states with the active left-handed neutrinos, models with sterile fermions are severely constrained from electroweak (EW) precision observables, laboratory data and cosmology.
In contrast with the huge experimental achievements in determining neutrino oscillation parameters, many questions remain to be answered concerning neutrino properties, as for instance the neutrino nature (Majorana or Dirac), the absolute neutrino mass scale and the hierarchy of the neutrino mass spectrum, which are not yet determined. Finally, and most importantly, one must unveil the neutrino mass generation mechanism at work and which new physics scales are required.
One of the most economical possibilities to account for massive neutrinos is to embed a standard seesaw mechanism [9,10,11] into the framework of the Standard Model (SM). The caveat of these scenarios is that, in order to have natural neutrino Yukawa couplings the typical scale of the extra particles (such as right handed neutrinos, scalar or fermionic isospin triplets) is in general very high, potentially very close to the gauge coupling unification (GUT) scale, thus implying that direct experimental tests of the seesaw hypothesis might be impossible. In contrast, low-scale seesaw mechanisms [12,13,14,15,16], in which sterile fermions are added to the SM particle content with masses around the electroweak scale or even lower, are very attractive from a phenomenological point of view since the new states can be produced in collider and/or low-energy experiments, and their contributions to physical processes can be sizeable.
In view of the strong potential of low-scale seesaw mechanisms, in this work we consider the inverse seesaw (ISS) mechanism [12,15,16] which requires the addition of both right-handed (RH) neutrinos and extra sterile fermions to the SM field content. The distinctive feature of the ISS is that an additional dimensionfull parameter (µ) allows to accommodate the smallness of the active neutrino masses m ν for a low seesaw scale, and still with natural Yukawa couplings (Y ν ∼ O(1)). In turn, this allows for sizeable mixings between the active and the additional sterile states. Such features are in clear contrast with, for instance, the canonical type-I seesaw [9], where O(1) Yukawa couplings require the mass of the right-handed neutrinos to be close to the GUT scale, M R ∼ 10 15 GeV, thus leading to extremely small active-sterile mixings.
Since both RH neutrinos and sterile states are gauge singlets, there is no requirement on their (generation) number from anomaly cancellation. Moreover, in view of the presence of two independent mass scales (the mass of the RH neutrinos and the Majorana mass of the sterile states), associated to gauge singlet fermions, it is only natural to investigate which is the minimal content of the ISS extension of the SM successfully accounting for neutrino data, while at the same time complying with all available experimental and observational constraints.
We thus embed the inverse seesaw mechanism into the SM, considering models with an arbitrary (and different) number of RH neutrinos and of additional sterile states, in order to establish which class of models provides a minimal 3-flavour and 3 + more-mixing schemes. The latter class of realisations (configurations) may offer an explanation to the reactor anomalies or, depending on the mass scales, a (partial) solution for the Dark Matter (DM) problem, in the form of a Warm DM (WDM) candidate [17]. In a first stage, we do not impose a particular mass scale for the (RH) Majorana states nor the hierarchy of the associated light spectrum; likewise, we do not specify a mass range for the sterile fields.
Our study has allowed to identify two classes of minimal ISS realisations that can successfully account for neutrino data: the first leads to a 3-flavour mixing scheme, and requires only two scales (that of light neutrino masses, m ν , and the mass of the RH neutrinos, M R ); the second corresponds to a 3 + 1-mixing scheme, and calls for an additional scale (µ ∈ [m ν , M R ]). For each of these minimal classes, we carried a numerical analysis taking into account all possible bounds associated to the presence of sterile fermions (which constrain the mixings between active and sterile neutrinos for different mass regimes). We also provide predictions regarding the hierarchy of the light neutrino spectrum (normal or inverted) and the effective mass in neutrinoless double beta decay, for each of the minimal realisations identified.
The paper is organised as follows: in Section 2, we briefly review the inverse seesaw mechanism and define the framework; we also determine the generic class of frameworks leading to 3-and to 3 + more-mixing schemes as well as their generic features concerning the different mass scales. In Section 3, we consider all the different constraints from neutrino data, electroweak observables and laboratory measurements applied in the analysis. Section 4 is devoted to the phenomenological analysis of the minimal ISS framework leading to the 3-flavour and to the 3 + 1-mixing schemes, respectively. Our final remarks are given in Section 5. For completeness, some technical details concerning the computation are included in the Appendices.

Towards the minimal inverse seesaw realisation
In this work we consider the inverse Seesaw mechanism [12,15,16] for the generation of neutrino masses and lepton mixings, with a minimal field content. We work in the framework of the SM extended by one or more generations of right-handed neutrinos ν R and additional fermionic singlets s.

The one generation case
We first consider the illustrative one generation case. In the basis n L ≡ (ν L , ν c R , s) T , the neutrino mass term reads: where C ≡ iγ 2 γ 0 is the charge conjugation matrix and the matrix M is given by We assume that there is no term mixing the left-handed neutrino with the fermionic singlet s (∼ ν c L s). In the above, d corresponds to the Dirac mass term. The matrix M also includes a Majorana mass term for the RH neutrino, The values of m and µ in Eq. (3) are arbitrary. However, accommodating neutrino masses of O(eV) implies that both must be very small in the case of the inverse seesaw framework. Assigning a leptonic charge to both ν R and s, with lepton number L = +1 [12,15,16] (such that the Dirac mass term −d * ν L ν R +h.c. preserves the leptonic number), the terms ν T R Cν R and s T Cs violate total leptonic number L by two units. Small values of m and µ are natural in the sense of 't Hooft [18] since in the limit where m, µ → 0, the total lepton number symmetry is restored. In the following, we assume for simplicity that µ and m are of the same order of magnitude.
In order to obtain the tree-level neutrino mass spectrum and the leptonic mixing, we diagonalize the matrix M as [19] U where U is a unitary matrix, and m 0,1,2 correspond to the physical neutrino masses. The mixing matrix is obtained from so that the matrix U diagonalizing M † M is the same as the one in Eq. (4). We determine the neutrino spectrum perturbatively: the perturbations correspond to taking into account the tiny effects of the lepton number violating diagonal entries, The lightest neutrino mass arises from perturbative corrections 1 to the zeroth order m 0 = 0 eigenvalue; the two other states are pseudo-Dirac heavy neutrinos, massive and degenerate. Concerning m 0 , the second order corrections m 2 0 (2) (the first order one gives vanishing contributions) are given by which reduces to the usual inverse seesaw expression once one assumes |d| |n|. The first order corrections to m 2 1,2 (0) = |d| 2 + |n| 2 lift the degeneracy: The corresponding eigenvectors allowing to build the leptonic mixing matrix can be found in Appendix A. Notice that in this approach, the only assumption on the magnitude of the physical parameters, i.e. |m|, |µ| |d|, |n|, is driven (and justified) by the naturalness criterium.

Minimal Inverse Seesaw realisations
In this section, we build the minimal ISS framework complying with experimental observations. The latter lead to the following requirements: • there are 3 generations of neutrino fields with SU (2) L ⊗U (1) Y gauge interactions (#ν L = 3); • there are at least 3 non-degenerate light mass eigenstates. 1 We denote by (n) superscript perturbative corrections of order n.
We extend the one generation matrix of Eq. (2) to the case of several generations of ν R and s fields, so that M now reads d, m, n, µ now being complex matrices. Since M is symmetric (due to the Majorana character of the fields), it follows that m and µ are also symmetric matrices.
In the following, we denote by #ν L , #ν R and #s the number of generations of left-handed, right-handed and sterile fields, respectively. The Dirac mass matrix d arises from the Yukawa couplings to the Higgs boson ( Φ = iσ 2 Φ), where Y is a complex matrix, l α L denotes the left-handed (LH) leptonic doublet, α and β being generation indices. After electroweak symmetry breaking (EWSB), the matrix d is given by and its dimension is The matrix n describes the lepton number conserving interactions involving ν c R and s fields, and its dimension is dim n = (#ν R × #s) .
Finally, the dimension of the (symmetric) Majorana mass matrices m and µ are given by Being gauge singlets, and since there is no direct evidence for their existence, the number of additional fermionic singlets #ν R and #s is unknown. In the following we determine their minimal values when accommodating either a 3-flavour or a 3 + 1 (or more) -flavour mixing schemes. The different possibilities are summarised in Table 1.
The first three columns of Table 1 indicate the total number of additional fermionic singlets #ν R + #s, #ν R and #s, respectively. The fourth column contains the number of massless eigenstates at zeroth order (in the absence of accidental cancellations between the a priori independent entries of the mass matrix). Always in the absence of accidental cancellations, the fifth column displays how many massless eigenstates acquire mass once higher order corrections from perturbations are taken into account (see Appendix A): although massive, these states remain light since the corresponding masses are proportional to entries of m and µ (this can be inferred from the one generation illustrative case, see Eq. (7)). It is important to notice that states which are already massive at zeroth order have masses proportional to the d and n matrix entries. Finally, the sixth column contains information on the number of non-degenerate light mass eigenstates predicted by each of the different ISS configurations considered.
The last two columns provide information on the phenomenological viability of the different ISS realisations. Firstly, neutrino oscillation experiments require at least two independent oscillation   Interestingly, models with 4 different light masses could potentially explain the (anti)neutrino anomalies reported by the short baseline experiments LSND [3] and MiniBooNE [4], the Gallium anomaly in radioactive source experiments [5] and the reactor antineutrino anomalies [2]. Such configurations, leading to a 3 + 1-mixing scheme (see for example [20]) are indicated by a ! (a) in the seventh column of Table 1.

!(s) !
For all cases with a viable mass spectrum -either (s) or (a) -we have then verified if the observed mixing pattern could be successfully reproduced. Should this be the case, a !is present in the eighth column of the table.
As can be seen from the information summarised on Table 1, the simplest model which could accommodate the observed neutrino spectrum is the one with (#ν R = 1, #s = 2), which will be here denoted as "(1,2) ISS". It predicts 4 light eigenstates, two of which are massive; provided that the latter are non-degenerate, one could have two independent mass squared differences (corresponding to the solar and atmospheric mass differences). Notice however that this model cannot provide the observed leptonic mixing matrix U PMNS . This is a consequence of having one of its light mass eigenstates dominated by sterile components, and as such it cannot be identified with a SM active neutrino. A similar problem is present for the "(1,3) ISS" configuration, which although in principle accommodating the correct neutrino mass spectrum fails to provide the observed mixings.
We thus conclude from this simple analysis that the most economical ISS realisation which can successfully accommodate low-energy neutrino data is that with (#ν R = 2, #s = 2), hereafter denoted by "(2,2) ISS" realisation. The latter does not provide an explanation for the reactor anomaly; the next-to-minimal ISS realisation which could explain such anomaly is the one with (#ν R = 2, #s = 3), which we denote by "(2,3) ISS" configuration.
Before addressing in detail the phenomenology of each minimal framework above identified, we will briefly comment on some aspects intrinsic to all ISS realisations.

Different neutrino mass scales
As a function of the number of generations for the sterile fields, the model always exhibits #ν L + (#s−#ν R ) light mass eigenstates. These states would be massless at zeroth order, and their masses arise from higher order corrections (in perturbation) due to the block-diagonal matrix which now generalizes ∆M , see Eq. (6). In addition, the full spectrum contains heavy states with masses In the limit where lepton number is conserved (i.e. ∆M = 0) these states become Dirac particles.
The low-energy phenomenology of these models is determined by two quantities: the scale of the Lepton Number Violating (LNV) parameters µ and the ratio between the scale of the Dirac mass terms d and that of the n mass matrix, denoted by k. To understand the key rôle of these quantities, let us consider again the illustrative one-generation model (i.e. #ν L = #ν R = #s = 1) of Section 2.1. The active neutrino mass of Eq. (7) can be rewritten as m ν = |µ|k 2 /(1 + k 2 ), with k = |d|/|n|. In the realistic case of several generations, d, n, µ are matrices, and these considerations loosely apply to the order of magnitude of their entries. The ratio k is directly related to deviations from unitarity of the leptonic mixing matrix, as shown in Appendix A, Eq. (A.11). Constraints on the non-unitarity of the PMNS matrix impose that k should not be too large; as we will discuss in the section devoted to the numerical analysis, solutions in agreement with experimental data can be found if, and only if, 10 −1 . These features are shared by the different realistic extensions presented in Table 1.
The mass spectrum of the ISS models is thus characterised by either 2 or 3 different mass scales (as illustrated in Fig. 1):

Removing unphysical parameters
The relevant leptonic terms of a general inverse seesaw Lagrangian can be written in the following compact form, where In the above equation α, β = 1, 2, 3, i, j = 1, . . . , #ν R and a, b = 1, . . . , #s. The total number n u of physical and non-physical parameters in the mass matrices present in the Lagrangian of Eq. (16) is equal to and detailed in Table 2.

Matrix
Total number of parameters #s × (#s + 1) Total 18 + #ν R (7 + #ν R + 2 #s) + #s(#s + 1) In order to determine the actual number of physical parameters, one has to identify all independent transformations under which the Lagrangian of Eq. (16) is invariant. One finds four classes of transformations with the following unitary matrices: The number of parameters defining the transformations of Eqs. (19 -22) is n t = 18+(#ν R ) 2 +(#s) 2 , as shown in Table 3, so that the number of physical parameters n p thus reduces to Matrix Number of free parameters Table 3: Number of parameters defining the transformations of Eqs. (19 -22).
Since L kin is invariant under each of the transformations of Eqs. (19 -22), we can use the latter to redefine the fields and cast the mass matrices only in terms of physical parameters. For instance, with the transformations of Eqs. (19,20), one can choose a basis in which the charged leptonic matrix m is real and diagonal, and similarly for the symmetric Majorana mass matrices m and µ (in this case using Eqs. (21,22)). Finally, one can eliminate three phases from the Dirac mass matrix d while keeping m real, via a combination of transformations of Eq. (19) and Eq. (20). In this simple example, there are exactly n p free parameters, as summarised in Table 4.

Effects of fermionic gauge singlets and constraints on the ISS parameters
In addition to succeeding in accommodating neutrino oscillation data, models with sterile fermions are severely constrained, since the mixings of the sterile neutrinos with the active left-handed states might induce contributions to several observables, leading to conflict with experimental data. The mixings of the sterile neutrinos with the active left-handed states imply a departure from unitarity of the 3 × 3 U PMNS matrix, which can have an impact on several observables, inducing deviations from the SM predictions. Bounds on the non-unitarity of the PMNS were derived in [21], using Non-Standard Interactions (NSI). These bounds are especially relevant in our analysis when the masses of the sterile states are heavier than the GeV, but some are still lighter than 174 GeV. If the sterile states are sufficiently light and have large mixings with the active neutrinos (as for example in the inverse seesaw [12], the νSM [13] and the low-scale type-I seesaw [14,22]), then the deviations from unitarity of the PMNS mixing matrix can be sizeable, and lead to (treelevel) corrections to the W ν vertex. This will have a significant impact to several observables, such as corrections to the invisible Z decay width [23], significant contributions to lepton flavour universality (LFU) violation observables [24,25,26], and new contributions to numerous low-energy rare decays.
Another important constraint concerns charged lepton flavour violation (cLFV) since the modified W ν vertex gives rise to cLFV processes, typically at rates higher than the current bounds unless the active-sterile mixings are small [12,15,16,27].
In the case of µ → eγ decays, the rate induced by the presence of the sterile states is given by [28]: where the index i runs over all neutrino states, U is the leptonic mixing matrix obtained after diagonalization of the mass matrix and G is the associated loop function. The current bound on this branching ratio is Br(µ → eγ) < 5.7 × 10 −13 at 90% C.L., as reported very recently by the MEG experiment [29]. This will prove to be the most relevant LFV bound in most of our scenarios with light sterile neutrinos.
Constraints arising from neutrinoless double beta (0ν2β) decay bounds can be particularly relevant, since in the ISS the heavy sterile states also contribute to the process. The effective neutrino mass m νe eff , to which the amplitude of the 0ν2β process is proportional, can receive further corrections with respect to the standard expression, 3 i=1 U 2 e,i m ν i . Since the heavy Majorana states mix to form pairs of pseudo-Dirac states, their contribution is proportional to their mass difference weighted by the ν e -sterile mixing. Each Majorana state thus contributes to the amplitude of a 0νββ decay as [30] A where M 0νββ (m i ) is the nuclear matrix element that characterises the process. The latter is a function of the neutrino mass m i and depends on the nucleus that undergoes the 0νββ transition. It can be satisfactorily approximated by the analytic expression where p 2 −(125 MeV) 2 is the virtual momentum of the neutrino. We will conduct a detailed analysis of the impact of two minimal ISS realisations, the "(2,2) ISS" and "(2,3) ISS", on the effective electron neutrino mass in Sections 4.1.4 and 4.2.5.
Moreover, if the typical scale of the new states is sufficiently light, they can be produced in collider or low-energy experiments, thus being subject to further constraints [31]. Robust laboratory bounds arise from direct searches for sterile neutrinos, which can be produced in meson decays such as π ± → µ ± ν, with rates that depend on their mixing with the active neutrinos. Therefore, negative searches for monochromatic lines in the muon spectrum can be translated into bounds on the active-sterile mixing [7,32].
All the above mentioned bounds will be taken into account in our subsequent numerical analysis of the two minimal ISS realisations.

Phenomenological analysis
Although it is possible to derive analytical expressions for the neutrino mass eigenvalues and leptonic mixing matrix (see Appendices A and B), these expressions are lengthy and involved, and do not easily convey the general features and behaviour of the ISS configurations investigated. We thus conduct a numerical analysis for each of the minimal "(2,2) ISS" and "(2,3) ISS" realisations. In order to unveil some interesting features, we performed a scan of the parameter space (corresponding to all the entries of the mass matrix; in our analysis we will not address the effect of CP violating phases, both Dirac and Majorana). This also allows to numerically compute interesting quantities, as for instance the effective mass in 0ν2β decay amplitude. All the constraints listed in Section 3 were implemented. We proceed to discuss the results in the following sections.

The "(2,2) ISS" realisation
Some aspects of this model have already been studied, in particular CP violation and Non Standard Interactions [33]. We have determined the neutrino spectrum and the leptonic mixing matrix using a perturbative approach, whose details are summarised in Appendix B. At second order in the perturbative expansion, the light neutrino spectrum is given by: where the parameters b and c are defined in terms of the entries of the (2,2) mass matrix; these expressions are lengthy, as explained in Appendix B. Notice that b and c do not depend on the submatrix m of the mass matrix of Eq. (9). Having one massless eigenstate (to all orders in perturbation theory) is a feature of this minimal "(2,2) ISS" model (see also Table 1). The expressions of Eq. (27) allow to easily understand why the "(2,2) ISS" model strongly prefers the normal hierarchy scheme. In order to accommodate an inverted hierarchy, i.e. m 2 2 m 2 3 10 −3 eV 2 and m 2 3 − m 2 2 10 −5 eV 2 , one would be led to comply with 10 −6 eV 4 + 4c 10 −10 eV 4 . This amounts to an extreme fine-tuning. Although some solutions can indeed be found (see the numerical studies of the following section), it should be stressed that accommodating a NH spectrum also requires a certain amount of fine-tuning.
Even if useful when addressing the issue of the hierarchy of the light neutrino spectrum, the analytical expressions we have derived for the neutrino masses and leptonic mixings cannot be used to extract general features, nor to infer the magnitude of the fundamental scales of the ISS model (i.e. the magnitude of the entries of the matrices µ, m, ...). To do so, we performed numerical scans of the "(2,2) ISS" parameter space, the result of which we proceed to report.

Mass hierarchy
As discussed in Section 2 and illustrated in Fig. 1, the low-energy phenomenology of a "(2,2) ISS" model is determined by two scales: that of the LNV parameter µ, and the ratio k between the magnitude of the entries of the d and n matrices, see Appendix A.
In our numerical analysis, we randomly scan over all parameters: the entries of the d and n submatrices are varied such that the obtained mixing matrix U PMNS is in agreement with oscillation data (global fits to both hierarchies, normal and inverted [1]) and the interval of variation for the entries of µ is chosen to ensure that the largest neutrino squared mass value ∼ 2.4 × 10 −3 eV 2 . We take all parameters to be real (leading to vanishing Dirac and Majorana phases, and hence no contributions to leptonic electric dipole moments).
The best fit values for the mass eigenvalues resulting from the global analysis of the oscillation experiments [1] are indicated in Fig. 2 by horizontal and vertical lines. This example clearly illustrates that the "(2,2) ISS" model favours a normal hierarchical scheme: as can be seen on the right panel of Fig. 2, no solutions can be encountered for an IH scheme (corresponding to ∆m 2 32 ∼ 10 −5 eV 2 together with m 2 2 ∼ m 2 3 ∼ 10 −3 eV 2 ). Moreover, for the NH scheme, finding solutions for the light neutrino masses in agreement with data -although possible -proves to be very difficult.

Constraints from unitarity
The non-observation of NSI in the leptonic sector as induced by the deviation from unitarity of the U PMNS matrix due to the presence of additional fermions puts stringent constraints [21] on the ISS parameter space.
The non-unitarity effects can be quantified by where N is the 3 × 3 submatrix encoding the mixings between the active neutrinos and the charged leptons, i.e. the PMNS matrix. Depending on the mass regime for the sterile fermions (above or below the EW scale) the constraints on N N † are different [21]. We thus identify the following mass regimes for our sample of "(2,2) ISS" mass matrices: • no (or only some) sterile states are above 1 GeV -implying that not all the extra states can be indeed integrated out; the NSI constraints of [21] do not apply in this case; • all sterile states are heavier than 1 GeV, but do not necessarily lie above the EW scale, Λ EW ∼174 GeV; • all sterile states are heavier than Λ EW .
When appropriate, we thus compute the amount of non-unitarity from Eq. (28), and apply the corresponding bounds, to further constrain the ISS parameter space. Notice that in the ISS models the non-unitarity effects are proportional to the ratio O(d)/O(n) (see for example the neutrino mass eigenvector expression for the one-generation model (Eq. A.11)).
We display on Fig. 3 the most constraining deviations from unitarity parametrised by = |1− N N † |, see Eq. (28), as a function of an effective factor k generalizing the one introduced in Section 2.3, which is defined as (see Eq. (B.4) in Appendix B): Each point is generated with random values for the entries of the d, n submatrices -but allowing the entries of each submatrix to vary at most over two orders of magnitude -, and such that the mass matrix would generate a PMNS matrix and a neutrino mass spectrum in agreement with experimental constraints (in the NH scheme). This would correspond to a tiny sub-sample of the points presented in Fig. 2 (close to the intersection of the green lines).  (29)). On the left, 22 , for a mass regime in which the sterile neutrino masses are between 1 GeV and Λ EW ; on the right, 12 , in the regime where all sterile states are heavier than Λ EW . The green lines indicate the corresponding upper bounds [21]. All points comply with oscillation data in the NH scheme.
As can be seen from both panels of Fig. 3, NSI constraints significantly reduce the number of otherwise phenomenologically viable solutions for the "(2,2) ISS" model.

Lepton number violating parameter space
From the numerous numerical scans we conducted, certain features of the "(2,2) ISS" model became apparent: • Low-energy neutrino data (i.e. masses and mixings) can be accommodated if the entries in each of the submatrices of Eq. (9) are allowed a strong hierarchy -varying at least over 2 orders of magnitude.
• The model leads to a strongly hierarchical light neutrino mass spectrum, with the second lightest neutrino mass being strongly suppressed with respect to the heaviest one (the first state being massless).
The size of the LNV parameters (i.e. the entries of the µ submatrix -recall that the LNV matrix m does not enter in the expression for the lightest neutrino mass eigenvalues, as derived in a perturbative approach -see for instance, Eq. (7)) is bounded from below by PMNS matrix constraints, and from above by the naturalness requirement. The lower limit is due to the fact that, to a good approximation, the entries of d must be at least one order of magnitude smaller than those of n (in order to accommodate oscillation data). In order to fulfill solar and atmospheric mass squared differences, and given that one typically has k < 10 −1 (see Eq. (7)), it follows that We have checked that the latter condition is indeed valid in the "(2,2) ISS" model; the lower values for the µ submatrix entries, in agreement with both U PMNS data and neutrino mass squared differences are: min |µ i,i | ∼ 0.13 eV, min |µ i =j | ∼ 5 × 10 −6 eV. The upper bound on the LNV parameters comes from 't Hooft naturalness criterium, even though a clear definition regarding the naturalness of a small dimensionful parameter breaking some SM accidental symmetries does not exist. In this study, we have posited a "naturalness" upper limit of 100 eV on the entries of the submatrix µ. This translates into a lower bound on the factor k (since m ν ≈ k 2 µ).

Neutrinoless double beta decay
When applied to the "(2,2) ISS" model, the effective neutrino mass m νe eff determining the amplitude of the neutrinoless double beta decay rate is given by (see Section 3) [30]: where p 2 −(125 MeV) 2 is the virtual momentum of the neutrino. From the analytical expressions derived in Appendix B.2, one can see that in the limit µ i,j , m i,j → 0, one has m 5 → m 4 , m 7 → m 6 , U 2 e,4 → U 2 e,5 , U 2 e,6 → U 2 e,7 , and thus the extra contribution vanishes. Our predictions for the effective electron neutrino mass are collected in Fig. 4, and displayed as a function of the mass of the lightest sterile state, m 4 . By defining an "average" effective sterile mass, m s = m 4 +m 5 +m 6 +m 7 4 , three distinct mass regimes for m s can be identified from Fig. 4, • m s |p|: in this regime the effective mass goes to zero, since from Eq. (31) one approximately has and one can write where M denotes the full neutrino mass matrix.
• m s ≈ |p|: the contribution of the pseudo-Dirac states becomes more important, and can induce sizeable effects to m νe eff . • m s |p|: in this regime the heavy states decouple, and the contributions to m νe eff only arise from the 3 light neutrino states.
Notice that the values of m νe eff displayed in Fig. 4 correspond to conservative (maximal) estimations; since in our scan all parameters are taken to be real, no cancellation due to possible (Majorana) phases can take place, and thus reduce the contributions of the "(2,2) ISS" model. It is important to stress that all points leading to Fig. 4 comply with all available low-energy constraints discussed in Section 3. The MEG bound on Br(µ → eγ) [29] and the constraints from laboratory experiments [32] are particularly important, and the latter are in fact responsible for the exclusion of a significant amount of points found (corresponding to the grey regions) in Fig. 4.  [34]; blue points pass all imposed constraints (oscillation data, NSI, Br(µ → eγ) and laboratory direct searches), while grey points are excluded by laboratory bounds.

The "(2,3) ISS" realisation
We now address the phenomenology of the next-to-minimal configuration, the "(2,3) ISS", where two generations of RH neutrinos and three sterile states are added to the SM content. In view of the degree of complexity of the analytical expressions derived for the simpler "(2,2) ISS", in this case we directly base our analysis on a numerical approach.

Allowed mass hierarchies
Concerning the neutrino spectra, the crucial difference of the "(2,2) ISS" and the "(2,3) ISS" configurations is that the latter contains four light states, one being dominantly sterile-like. Its mass typically lies below the GeV (in the analysis we have explored the interval [0, 100] keV for all the entries of the µ submatrix), recall that the four remaining states are heavy, pseudo-Dirac pairs. As can be seen in Table 1, and similar to what occurred for the "(2,2) ISS", the lightest neutrino is also massless in the "(2,3) ISS" configurations. Thus, bounds on squared mass differences also translate into bounds for the masses themselves.
Our study reveals that the "(2,3) ISS" model is not as fine-tuned as the "(2,2) ISS" one. Allowing the entries of each submatrix of Eq. (9) to vary over one order of magnitude leads to abundant solutions in agreement with low-energy neutrino data. Concerning the hierarchy of the light neutrino spectrum, we have verified that both NH and IH spectra are possible in the explored "(2,3) ISS" parameter space, although IH tends to be only marginally allowed, as is illustrated on Fig. 5. For the left panel (NH), the parameters were varied as d i,j ∈ [10 6 , 10 7 ] eV, n i,j ∈ [10 7 , 10 8 ] eV, m i,j , µ i,j ∈ [10 −1 , 10] eV, while leading to the right plot (IH) we considered d i,j ∈ [10 6 , 10 7 ] eV, n i,j ∈ [10 8 , 10 9 ] eV, m i,j , µ i,j ∈ [10, 10 3 ] eV. Figure 5: Squared masses of the active neutrinos for the "(2,3) ISS" model (the lightest neutrino is massless). All points displayed fulfill the experimental constraints on the PMNS entries in the NH (left) and IH (right) schemes. The green lines denote the experimental best fit values [1] in the NH or IH schemes. The scan details are summarised in the text.

Constraints from non-unitarity
Similar to what was previously discussed for the "(2,2) ISS" configuration, the constraints coming from the non-observation of NSI (see Section 3) also apply to "(2,3) ISS" models. We conducted here an analogous study: the formulae and notations are simple generalizations of those introduced in Section 4.1.2, the only difference being that in the present case the index i in Eq. (28) runs over the states that are integrated out ( 1 GeV), i.e., i = 5, . . . , 8. Moreover and since we are interested in a potential "Warm" DM candidate, we consider realisations of the "(2,3) ISS" model in which only the lightest sterile state lies below 100 keV (i.e. µ ∈ [0, 100] keV).
In Figure 6 we display two examples of deviations from unitarity as parametrised by αβ ≡ 8 i=5 U α,i U † i,β as a function of an effective factor k. We notice that the relative density of points in the figure confirms that the "(2,3) ISS" allows for both spectra, although with a clear preference for NH. As in the previous "(2,2) ISS" model, we again verify that NSI constraints significantly reduce the number of viable solutions for a "(2,3) ISS" configuration.

LFV constraints: Br(µ → eγ)
For completeness, we illustrate the contributions of the new sterile states to rare LFV processes, in particular considering Br(µ → eγ), see Eq. (24). In Fig. 7, we display this observable as a function i,β entries, as a function of an effective factor k (generalization of Eq. (29) for the "(2,3) ISS" model). On the left, 22 , for a mass regime in which the sterile neutrino masses are between 1 GeV and Λ EW ; on the right, 12 , in the regime where all sterile states are heavier than Λ EW . The green lines indicate the corresponding upper bounds [21]. Blue (red) points comply with oscillation data in the NH (IH) scheme.
of the mass of the next-to-lightest sterile state, m 5 . The investigated parameter space leads to contributions typically below the future experimental sensitivity. However, for m 5 in the range [10 2 , 10 4 ] GeV, one might observe a cLFV signal of the "(2,3) ISS" at MEG.  [29] (future sensitivity [35]); blue and red points correspond to NH and IH solutions, respectively, and pass all imposed constraints (oscillation data and NSI).

An intermediate sterile scale
A fundamental difference between the "(2,2)" and the "(2,3) ISS" models is that, since in the latter case #s − #ν R = 1 (see Section 2.3), the model has a third intermediate energy scale O(µ), which corresponds to the mass of a sterile state. It follows that if µ ≈ eV this model can accommodate a 3 + 1-scheme that can potentially explain the (anti)-neutrino anomalies in the short baseline, Gallium and reactor experiments. Should µ ≈ keV, then the model can potentially provide a WDM candidate (see for example the analysis of [17]).
In Figure 8 we display the mixings of the light sterile state with ν e , as a function of m 2 4 . All points are in agreement with constraints from oscillation data, NSI, laboratory and LFV constraints. As is clear from Fig. 8, the parameter space of the "(2,3) ISS" can provide solutions to either reactor anomaly. It can also provide a WDM candidate in the form of a sterile state of mass ∼ 1 keV. Figure 8: Mixings between the electron neutrino and the lightest sterile state, as a function of the sterile squared mass m 2 4 . The green lines indicate the best fit values of (∆m 2 41 , |U e4 |) for the 3 + 1-scheme [6], while the purple vertical line indicates the value m 2 4 = (2 keV) 2 , corresponding to the mass of the (warm) dark matter candidate suggested in [17]. Blue and red points correspond to NH and IH solutions, respectively. The points displayed comply with all imposed constraints (oscillation data, laboratory, NSI and Br(µ → eγ)).

Neutrinoless double beta decay
Due to the presence of the extra light sterile state, in the "(2,3) ISS" model there is an additional contribution to the effective mass derived in Eq. (31). In our analysis we assumed the lightest sterile state to have a mass m 4 < 100 keV |p| ≈ 125 MeV, it contributes to the neutrinoless double beta decay effective mass as trivially generalising Eq. (31) and where above, p 2 is again the virtual momentum of the propagating neutrino.
In Figure 9 we summarise our predictions for the effective electron neutrino mass as a function of m 5 . Like in the previous case, by defining an "average" heavy sterile mass m s = m 5 +m 6 +m 7 +m 8 scenario. Especially in regimes of heavier sterile masses (i.e., m 5 1 GeV), the model is fairly predictive regarding the 0ν2β decays: the value of the effective mass in "(2,3) ISS" scenario lies just below the current experimental bound and within the future sensitivity of ongoing experiments [34]. Somewhat lighter sterile masses could also account for an effective mass within experimental reach, but these solutions are already excluded by the recent MEG bound and by laboratory constraints. Figure 9: Effective electron neutrino mass, m νe eff , as a function of m 5 . The green full and dashed horizontal lines denote the current upper bound and the expected future sensitivity [34]; blue and red points correspond to NH and IH solutions, respectively, and pass all imposed constraints (oscillation data, NSI, Br(µ → eγ) and laboratory direct searches), while grey points are excluded by laboratory bounds.

Conclusions and future prospects
In this work we proposed a methodological approach to identify the most minimal Inverse Seesaw realisations fulfilling all phenomenological requirements. By adding extra sterile fermions to the SM (right-handed neutrinos, ν R , and sterile singlets, s) whose number of generations were not fixed (#ν R not necessarily equal to #s), we showed that it is possible to construct several distinct ISS models that can reproduce the correct neutrino mass spectrum.
Our general analysis has shown that the mass spectrum of an ISS model is characterised by either 2 or 3 different mass scales, corresponding to the one of the light active neutrinos, that corresponding to the heavy states, and an intermediate scale associated to #s − #ν R sterile states (only relevant when #s > #ν R ).
The approach we followed was based on time-independent perturbation theory for linear operators, which allowed to diagonalize the neutrino mass matrix analytically. One can thus obtain analytic expressions for the neutrino eigenstates and associated masses as a power series of the small parameters that violate the total lepton number.
As a result, we were able to identify two classes of truly minimal ISS realisations that can successfully account for neutrino data. The first, here denoted "(2,2) ISS" model, corresponds to the SM extended by two RH neutrinos and two sterile states. It leads to a 3-flavour mixing scheme, and requires only two scales (the light neutrino masses, m ν and the RH neutrino masses, M R ). Although considerably fine tuned, this ISS configuration still complies with all phenomenological constraints, and systematically leads to a Normal Hierarchy for the light neutrinos. The model could marginally give rise to an effective mass for 0ν2β within experimental reach, but all these regions turn out to be excluded by current laboratory constraints and MEG bounds on µ → eγ decays.
The second, the "(2,3) ISS" realisation, corresponds to an extension of the SM by two RH neutrinos and three sterile states. This class allows to accommodate both hierarchies for the light spectrum (although the IH is only marginally allowed), in a 3 + 1-mixing scheme. The mass of the lightest sterile neutrino can vary over a large interval: depending on its regime, the "(2,3) ISS" realisation can offer an explanation for the reactor anomaly (in this case, m 4 ∼ eV), or provide a Warm Dark Matter candidate (for a mass of the lightest sterile state around the keV). However, the detailed study of the latter possibility is beyond the scope of this work and requires a complete and comprehensive analysis that will be conducted in a subsequent study. Finally, concerning 0ν2β decays, the "(2,3) ISS" scenario leads to effective masses close to the current experimental bound and within future sensitivity of coming experiments [34].
In this work, we have focused on the determination of the truly minimal inverse seesaw realisations. Our approach can be easily generalised to probe the phenomenological viability and impact of any ISS extension of the SM (for an arbitrary number of RH states and sterile fermions).

A Perturbative determination of the neutrino masses and of the leptonic mixing matrix
In the one generation ISS model, and in the basis defined by n L ≡ (ν L , ν c R , s) T , the neutrino mass matrix can be written as where d, m, n, µ are complex numbers. This symmetric matrix can be diagonalized via [19] where U is a unitary matrix and m 0,1,2 are the physical masses. To obtain U , we use the hermitian so that the matrix U diagonalizing M † M is the same as the one in Eq. (A.2).
In the following, we proceed to diagonalize the one-generation squared mass matrix M † M of Eq. (A.1), using perturbation theory for linear operators. We also discuss the validity of the perturbative approach. The mass matrix M can be decomposed as where M 0 is the zeroth order matrix and ∆M is the perturbation (which violates lepton number by two units). One can write M † M as where M 2 I and M 2 II are the components of the perturbation that are homogenous functions of first and second order in the small parameters m and µ (|m|, |µ| |d|, |n|). The perturbativity condition ||∆M || ||M 0 || translates into conditions for the M 2 0 , M 2 I and M 2 II matrices The perturbative determination of the mass eigenvalues is thus ensuring , provided that |m|, |µ| |n|. 2 , an orthonormal combination of eigenvectors associated to the degenerate eigenvalue |d| 2 + |n| 2 , the first order correction to x (0) 0 is given by Since |µ|, |m| |n|, the coefficients in Eq. (A.8) verify Similar arguments apply to the first order corrections to x (0) j=1,2 ; the second order eigenvector corrections are still subdominant, thus confirming the validity of the perturbative approach.
The lightest neutrino mass arises from perturbative corrections to the m = 0 eigenvalue, while the two other states are massive and degenerate (pseudo-Dirac heavy neutrinos). The correction to m 2 0 (0) at second order is which reduces to the usual inverse seesaw result once the condition |d| |n| is assumed. As discussed in Section 2, in this approach the only assumption on the magnitude of the physical parameters is driven by the naturalness requirement, i.e. |m|, |µ| |d|, |n|. The eigenvector associated to m 2 0 (2) is given at zeroth order in the perturbative expansion by 2 Using Eq. (23), the number n p of physical parameters is 24. In the following we choose 3 a basis in which one has exactly 24 free parameters, as shown in Table 5.

B.1 Massless eigenstate
Having a massless eigenstate is an unavoidable feature of the minimal "(2,2) ISS" and "(2,3) ISS" realisations. In the minimal "(2,2) ISS" realisation, the massless eigenstate is given by which is compatible with the constraints on the U PMNS matrix, in both cases of normal and inverted hierarchy.