Minimal inverse-seesaw mechanism with Abelian flavour symmetries

We study the phenomenology of the minimal $(2,2)$ inverse-seesaw model supplemented with Abelian flavour symmetries. To ensure maximal predictability, we establish the most restrictive flavour patterns which can be realised by those symmetries. This setup requires adding an extra scalar doublet and two complex scalar singlets to the Standard Model, paving the way to implement spontaneous CP violation. It is shown that such CP-violating effects can be successfully communicated to the lepton sector through couplings of the scalar singlets to the new sterile fermions. The Majorana and Dirac CP phases turn out to be related, and the active-sterile neutrino mixing is determined by the active neutrino masses, mixing angles and CP phases. We investigate the constraints imposed on the model by the current experimental limits on lepton flavour-violating decays, especially those on the branching ratio $BR(\mu\rightarrow e \gamma)$ and the capture rate $CR(\mu-e,{\rm Au})$. The prospects to further test the framework put forward in this work are also discussed in view of the projected sensitivities of future experimental searches sensitive to the presence of heavy sterile neutrinos. Namely, we investigate at which extent upcoming searches for $\mu\rightarrow e \gamma$, $\mu \rightarrow 3e$ and $\mu-e$ conversion in nuclei will be able to test our model, and how complementary will future high-energy collider and beam-dump experiments be in that task.


I. INTRODUCTION
The discovery of neutrino oscillations has not only established that neutrinos are massive particles, but also opened the way to a whole new domain of searches for physics beyond the Standard Model (SM). Although such a milestone represents a major achievement in particle physics, our current understanding of neutrino properties leaves unanswered many questions concerning the fundamental properties of those particles, namely the ones related to the origin of their masses. From a theoretical perspective, the seesaw mechanism [1][2][3][4][5][6] offers an elegant framework for the explanation of neutrino masses and lepton mixing in extensions of the SM. As opposed to the canonical type-I seesaw, where very heavy "right-handed" (RH) neutrinos or tiny Yukawa couplings are required to generate small neutrino masses, in the so-called inverse seesaw (ISS) [7][8][9] neutrino mass suppression is triggered by small lepton-number violating (LNV) mass parameters. In this case, the lightness of neutrinos stems from an approximate lepton-number symmetry which is restored when those parameters are set to zero. Therefore, the ISS provides a natural neutrino-mass generation mechanism in the t'Hooft sense [10].
A crucial feature of the ISS (not shared by the canonical type I seesaw) is that small Majorana neutrino masses can be generated with RH neutrino masses at the TeV scale (or below) and O(1) Yukawa coupling parameters. As a result, the mixing between the (active) light neutrinos and the new (sterile) states can be sizeable for sterile neutrino masses lying not far from the electroweak scale. The presence of new neutral fermions interacting with SM leptons and gauge bosons motivates phenomenological studies beyond the SM, making the ISS a perfect theoretical framework to guide new physics probes. In particular, experimental searches for charged lepton flavour violating (cLFV) processes like µ → eγ [11,12], µ → eee [13] and µ − e conversion in heavy nuclei [14][15][16][17][18] have been studied in the ISS framework [19][20][21][22][23][24][25][26] with the purpose of understanding at which extent our current knowledge on those processes is able to constrain the ISS parameter space. Depending on their masses and mixing with the SM degrees of freedom, sterile neutrinos may also lead to interesting signals potentially observable at the Large Hadron Collider (LHC), as well as at other experiments sensitive to new physics effects induced by the presence of those particles [27][28][29][30][31][32][33].
Global analyses of neutrino oscillation data in a three-neutrino mixing scheme [34][35][36] imply the existence of at least two massive active neutrinos. Yet, controversial results from current neutrino oscillation experiments like LSND [37] and MiniBooNE [38] may be hinting at the existence of sterile neutrinos with masses in the eV range. To accommodate all data, more general active-sterile neutrino mixing schemes would then be required. This inspired several analyses of oscillation data beyond the three-neutrino paradigm (see, for instance, Refs. [39][40][41][42][43]). As turns out, it is possible to construct several ISS low-scale models that are compatible with neutrino oscillation data and, simultaneously, satisfy all phenomenological constraints. In particular, it has been shown that the minimal ISS realisation corresponds to extending the SM with two RH neutrinos and two sterile singlet fermions [20], to which we will refer as the ISS(2,2) model.
A longstanding and challenging issue in particle physics is the lack of a guiding principle to explain the flavour arXiv:2012.04557v2 [hep-ph] 21 May 2021 structure of the SM, i.e., the observed fermion mass spectra and mixing patterns. This flavour puzzle provides a strong motivation for building models with additional particle content and extended continuous and/or discrete symmetries. Once such symmetries are explicitly or spontaneously broken, they will lead to the required fermion mass and mixing structures. Several frameworks have been put forward to tackle this puzzle [44]. One of the simplest approaches consists on the implementation of texture zeros in the Yukawa coupling and mass matrices, imposed by continuous U(1) and/or discrete Z N transformations (see, for instance, Refs. [45][46][47][48][49][50][51][52][53][54][55][56][57]). In the SM extended with RH neutrinos, the realisation of texture zeros with such symmetries is not compatible with data since, in general, they lead to massless charged leptons, massless neutrinos or vanishing lepton mixing angles [57,58]. This is due to the fact that all fermions in the SM couple to the same Higgs field. Thus, enlarging the Higgs sector is a viable solution to surmount this difficulty, being the two-Higgs doublet model (2HDM) [59] the most economical one. Inspired by the above ideas, in this work we consider the ISS(2,2) within the 2HDM supplemented with Abelian symmetries to ensure maximal predictability, i.e., to impose the most constraining flavour structure, so that the charged-lepton masses and current neutrino data can be accommodated, while fulfilling all relevant phenomenological constraints. This can be realised by adding to the scalar sector of the SM another scalar doublet and two complex scalar singlets which, upon spontaneous symmetry breaking, generate all relevant mass terms required to implement the ISS (2,2). Moreover, we will show that CP can be spontaneously broken by the complex vacuum expectation value (VEV) of one of the singlets, and that such CP violation (CPV) can be communicated to the neutrino sector via neutrino-scalar interactions.
This paper is organised as follows. In Section II, general aspects of the ISS mechanism are reviewed paying special attention to the comparison between the effective and full treatment of neutrino masses and mixing. The most restrictive flavour structures for the mass matrices in the ISS(2,2) framework are then identified in Section III by performing a systematic search of all possible texture-zero combinations leading to low-energy neutrino parameters compatible with global analyses of neutrino oscillation data. After setting the successful cases, in Section IV we select those which can be realised by Abelian horizontal symmetries. The phenomenological analysis starts in Section V, where spontaneous CPV (SCPV) is considered and the relation between the Dirac and Majorana phases is established in light of present neutrino data. Predictions for the effective neutrino mass parameter relevant for neutrinoless double beta decay are also discussed in that section. The impact of radiative corrections on light-neutrino masses is analysed in Section VI, while the constraints imposed by cLFV decays on the model parameter space are investigated in Section VII. Possibilities of testing the ISS(2,2) with Abelian flavour symmetries at other experiments as, for instance, the LHC, future colliders, beam-dump experiments and cLFV searches are discussed in Section VIII. Finally, our concluding remarks are presented in Section IX. Details regarding the scalar sector and the computations of cLFV decay rates are collected in the appendices.

II. GENERAL ASPECTS OF THE INVERSE SEESAW MECHANISM
The ISS mechanism can be implemented by extending the SM particle content with n R RH neutrinos ν R and n s sterile fermion singlets s, leading to what we denote as ISS(n R , n s ). 1 In this framework, the generic mass Lagrangian for leptons is given in the flavour basis by where ν L = (ν eL , ν µL , ν τ L ) T , ν R = (ν R1 , ... , ν Rn R ) T , s = (s 1 , ... , s ns ) T . For a fermion field ψ we have ψ c ≡ Cψ T with C denoting the charge conjugation matrix. In the above equation, M is the 3 × 3 charged-lepton mass matrix, M D is a 3 × n R Dirac-type mass matrix, M R is a n R × n s matrix, and M s is a LNV n s × n s Majorana mass matrix. The latter can be naturally small in the t'Hooft [10] sense, since lepton number conservation is restored in the limit where the last term in Eq. (1) vanishes. Defining N L = (ν L , ν c R , s) T of dimension n f = 3 + n R + n s , we can write L mass in the compact form where M is the full n f × n f neutrino mass matrix. The charged-lepton mass matrix is bidiagonalised through the unitary transformations e L → V L e L , e R → V R e R , so that 3 with m e,µ,τ denoting the physical charged-lepton masses. For a given M , V L and V R are determined by diagonalising the Hermitian matrices H = M M † and H = M † M as The weak-basis states N L,R are related to the mass eigenstates (ν 1 , ..., ν n f ) T by a n f × n f unitary matrix U such that the full neutrino mass matrix is diagonalised as where m 1,...,n f are the n f (real and positive) Majorana neutrino masses. Notice that, in general, the light-active (heavy-sterile) neutrino masses are labelled as m 1,2,3 (m 4,...,n f ). In the ISS approximation limit (M s , M D M R ), the neutrino mass matrix M of Eq. (2) can be block-diagonalised by writing it in the form The full unitary matrix U of Eq. (6) can then be parametrised as [60] U = so that where U ν and U s are 3 × 3 and (n R + n s ) × (n R + n s ) unitary matrices, respectively; M eff and M heavy are the effective light and heavy-neutrino mass matrices. At leading order, M heavy M R yielding n R + n s heavy neutrinos. In Eq. (8), F is a 3 × (n R + n s ) matrix given at first order in the seesaw approximation by This leads to the 3 × 3 effective light-neutrino mass matrix which can be diagonalised through a unitary rotation of the active neutrino fields, ν L → U ν ν L , satisfying wherem 1,2,3 are the real and positive light neutrino masses in the ISS approximation. The unitary matrix U ν is obtained from the diagonalisation of the Hermitian matrix H eff = M eff M † eff , yielding the unitary lepton mixing matrix after performing the rotation to the charged-lepton mass basis. In general, for massive Majorana neutrinos, U can be parametrised by three mixing angles θ 12 , θ 23 , and θ 13 , and three CP-violating phases: the Dirac-type phase δ and two Majorana-type phases α 21 and α 31 [61],  where c ij ≡ cos θ ij and s ij ≡ sin θ ij . Several neutrino oscillation experiments have been constraining neutrino mass and mixing parameters, namely ∆m 2 21 = m 2 2 − m 2 1 , ∆m 2 31 = m 2 3 − m 2 1 , θ 12 , θ 23 , θ 13 and δ. We present in Table I the results obtained from the most recent global fit of neutrino oscillation parameters [35]. Both mass orderings are considered: normal ordering (NO) where m 1 < m 2 < m 3 , and inverted ordering (IO) where m 3 < m 1 < m 2 . Since no experimental information has been obtained so far on CPV Majorana phases, they remain unconstrained in our analysis. Notice that within the three-neutrino paradigm, the lepton sector is described by a total of twelve parameters: three charged lepton masses, three light neutrino masses, three mixing angles and three phases. We also remark that the existence of a massless neutrino is presently not excluded by data. In such case, there is only one physical Majorana phase and, thus, the total number of physical parameters in the lepton sector is reduced to ten.
We now characterise active and active-sterile mixing considering the full mixing matrix U or, more specifically, the rectangular 3 × n f matrix W αj ≡ U αj (α = e, µ, τ , j = 1, . . . , n f ) which, according to Eq. (8), can be decomposed in the form where W ν and W s are 3 × 3 and 3 × (n R + n s ) matrices, respectively. It is clear that V † L W s defines the mixing between the three active neutrinos and the n R + n s sterile states in the physical charged-lepton basis. Due to the additional fermion states, active-neutrino mixing is determined by the non-unitary matrix [62] where U is the unitary mixing matrix given in Eq. (14) and η is an Hermitian matrix encoding deviations from unitarity of U. Expanding Eq. (8) up to second order in F, one has √ 1 − FF † 1 − 1 2 FF † which, together with Eqs. (10) and (17), leads to Active-sterile neutrino mixing is described by W s given in Eq. (16), which at first order in F is in the basis where M is diagonal. For the analysis that follows, it is convenient to define the following matrices [63], which obey the equalities Note that the mixing between the light and sterile neutrinos is given by the matrix elements B αj for α = e, µ, τ and j = 4, . . . , n f , in the charged-lepton physical basis. Furthermore, the parameters η αβ encoding deviations from unitarity [see Eq. (17)] can be expressed in terms of B through the relation where we used the first relation in Eq. (21) in order to write η αβ solely in terms of the active-sterile mixing.

III. MAXIMALLY-RESTRICTIVE TEXTURES FOR LEPTONS
In this section, we identify the maximally-restrictive textures for the set of matrices (M , M D , M R , M s ) compatible with neutrino oscillation data within the minimal ISS(2,2) framework, 2 where two ν R and two s fermion singlets are added to the SM particle content, i.e. n R = n s = 2 and n f = 7. Our texture-zero analysis is performed assuming the seesaw approximation given in Eq. (11). Later on, we will comment on the validity of this approximation when comparing with the results obtained with the full neutrino mass matrix M. The identification of the compatible textures is based on a standard χ 2 -analysis, using the function where x denotes the input parameters, i.e., the matrix elements of M , M D , M R and M s ; P i (x) is the model prediction for a given observable with best-fit (b.f.) value O i , and σ i denotes its 1σ experimental uncertainty. In our search for viable sets (M , M D , M R , M s ), we require the charged-lepton masses to be at their central values [64], such that the χ 2 -function is minimised only with respect to the six neutrino observables, namely the two neutrino mass squared differences ∆m 2 21 , ∆m 2 31 , the three mixing angles θ 12 , θ 23 , θ 13 and the Dirac CPV phase δ, using the current data reported in Table I [35]. Notice that, in the ISS(2,2) framework, there is always a massless neutrino (m 1 = 0 for NO orm 3 = 0 for IO).
For a given set of input matrices, we consider compatibility with data if the deviation of each neutrino observable from its experimental value is at most 3σ at the χ 2 -minimum [50,57,65,66]. If this is the case, we also test the compatibility of the textures at 1σ. For the sake of simplicity, we shall use the following sequential notation to label the position of the matrix elements of a given 3 × 2 and 2 × 2 texture T, respectively, where we denote the position of any vanishing element labelled i with a subscript, i.e., T i . For instance, in this notation, a matrix with vanishing 11 and 22 elements would be labelled as T 14 .

IV. ABELIAN SYMMETRY REALISATION OF COMPATIBLE TEXTURES
We start this section by specifying the scalar sector of the model. As mentioned before, maximally-restrictive texture zeros in Yukawa coupling matrices cannot be implemented in the SM with Abelian symmetries, since all fermion fields couple to the same Higgs doublet. Hence, to realise such textures, our minimal setup will require the presence of at least two Higgs doublets Φ a (a = 1, 2). Furthermore, to avoid bare mass terms in the Lagrangian, we also add two complex scalar fields S a (a = 1, 2), so that M s and M R are dynamically generated through couplings of S 1 and S 2 with s T Cs and ν R s, respectively. We parameterise Φ a and S a as where v a and u a are the VEVs of the neutral components of Higgs doublets φ 0 a and the scalar singlet fields, respectively. Note that only the phase difference θ = θ 2 − θ 1 is physical (a more detailed analysis of the scalar sector can be found in Appendix A).
Given the minimal fermion and scalar contents described above, the Yukawa Lagrangian relevant for our work is 3 Upon spontaneous symmetry breaking, the scalar fields acquire non-zero VEVs and the above Yukawa interactions yield the generic mass Lagrangian of Eq. (1) for the ISS (2,2). The corresponding mass matrices are given by To implement Abelian flavour symmetries, we require the full Lagrangian to be invariant under the field transformations where, for each field component F , X F denotes a phase of the form e ix F . This invariance requirement yields the following constraints on the Yukawa matrices of Eq. (29): which can be translated into relations among the various field-transformation phases x F . To implement a texture-zero entry in one of the above Yukawa matrices we require that the corresponding phase relation is not fulfilled. We now proceed to identify which of the maximally-restrictive texture sets compatible with neutrino data (see previous section) can be realised by imposing discrete or continuous Abelian symmetries. To this end, we will apply two methods that complement each other, namely the canonical [68,69] and the Smith normal form (SNF) [70,71] methods. Our methodology follows closely the one employed in Refs. [57,66]. We start with the canonical approach applied to the maximally-restrictive textures presented in Table II to reduce the scope of realisable textures before employing the SNF method. We recall that the charged-lepton textures 4 1 and 4 2 cannot be realised through Abelian symmetries in the 2HDM [57]. For the remaining cases, we first write all possible decompositions of the mass matrix textures into the corresponding two Yukawa matrices defined in Eq. (30). Afterwards, for a given (M , M D , M R , M s ) combination, and for all decompositions of its matrices, we solve the corresponding system of algebraic relations for the field phases (or charges) stemming from Eq. (32). If a solution exists for a set of charges, then that specific (M , M D , M R , M s ) is realisable by Abelian symmetries with the fields carrying those charges. In Table V, we present the realisable mass matrix textures and their corresponding Yukawa decompositions, respectively. Notice that although in some cases two decompositions are possible for a given mass matrix, only one is realisable. We set the ordering for Y 1,2 and Y 1,2 D as the one given in Table V. Also, we use the notation R is forbidden by the symmetries as we shall see promptly. Hence, since for all realisable cases M R and M s are fixed by the textures T 14 and T 23 , respectively, from now on we will refer to each case just through the pair notation (M , M D ).
Applying the SNF method to the texture sets passing the canonical method test, we find that the minimal Abelian symmetry group G realising such textures is Irrespective of the type of Yukawa textures, the Lagrangian is invariant under the global continuous symmetry U(1) Y , Y being the SM hypercharge. Since, obviously, that U(1) Y does not impose any texture zero, the actual flavour symmetry group  Furthermore, the Yukawa Lagrangian (29) is also invariant under the following U(1) global symmetry, with the other fields remaining invariant. Although this symmetry does not impose any texture zero on the mass matrices, it restricts the possible coupling terms that can appear in the Lagrangian. Besides the ones given in Eq. (29), only a bare Majorana mass term of the form s T Cs is allowed by this symmetry. Nevertheless, it can be shown that the remaining Abelian symmetries forbid such a bare term. The minimal group that discretises this U(1) symmetry is Z 4 , with q 1 = π/2. Thus, the actual flavour symmetry group is The maximally restrictive texture sets 5 1,I , T 45 , 4 3 , T 124 and 4 3 , T 456 are realised by the flavour symmetry for n = 2, while 4 3 , T 136,I and 4 3 , T 146,I are realised by the Abelian symmetry group for n = 4. The corresponding U(1) F charges can be determined through the canonical method, while the discrete group charges are obtained resorting to the SNF method. We present in Table VI, for each texture set, the Abelian symmetry group that realises the set and the associated transformation charges for each field. In all cases, the full texture decomposition is imposed by the U(1) F symmetry alone. The discrete groups, Z 2 or Z 4 , only preserve some of the texture zeros but ultimately fail in imposing them totally. Yet, they are crucial in forbidding the bare Majorana mass term for sterile singlet fermions s. In fact, the U(1) F charges alone only forbid the diagonal elements of the bare mass term, while the charges of the discrete groups forbid the remaining off diagonal elements. Therefore, the Yukawa Lagrangian remains (5 1,I , T45) (4 3 , T124) (4 3 , T456) (4 3 , T136,I) (4 3 , T146,I) restricted to the form given in Eq. (29) where the term with Y 2 R is forbidden and Y R ≡ Y 1 R . As a final comment, let us note that for the realisable cases the U(1) F symmetry can be discretised into a minimal set of charges corresponding to a Z 5 symmetry.

V. LEPTON MASSES, MIXING AND LEPTONIC CPV
In the previous section we have identified which of the maximally-restrictive texture sets can be realised through Abelian symmetries in the ISS(2,2) context. As will become clear later (see Section VII), throughout the rest of this paper we restrict our phenomenological analysis to the combination 5 1,I , T 45 . In this case, the charged-lepton mass matrix can be parameterised as where a i can always be made real by phase field redefinitions, and m 1,2,3 are the charged lepton masses. Note that the charged-lepton state 1 is decoupled from the remaining ones. The unitary matrices V L,R that diagonalise the Hermitian matrices H = M M † and H = M † M are given by Here, we have used the compact notation c L,R ≡ cos θ L,R and s L,R ≡ sin θ L,R with the angles θ L,R given by We consider the three distinct cases of 5 1 1 textures with 1 = e, µ, τ , labelled as 5 e,µ,τ 1 . Since after the diagonalisation of the charged-lepton mass matrix the unitary matrices V L,R are such that the correct mass ordering is obtained as in Eq. (3), we have with As seen in the previous section, for M of the type 5 1 , only the M D , M R and M s matrices of the type T 45 , T 14 and T 23 , respectively, lead to maximally-restricted neutrino mass matrices. In terms of the Yukawa matrices given in Eq. (29), those mass matrices are realised through Abelian symmetries for the decompositions (see Table V) Notice that we have rephased the fields to remove the unphysical phases such that b i , d i and f i are real and only the phases β 1,2 remain in all Yukawa couplings of Eq. (29). From now on, instead of considering the general case of complex Yukawa couplings, we will consider the scenario in which CP is imposed at the Lagrangian level and, thus, β 1,2 = 0. As shown in Appendix A, the scalar potential of the fields Φ and S 1,2 , with the soft breaking of the U(1) × Z 2 × U(1) F symmetry, allows for a CP-violating vacuum configuration with 4 Together with Eqs. (42), these VEVs lead to the mass matrices being the matrix elements defined as where p and q are rescaling adimensional parameters, which will be useful for later discussions. Taking into account Eqs. (11) and (44), the effective neutrino mass matrix in the original symmetry basis reads with all parameters real. Performing the rotation to the charged-lepton mass basis with one of the unitary matrices V L given in Eq. (40), we obtain for the 5 e 1 case: 4 Spontaneous CP violation can also be successfully communicated to the lepton sector for the texture combinations 4 3 , T 124 and 4 3 , T 456 . This is not the case for the pairs 4 3 , T 136,I and 4 3 , T 146,I for which M D has one texture zero per row (see Table V) and the complex phase ξ can be rephased away.
while for 5 µ 1 and 5 τ 1 the permutations P 12 and P 12 P 23 of Eq. (41) have to be applied on the left and right. The above matrix must be matched with the one defined in terms of the physical low-energy parameters according to Eqs. (12)- (15), for which the matrix elements are for both NO and IO neutrino masses. The lepton mixing matrix U is parametrised as in Eq. (15), with α 31 = 0 and α 21 ≡ α (since one neutrino is massless in our framework -see Section III). It is clear from Eq. (47) that there are six relevant effective parameters (x, y, z, w, ξ and θ L ) that determine neutrino masses and mixings. These are to be compared with seven low-energy physical parameters which define M ij , namely three mixing angles θ ij , two neutrino masses and two CPV phases (the Dirac and Majorana phases δ and α, respectively). Thus, there is a relation among the elements of the effective neutrino mass matrix, which results in a correlation between two low-energy parameters. It can be shown that, for the 5 e 1 case, the said relation is while the corresponding ones for the 5 µ 1 and 5 τ 1 cases can be obtained by performing the index replacements (11 → 12, 13 → 23) and (11 → 13, 13 → 33), respectively. For a given set of θ ij and ∆m 2 21,31 values, the above equations establish how the two CP-phases δ and α are correlated. Moreover, we notice that all parameters in Eq. (47) can be expressed in terms of low-energy neutrino observables. Indeed, for 5 e 1 we have The equivalent expressions for 5 µ 1 and 5 τ 1 can be obtained performing the replacements (11 → 22, 13 ↔ 23) and (11 → 33, 13 → 23, 12 → 13, 23 → 12), respectively.
In order to establish numerically how δ and α are related, we vary the mixing angles θ ij and the neutrino masssquared differences ∆m 2 21,31 within their experimental 1σ and 3σ ranges (see Table I), while changing δ from 0 to 2π. Then, for both NO and IO cases, we compute M ij through Eqs. (48) and (49), respectively. The Majorana phase α is obtained by solving Eq. (50) for 5 e,µ,τ 1 , leading to the results presented in Fig. 1, where in the left (right) column we show the allowed regions in the (δ, α) plane for the NO e,µ,τ (IO e,µ,τ ) corresponding to 5 e,µ,τ 1 with a NO (IO) neutrino mass spectrum. The blue (magenta) regions were obtained taking the 1σ (3σ) intervals for θ ij and ∆m 2 21,31 , while the vertical dark (light) grey band marks the current experimentally allowed region for the Dirac CP phase δ at 1σ (3σ). The results show that there is a strong correlation between α and δ. For both NO and IO mass ordering, the plots exhibit an approximate symmetry under the shift δ → δ + π, which is due to the fact that Eq. (50) is nearly invariant under that transformation at zeroth order in the smallest mixing angle θ 13 . Note that the absence of Dirac-type CP violation (δ = 0, π) implies α = kπ (k ∈ Z). This can be confirmed analytically by evaluating the Dirac and Majorana CP weak basis invariants, J CP Dirac and J CP Maj [72], which are both proportional to sin(2ξ). Notice also that a future measurement of δ in the intervals [45 • , 135 • ] and [225 • , 315 • ] would exclude the NO µ and NO τ cases since, in these ranges, Eq. (50) has no solution.
The fact that the Majorana and Dirac CP-violating phases are related brings up an interesting connection between neutrinoless double beta decay (ββ 0ν ) and neutrino oscillations. With one massless neutrino, the ββ 0ν widths are proportional to the effective neutrino mass parameter m ββ that, in the most general case, should only be sensitive to Majorana phases. Namely, for the NO and IO cases [73] one has NO: m ββ = ∆m 2 21 c 2 13 s 2 12 e −2iα + ∆m 2 31 s 2 13 , Given that α is related to δ through Eqs. (50), a dependence of m ββ on δ can be established, as shown in Fig. 2 for the same cases treated in Fig. 1. The present upper limits on m ββ reported by the KamLAND-Zen [74], GERDA [75],  2 31 in the 1σ (blue) and 3σ (magenta) allowed ranges given in Table I. The dark (light) grey vertical band marks the 1σ (3σ) range for δ shown in the same table, while the vertical dashed line is at the δ best-fit value. The left (right) column corresponds to the cases with NO (IO) neutrino mass spectrum and a charged lepton mass matrix of the 5 e 1 (top), 5 µ 1 (middle) and 5 τ 1 (bottom) type. Hereafter, these different possibilities will be labelled as NOe,µ,τ and IOe,µ,τ . (top-right panel) and 5 τ 1 (bottom) cases (the colour codes are the same as in Fig. 1). In each panel we show the NO and IO results in the lower and upper fraction of the vertical scale, respectively. We also show the upper bounds on m ββ reported by the KamLAND-Zen [74], GERDA [75], CUORE [76] and EXO-200 [77] collaborations (the height of the corresponding rectangles represent the uncertainty in the upper bounds). On the right the sensitivity of several future experiments is shown.
CUORE [76] and EXO-200 [77] collaborations are also shown (the height of the bars reflects the uncertainties in the nuclear matrix elements relevant for the computation of the decay rates). We also show the future m ββ sensitivities planned by the projects AMORE II [78], CUPID [79], LEGEND [80], SNO+ I [81], KamLAND2-Zen [74], nEXO [82] and PandaX-III [83]. Our results show that, although m ββ is always below the current bounds (even taking the less conservative limits), several future experiments will be able to probe the whole m ββ for IO neutrino masses. In particular, for the current best-fit values given in Table I, we have m ββ 40 meV, which is on the upper IO region. As usual, future experiments with sensitivities of (1.0 -4.5) meV will be needed to probe the NO regime (see e.g. Refs. [84] and [85] for reviews on future prospects of ββ 0ν experiments).
We will now study the properties of the sterile neutrino sector, particularly its spectrum and mixing with active light neutrinos (heavy-light mixing). Although for our phenomenological analysis a full numerical computation will be carried out, we first provide an analytical insight following the discussion in Section II -see Eqs. (16)- (20). According to the definition of the effective heavy neutrino mass matrix M heavy given in Eq. (9), and taking into account M R and M s from Eq. (44), we have for the 4 × 4 unitary matrix U s : where c 1,2 ≡ cos ϕ 1,2 , s 1,2 ≡ sin ϕ 1,2 and  In the ISS approximation, ϕ 1 ϕ 2 π/4 leading to c 1 s 1 c 2 s 2 1/ √ 2, which characterises the two pairs of pseudo-Dirac neutrinos with masses 5 and mass differencesm 5 −m 4 = µ s ,m 7 −m 6 = pµ s . As expected, the degree of degeneracy between the masses is controlled by the small LNV parameter µ s . From Eqs. (19), (40) and (54), we obtain for the heavy-light neutrino mixing defined in Eq. (20): in the 5 e 1 case. The corresponding results for 5 µ 1 and 5 τ 1 are obtained applying on the left of the above matrix the permutation P 12 and P 12 P 23 of Eq. (41), respectively. Notice that each neutrino in a quasi-Dirac pair couples similarly to each lepton flavour α, i.e. B α4 B α5 and B α6 B α7 . Since we are working in the charged-lepton mass basis, the above matrix provides an approximation to the mixing among the charged lepton e α and the sterile neutrinos ν 4−7 . It turns out that, due to the Abelian symmetries imposed to realise the maximally-restricted textures, the charged-current (CC) flavour mixings for distinct lepton flavours are related by low-energy neutrino parameters. Indeed, for the 5 e 1 case, where all parameters involved depend on the neutrino observables as shown in Eq. (51). The corresponding relations for the 5 µ 1 and 5 τ 1 textures are obtained by performing the replacements (e ↔ µ) and (e ↔ τ ), respectively. In Table VII we show the numerical values of some heavy-light mixing ratios using the best-fit values for the low-energy neutrino parameters given in Table I.

VI. RADIATIVE CORRECTIONS TO NEUTRINO MASSES
The analysis presented in the previous section was based on the assumption that the (tree-level) ISS approximation for the neutrino mass matrix given in Eq. (11) is valid, i.e., the parameters in Eq. (44) are such that µ s , m Di M . However, the presence of new fermions and scalars may induce relevant corrections to the light-neutrino masses that should not be overlooked. This matter becomes even more important when considering the high precision achieved in the determination of the oscillation parameters, which is currently at the level of a few per cent for some of those observables (see Table I). The one-loop radiative corrections to neutrino masses have been computed in several works [86][87][88][89]. Here we revisit the calculation of the one-loop corrections to the light neutrino mass matrix, and adapt it to our case (the one-loop neutrino self-energy diagrams are succinctly presented in Fig. 3). We compute the 3 × 3 one-loop correction matrix δM given by [86,87] where Σ( / p) are the (active) neutrino self-energies, evaluated at / p = 0 since the tree-level light neutrino masses are extremely small. The number of neutral (charged) scalar mass eigenstates S 0 a (S ± a ) is denoted by n 0 − 1 and (n ± − 1), respectively. First, note that the W ± -boson contribution vanishes since the self-energies are evaluated at / p = 0. Furthermore, by using the unitarity of the full n f × n f matrix U defined in Eq. (6), it can be shown that the G ± and S ± a contributions to δM also vanish [88]. Thus, only δM Z , δM G 0 and δM S 0 a will be relevant. Performing the computation in the basis where the tree-level neutrino mass matrix M given in Eq. (6) is diagonal, we obtain The couplings C ij are defined in Eq. (20), m k are the tree-level neutrino masses, and m S 0 a are the neutral scalar masses. The Yukawa couplings matrices ∆ a ν appear in the interaction terms between neutrinos and neutral scalars S 0 a -see Eq. (B6) of Appendix B. The loop function in the above expressions is In our framework, the neutral scalar degrees of freedom originate from the two Higgs doublets Φ 1,2 and the two scalar singlets S 1,2 (see Table VI), which results in a total of seven scalar mass eigenstates (S 1−7 ) (see Appendix A for more details on the scalar sector). In the exact alignment limit with no mixing among scalar doublets and singlets, and following the same argument as before regarding the unitarity of U [88], it can be shown that the S 0 3,4,6,7 contributions to δM vanish. Hence, only the (SM Higgs boson) H 0 , R and I will be relevant, being their contributions given by where the matrix N ν is defined in Eq. (B9). Notice that, assuming m R = m I , the R and I contributions cancel each other. This is what we will assume in our analysis (we refer the reader to Ref. [90] for a model where these neutral scalar contributions are taken into account). As for divergences, those coming from the G 0 and H 0 loops cancel each other, likewise for the R and I scalars. The Z-boson contribution is finite by itself due to the first relation in Eq. (22). We now wish to evaluate how the results of the previous section are affected when performing a numerical analysis to compute neutrino mass and mixing observables, including the loop corrections discussed above. In particular, we are interested in • Comparing the tree-level light neutrino parameters obtained using the seesaw-approximated M eff of Eq. (11) with those stemming from the full neutrino mass matrix M in Eq. (2). This will not only provide an insight about how the results are affected by neglecting higher-order terms in the seesaw expansion, but will also set limits on the parameters, above which the approximation holds up to a certain precision. As a reference observable we choose the neutrino mass-squared difference ∆m 2 31 , which turns out to be the most sensitive one to such corrections. To quantify the effect of considering the ISS approximation at lowest order, we define the parameter where the light-neutrino masses m i andm i are determined using Eqs. (6) and (12), respectively. For instance, ∆ ISS = 0.1 indicates that the value of ∆m 2 31 determined from M differs from that computed with M eff by 10%. • Evaluating the impact of the one-loop δM on the determination of low-energy neutrino parameters. For that, we will consider the neutrino mass squared differences ∆m 2 21 and ∆m 2 31 , since for the mixing angles and CP phases the corrections are negligible when compared with the current experimental precision. Likewise the previous case, we define the parameter where m i are the one-loop corrected neutrino masses and ∆m 2 ij are the tree-level neutrino mass-squared differences computed with the full M. Notice that, since we evaluate δM in the basis where M is diagonal, the light neutrino masses m i are determined diagonalising the matrix  43)] and consider all physical neutral and charged scalar masses to be 1 TeV, except for the SM Higgs boson with mass m H 0 = 125 GeV. Under these premises we span the parameter space in the following way: • The low-energy neutrino parameters are fixed to their best-fit values given in Table I Choosing the initial values M = 100 GeV and µ s = 10 eV, we determine m Di as shown in Eq. (46). In order to probe a wide range of scales, we vary M and µ s in the intervals [1, 10 4 ] GeV and [1, 10 11 ] eV. For an arbitrary pair (M, µ s ), we set the rescaling parameters with respect to the initial values, namely a = M/100 and b = µ s /10. The corresponding m Di are obtained using Eq. (69). This procedure leaves invariant the effective neutrino mass matrix (in the ISS approximation) and guarantees that the low-energy parameters are always those corresponding to the experimental best-fit values. Notice that such procedure may lead to very large Dirac masses m Di . Thus, in order to ensure perturbativity of the Yukawa couplings b i in Eq. (42), we require It is important to stress that rescaling M , µ s and m Di is the only way to probe the parameter space of our model since ratios among different m Di are determined by low-energy parameters which are fixed.
• For each set of (M, µ s , m Di ), the full 7 × 7 neutrino mass matrix M is defined using Eqs. (2) and (44), and then diagonalised as indicated in Eq. (6) to determine U and m 1−7 . The active neutrino mixing is characterised by the 3 × 3 (non-unitary) matrix U of Eq. (17), which is the upper-left 3 × 3 block of U . The non-unitarity effects are parameterised by the matrix η of Eq. (17). Finally, we compute the one-loop corrections to the light neutrino masses as explained above.
We remark that the CC mixing between the charged lepton with flavour α = e, µ, τ and the heavy sterile neutrino ν j (j = 4 − 7) is set by B αj as defined in Eq. (20)-see also Eq. (B1). In practice, B αj is the (α, j) element of U computed above and defined in the charged-lepton mass basis. Notice also that we will be able to cover wide ranges of B αj since these elements scale as m D /M . Throughout the remaining of this work, we will use as reference parameters the average mass of the lightest sterile neutrino pair m 45 , a degeneracy parameter r N and the mixing of the electron with the lightest sterile neutrino V eN , defined as In Fig. 4 we show the contour-level plots for ∆ 21 1L (left panels) and ∆ 32 1L (right panels) defined in Eq. (68) for the NO e (upper panels) and IO e (lower panels) cases. The contours of V 2 eN , r N and ∆ ISS are also shown by solid, dash-dotted and dotted lines, respectively. Within the grey shaded region b max i > 5. From the inspection of these plots we conclude the following: • The validity of the (tree-level) inverse-seesaw approximation is verified at less than the percent level for µ s 10 − 20 eV (region above the ∆ ISS = 1% horizontal line). The reason why ∆ ISS does not depend on m 45 for a given µ s , can be understood taking into account that the next-to-leading-order approximation of the inverse seesaw scales as µ s m 4 D /M 4 and, consequently, (69)), m Di scales as m 45 leaving ∆ ISS invariant. Moreover, for a given m 45 , a rescaling of µ s → bµ s leads to a rescaling of ∆ ISS → ∆ ISS /b, as can be seen in Fig. 4.
• The chosen intervals for µ s and M allows us to swipe a wide range for the heavy-light mixing parameter, namely V 2 eN lies between 10 −13 (or even less) and 10 −3 in the parameter-space region where the inverse seesaw approximation holds up to 1%. The  Fig. 4) since in this case the leading-order ISS approximation for V eN is no longer valid. For the degeneracy parameter r N , the contours follow the linear relation r N µ s /m 45 already shown in Eq. (70) (dash-dotted lines in the right panels of Fig. 4).
• The nearly-vertical levels of ∆ ij 1L indicate that radiative corrections to neutrino masses depend mainly on m 45 . For NO e , the one-loop corrected ∆m 2 21 deviates from the tree-level result by less than 1% for m 45 700 GeV, while for ∆m 2 31 that threshold is at m 45 20 GeV. Instead, the corresponding upper limits on m 45 for IO e are at 3 GeV and 100 GeV, respectively. Notice that, being the one-loop corrections typically larger than the experimental uncertainty (see Table I) for m 45 above those values, this does not mean that the symmetry realisation presented in Section IV is no longer valid. It just signals the fact that the tree-level relations (51) start failing. The take-home message to learn from this analysis is that the scale invariance of the tree-level inverse-seesaw approximation holds up to a certain level depending on the sterile neutrino mass.

VII. CHARGED LEPTON FLAVOUR VIOLATION
Lepton flavour violating processes are in the front line of experimental searches for physics beyond the SM [91]. The ISS model, being a paradigm for low-scale neutrino mass generation, provides a natural scenario for the observation of flavour transitions beyond neutrino oscillations. The charged and neutral current interactions with heavy sterile neutrinos induce new phenomena which are strongly suppressed in the SM extended with light neutrinos only. Nevertheless, the predictive power of the general ISS is limited by the arbitrariness of its parameter values. In the present framework, the symmetries discussed in Section IV provide the ground for a testable scenario in the light of present and future experimental probes on LFV processes. Thus, we now study LFV in the context of the scenarios set in the previous sections. Our attention will be focused on the cLFV processes listed in Table VIII, where the present experimental upper limits and future sensitivities for the BRs and CRs are shown. We have revisited the computation of the rates taking into account the field content in our framework. The details are shown in Appendices B-D.
As stated in Section V, we focus our study on the case of the charged lepton matrix texture 5 1 . To justify this choice, let us look at the cLFV decay − α → − β + γ − δ , which is mediated at tree level by the neutral scalars R and I coming from the two Higgs doublets (in the alignment limit). The corresponding BR for this process reads [57] where g αβ,γδ being m R,I the neutral scalar masses. The matrix N e dictates the interactions between charged leptons and the neutral scalars as shown in Eq. (B13). For the 5 1 case, the structure of N e is The presence of zeros in N e imposed by the flavour symmetry leads to a natural suppression of the above BRs. In fact, considering the most constraining three-body decay µ → 3e (see Table VIII), the tree-level contribution vanishes for the textures 5 e,µ 1 (decoupled electron or muon) irrespective of m R,I . Although this does not hold for 5 τ 1 , in this case the BRs are strongly suppressed by the tiny couplings in the µ − e sector. The same conclusion cannot be drawn for the texture 4 3 since the flavour symmetry does not yield the charged-lepton decoupling feature. To suppress the decay rates in this case, not only very large m R,I masses are required but also they must be extremely fine tuned [57], as can be readily seen from Eq. (72). Similarly, the analysis of the one-loop contribution of the neutral scalars R and I to the decay α → β γ reveals that in the 4 3 case requires fine-tuned scalar masses. On the other hand, for 5 e,µ  Table VIII).
In the blue shaded region CR(µ−e, Ti) < 10 −18 . Limits on b max i and ∆ISS are also shown (grey, green and cyan shaded regions). The results are shown for all cases found to be realisable through Abelian symmetries, i.e. for the 5 e,µ,τ 1 cases with NO (left panels) and IO (right panels) neutrino mass spectra. scalar contribution to the µ → eγ amplitude vanishes.
We now analyse the constraints imposed on the parameter space taking into account the present limits and future sensitivities for the processes indicated in Table VIII. We focus on µ → eγ, µ → 3e and µ − e conversion in Au, Ti and Al, since these are either the most constraining at present or the ones for which the projected sensitivity is higher. Later on, we will comment on other LFV three-body and Z decays. The (m 45 , µ s ) parameter space is covered by means of a full numerical analysis as described in the previous section, using the results given in Appendices B-D. For simplicity, we work in the alignment limit for the two Higgs doublets and in the decoupling limit for the two neutral scalar singlets, as explained in Appendix A. The interaction terms between fermions and scalar eigenstates are those given in Section B 2. Since the analysis of textures and symmetries for the quark sector is out of the scope of this work, the up-quark and down-quark mass matrices are taken to be diagonal, i.e., M u = diag(m u , m c , m t ) and M d = diag(m d , m s , m b ) [64]. Consequently, the CKM matrix V is the identity matrix, and the flavour-changing matrices N d and N u , defined in Eqs. (B10) and (B11), vanish.
The results of our numerical analysis are shown in Fig. 5 for the six cases NO e,µ,τ (left panels) and IO e,µ,τ (right panels). The colour codes in the legend of the upper-left panel apply to the whole figure. By inspecting these results we conclude the following: • The validity of the inverse-seesaw approximation up to 1% level, i.e. ∆ ISS < 1%, imposes lower bounds on the LNV parameter µ s > 10 − 20 eV (cyan shaded regions), which correspond to upper bounds on the mixing V 2 eN 10 −4 − 10 −3 (see Fig. 4). The light (dark) grey regions show that a considerable fraction of the parameter space is excluded if one takes into account b max i < 1 (5) as a perturbativity requirement.
Moreover, the improvement on BR(µ → eγ) foreseen by MEG II (solid orange contour) would have a marginal impact in covering the parameter space in our framework. On the other hand, reaching a sensitivity of BR(µ → 3e) at the 10 −16 level would be more relevant in constraining the parameter space, especially for heavier sterile neutrinos, i.e., for larger m 45 .
• The COMET and PRISM/PRIME projected sensitivities for CR(µ − e, Al) and CR(µ − e, Ti), represented by black and blue dash-dotted contours, respectively, cover a considerable part of the parameter space, leaving unprobed the regions in shaded blue where CR(µ − e, Ti) < 10 −18 . In the best-case scenario (NO e ), probing CR(µ − e, Ti) down to 10 −18 would cover the whole parameter space, as can be seen in the upper left panel.
The above results provide a general idea regarding how present experimental data constrain the minimal ISS with Abelian symmetries, and how future experiments would further probe its parameter space. However, it is interesting to investigate possible relations among different processes and, in particular, to ask whether the observation of a particular cLFV decay would allow us to draw conclusions regarding others. Notice that, in general, this is only possible when there is some relation among the LFV parameters and/or masses, as it is our case (see Table VII). With this purpose, we compare BR(µ → eγ) with BR(τ → α γ), and CR(µ − e, Ti) with BR(µ → 3e). The results are shown in Figs. 6 and 7 for NO and IO, respectively, where all points respect the present limits shown in Table VIII. The colour of each point is linked to the corresponding value of m 45 following the colour map shown in the middle of the figure. Some important conclusions stand out from the observation of these results. Namely, it is clear that any future observation of a radiative (or three-body) τ decay with a BR down to 10 −9 would exclude all scenarios presented in this work. In fact, for both τ → µγ and τ → eγ the BRs are 10 −11 at most, well below the Belle II sensitivity (left panels). In general, the spreading of the scatter points is due to variations in heavy-light mixing and heavy neutrino masses. However, in some cases, we observe a linear relation between BR(τ → µγ) or BR(τ → eγ) and BR(µ → eγ). For instance, using Eqs. (C2) and (C3) together with the approximate results of Table VII, it can be shown that for NO e the relation holds for the whole range of m 45 , in perfect agreement with the numerical results shown in the inner plot of the upper-left panel in Fig. 6. From the comparison of CR(µ − e, Ti) with BR(µ → 3e) (right panels of Figs. 6 and 7) we see that in all cases BR(µ → 3e) is at most 10 −13 . Notice also that, for m 45 100 GeV, there is an approximate linear relation between CR(µ − e, Ti) with BR(µ → 3e), which is no longer valid for higher masses due to cancellations among the various contributions to µ − e conversion amplitudes. We have also investigated whether our framework could lead to observable signals in LVF Z decays by computing BR(Z → α β ) with α = β = e, µ, τ . We have concluded that, after imposing the present constraints on µ → eγ and µ − e conversion in nuclei, the Z → α β rates are well below the future sensitivities given in Table VIII Table VIII). The vertical dashed line corresponds to the Mu3e projected sensitivity for µ → 3e. In all panels, the scatter points obey all constraints shown in Table VIII, being their colour linked to the value of m45 according to the colourmap shown in the middle.

VIII. CONSTRAINTS ON HEAVY STERILE NEUTRINOS AND FUTURE PROSPECTS
In the previous section we have analysed the constraints imposed by LFV experimental searches on the minimal inverse seesaw model with Abelian flavour symmetries, adopting as reference parameters the LNV parameter µ s and the average sterile neutrino mass m 45 . Although this provides a clear understanding on how constrained the original scales in the Lagrangian are, it is convenient to look at the problem from a perspective where µ s is replaced by the active-sterile mixing parameters B αj . Notice that, in our framework, we only need to consider one of these quantities since, as seen in Section V and summarised in Table VII, they are all correlated. From now on, we will take as constrained parameters m 45 and V 2 eN defined in Eq. (70). For each of the scenarios analysed in this work, the correspondence between the (m 45 , µ s ) and (m 45 , V 2 eN ) parameter spaces is established by figures like Fig. 4 (left panels). Following this approach, we will be able to compare the constraining power of the cLFV processes discussed in the previous section with other experimental searches which are usually translated into constraints on mass and mixing parameters (see e.g. Refs. [22,29,31,[103][104][105][106]). We are not interested in carrying out an exhaustive analysis of all sensitivity studies performed so far. Instead, we will consider the following searches: • Beam-dump experiments: In a beam-dump experiment a primary beam strikes a high-density target and produces a large number of secondary heavy mesons which, in the presence of active-sterile mixing, can decay to final states with sterile neutrinos. Part of them fly towards a detector, decaying inside its volume. The NA3 experiment used a 300 GeV π − beam incident on 2 meter long beam dump to, among other purposes, look for the decays of heavy neutrinos N into ν and π + − ( = e, µ) final states. In such experiment, heavy neutrinos may originate from the rare π and K meson decays, or in semi-leptonic decays of charm D and F or beauty B mesons. At NA3 the limit V 2 eN 10 −4 has been set within a heavy neutrino mass region 1 − 2 GeV. The CHARM experiment [108] has conducted similar N → ν searches using a prompt neutrino beam produced by dumping 400 GeV protons on a copper beam dump, setting V 2 eN 10 −7 for 1 − 2 GeV heavy neutrinos. We will consider the NA3 and CHARM exclusion regions reported in Figs. 8a and 2 of Refs. [107] and [108], respectively.

Future:
The beam-dump experiment SHIP will use 400 GeV protons extracted from the CERN Super Proton Synchrotron (SPS) accelerator, dumped on a high-density target. The experiment will search for heavy neutrinos with mass up to ∼ 6 GeV produced in the decays of D and B mesons. As for active-sterile mixing, SHIP will be able to probe it down to V 2 eN ∼ 10 −10 for 1.6 GeV heavy neutrinos. In this work, we will consider the SHIP exclusion region given in Fig. 3 of Ref. [109]. Although designed to measure active neutrino oscillation parameters with high precision [110], the DUNE experiment will also be able to search for heavy sterile neutrinos [111]. This will be achieved by striking a target with a very high-energy proton beam (up to 120 GeV), leading to the production of mainly pions and kaons which may produce sterile neutrinos in their decays. For illustration, we will consider the results presented in Fig. 2 of Ref. [111].
• High-energy colliders: At the Large Electron-Positron (LEP) collider, the L3 and DELPHI collaborations have looked for heavy neutrinos N produced via on-shell Z boson decays e + e − → Z → N ν. Several N decay modes were considered, namely N → Z * ν (Z * → , νν, jj) and N → W * (W * → ν , jj ). The L3 results were able to probe V 2 eN down to 10 −4 for ∼ 20 GeV heavy neutrinos, while the DELPHI collaboration conducted similar searches excluding V 2 eN 3 × 10 −5 for masses in the range 3 − 50 GeV. We will consider the L3 and DELPHI exclusion regions given in Figs. 6 and 10 of Refs. [112] and [113], respectively.
At the LHC, the ATLAS and CMS collaborations are also looking for heavy neutrino signals. Both collaborations search for N production in W ± → ± N followed by subsequent decays N → W ± * ∓ (W ± * → ± ν ) with = e, µ. ATLAS has explored event signatures consisting of three charged leptons (electrons and muons) with same-sign dileptons of the same flavour (LNV mode). CMS has extended the search to include events with lepton number conservation, thus being sensitive to displaced decays. Overall, the ATLAS and CMS analyses on trilepton signatures excluded V 2 eN 10 −4 − 10 −5 for heavy neutrino masses in the 5 − 50 GeV range. Both collaborations have also searched for the decays of heavy neutrinos produced in pp → W ± * → ± N into same-sign dileptons and jets N → W ± ± (W ± → jj ). For the ATLAS exclusion regions we will consider the results given in Figs. 6 and 8 of Refs. [114] and [115], while for the CMS ones we take the results of Figs. 2 and 4 of Refs. [116] and [117], respectively. In the presence of active-sterile neutrino mixing, new interactions of the SM Higgs boson may arise, opening the H 0 → N ν decay channel (if kinematically allowed). The subsequent decays N → W * (W * → ν) and N → Z * (Z * → + − ) at the LHC have been studied in Ref. [118] to constrain the mixing-mass parameter space as shown in Fig. 3 of that reference.
Future: Future high-energy colliders will play a crucial role in searching for heavy sterile neutrinos. In particular, during the high-luminosity LHC phase (HL-LHC), ATLAS and CMS will be able to cover masses up to 2 − 3 TeV. Sensitivity studies have also been performed for a Future Circular Hadron Collider (FCC-hh) at a 100 TeV centre-of-mass energy [30,119]. For the HL-LHC and FCC-hh cases we will consider the exclusion regions given in Fig. 25 of Ref. [119] corresponding to LHC14 and LHC100 with integrated luminosities L = 3 ab −1 and 15 ab −1 , respectively. Heavy-neutrino searches performed at a future high-luminosity e + e − storage ring collider (FCC-ee) can drastically improve the limits on active-sterile mixing down to 10 −11 for ∼ 60 GeV neutrinos (Fig. 8 of Ref. [120] ). In a future e + e − linear collider, as the International Linear Collider (ILC), the sensitivity on heavy-light neutrino mixing can reach values down to 10 −4 for a 500 GeV centre-of-mass energy and an integrated luminosity of 100 fb −1 (Fig. 15  operating at 3 TeV and with L = 1 ab −1 , values of V eN ∼ 10 −5 − 10 −4 can be probed for a 600 GeV − 2.5 TeV mass range (Fig. 24 of Ref. [122]).
Detectors placed near LHC interaction points would allow for searches of heavy-sterile neutrinos produced in pp collisions through the reconstruction of displaced vertices in a low-background environment. Several proposals have been been put forward to conduct this kind of analyses, namely the AL3X [123], CODEX-b [124], FASER2 [125], MATHUSLA [126] and MoEDAL [129] detectors. The sensitivity improvement with respect to that achieved by the main detectors (ATLAS, CMS and LHCb) could be of several orders of magnitude in the low-mass regime. In this work we will consider the FASER2 and MATHUSLA exclusion regions given in Figs. 5 and 37 of Refs. [130] and [126], respectively.
Given that we are dealing with nearly-degenerate sterile neutrinos due to the smallness of µ s , LNV decay modes are expected to be suppressed as a result of the quasi-Dirac nature of the heavy sterile neutrinos. This is, however, not the case if the average decay width is of the order of the mass splitting, i.e. Γ N = (Γ N1 + Γ N2 )/2 ∆m N . This issue is especially relevant when looking for the same-sign dilepton signatures discussed above. In order to provide an insight regarding whether LNV decays are suppressed or not, we will use the same-sign to oppositesign ratio [29,127] such that R ll ≥ 1/3 can be adopted as a criterion to identify the regions of the parameter space where LNV decays are unsuppressed [128]. To compute Γ N = (Γ 4 + Γ 5 )/2 in terms of the sterile neutrino masses and mixings we use the results of Ref. [103].
• Electroweak precision data (EWPD): As already mentioned, in the presence of sterile neutrinos, the active neutrino mixing matrix U relevant for neutrino oscillations [cf. Eq. (17)] is no longer unitary. Deviations from unitarity are constrained by neutrino oscillation data, electroweak precision tests and lepton flavour violating decays [105,[131][132][133][134][135][136][137][138]. In fact, the off-diagonal elements of η αβ defined in Eq. (23) are mainly restricted by the LFV decays studied in the previous section. On the other hand, η αα are restricted by SM gauge boson decays, namely W → α ν α and Z → νν, and universality tests in W and π decays. We will use here the limits for |η αα | obtained in Ref. [105], namely: Notice that, in our framework, the η αβ are not independent since the B αi are related to each other as a result of the Abelian flavour symmetries. These relations are written in terms of low-energy neutrino observables as shown in Eq. (57), implying that η αβ can be expressed by a single mixing parameter which, for convenience of our analysis, we choose to be V 2 eN = |B e4 | 2 . In this case, we have where we have used Eq. (23) and the fact that |B α4 | |B α5 | and |B α6 | |B α7 |. It is then possible to use the above equations together with Table VII to compute the x αi factors and extract upper bounds on V 2 eN from the limits given in Eq. (76) for |η αα | (see Table IX).
• Neutrinoless double beta decay: In the presence of sterile neutrinos, the effective neutrino mass parameter m ββ , relevant for neutrinoless double beta decay, is [139] m ββ where p 2 −(100 MeV) 2 is the virtual momentum of the neutrinos and m i are the physical neutrino masses. The first and second sums run over the number of light and heavy neutrinos which, in the present case, is three and four, respectively. For the 1 GeV-10 TeV mass range studied in this work, the contributions of the second term in Eq. (78) are negligible and, thus, the results in Fig. 2 remain valid (for neutrinoless double beta decay studies in the presence of sterile neutrinos see e.g. Refs. [33,[140][141][142][143][144][145]).
In the left panels of Figs. 8 and 9 we present a summary of all the current constraints discussed above, together with those stemming from µ → eγ (MEG) and µ − e conversion in Au (SINDRUM) searches (see Fig. 6), now shown in the (m 45 , V 2 eN ) plane. For the EWPD exclusion regions we consider the most restrictive V 2 eN limits given in Table IX, i.e. those extracted from |η µµ | (third column). On the right of the same figures, the projected sensitivities of the several experiments enumerated above are shown, including the cLFV ones already presented in Fig. 6 in the (m 45 , µ s ) plane. For all cases, the overlap of the current exclusion regions (left panels) is shown in light yellow. By looking at these two figures one can conclude that: • For m 45 2 GeV, the strongest constraints are typically those imposed by the SINDRUM and MEG limits on BR(µ → eγ) and CR(µ − e, Au), respectively, and by EWPD (left panels). One exception is the IO e case where, for 2 GeV m 45 50 GeV, the DELPHI, ATLAS and CMS limits are stronger. In all situations, the CHARM exclusion region is more constraining when m 45 = 1 − 2 GeV. Notice that the EWPD exclusion regions are not the same for the different NO and IO scenarios since the U(1) flavour symmetries, together with present neutrino data, impose different relations among the B αj . Thus, the limits on |η αα | cannot be directly translated into limits of a single B αj by neglecting the remaining B αk with different k = j.
• Any signal of sterile neutrinos with V 2 eN 10 −4 at future hadron or linear colliders (HL-LHC, FCC-hh, CLIC and ILC regions) would not be compatible with the limits already imposed by current constraints from LFV searches and EWPD (see right panels). Therefore, in the context of the present work, high-energy collider probes conducted at the FCC-ee and at experiments like SHIP, MATHUSLA, DUNE and FASER2 turn out to be of utmost importance (obviously, taking into account the considered sensitivity studies).
• While for NO e cLFV indirect searches are fully complementary to the aforementioned direct ones, this is not the case for the remaining NO and IO scenarios. In particular, for inverted neutrino masses, the region with V 2 eN 10 −9 − 10 −8 cannot be probed by future µ − e conversion experiments. In this case, such mixing regimes can be covered by displaced-vertex experiments and by a high-luminosity Z factory like the FCC-ee. Notice that R ll ≥ 1/3 within the sensitivity regions of those searches (see the brown solid lines in the left panels), indicating that LNV sterile neutrino decays are not suppressed by their quasi-Dirac nature. It should also be mentioned that in the absence of a positive µ → eγ signal, the impact of MEG II data would be mild (compare the yellow regions with the solid orange line in the right panels). Instead, if that decay is observed, ranges for m 45 and V 2 eN can be set in most of the cases, being the latter relatively narrow. As for µ → 3e, future probes conducted by the Mu3e collaboration will be able to cover V 2 eN down to 10 −6 − 10 −7 for wide ranges of sterile neutrino masses.
To conclude this section, we remark that although we have shown the results of our analysis in the (m 45 , V 2 eN ) plane, it is relatively easy to infer how the obtained exclusion regions and sensitivity contours would appear choosing a different mixing parameter by taking into account the relations in Table VII.

IX. CONCLUDING REMARKS
In this paper we have thoroughly investigated the minimal inverse-seesaw mechanism with couplings constrained by U(1) flavour symmetries, and with all fermion masses generated via spontaneous symmetry breaking through VEVs of doublet and singlet scalar fields. After finding the maximally-restrictive mass matrices compatible with present neutrino data, we have identified all possible U(1) symmetry realisations and concluded that at least two Higgs doublets and two complex scalar singlets are required to successfully implement those symmetries. The presence of such singlets opens up the possibility for SCPV, which turns out to be successfully communicated to the lepton sector via their couplings to the new sterile fermions. As a result of SCPV and the Abelian symmetries, the low-energy Majorana and Dirac CP phases turn out to be related to each other. We have also shown that, including one-loop corrections to neutrino masses and requiring them to be at the one percent level, sterile-neutrino mass ranges are established, within which the tree-level results are still valid in light of the present experimental precision in the determination of the oscillation parameters. Due to the flavour symmetries, the heavy-light mixings are not independent, being their ratios entirely determined by the lepton masses, mixing angles and CP phases. This provides a very constrained setup for phenomenological studies in the framework of current and future probes sensitive to the presence of sterile neutrino states. We have studied several cLFV decays and obtained the exclusion regions set by the experimental limits on BR(µ → eγ) and CR(µ − e, Au). These results establish upper bounds on V 2 eN of about 10 −4 − 10 −5 . The prospects to further explore the parameter space were discussed in view of the projected sensitivities of future LFV searches, especially those dedicated to µ → eγ, µ → 3e and µ − e conversion in nuclei.
After analysing the constraining power of cLFV processes, we focused our analysis on alternative probes, namely collider and beam-dump experimental searches that are sensitive to the presence of sterile neutrinos. With minor exceptions, we concluded that the HL-LHC, FCC-hh, ILC and CLIC sensitivity regions are already excluded by current LFV and EWPD constraints for all possible U(1) symmetry realisations. On the other hand, searches at a high-luminosity Z factory as the FCC-ee and at experiments like SHIP, MATHUSLA and FASER2 would be highly complementary to the Mu3e, Mu2e, COMET and PRISM/PRIME projects. Although we have not explored all possible future scenarios which could arise from independent results of different searches, it is clear that a single positive signal in any of those experiments would definitely put at test the scenarios studied in this paper. In this sense, further symmetry-motivated studies performed in the context of sterile neutrino searches are most welcome.
where v a and u a are the VEVs of the Φ a neutral component and of S a , respectively. In the present case, only the phase difference θ = θ 2 − θ 1 is relevant. The above fields transform under the Abelian symmetries as shown in Table VI. The most general potential invariant under those symmetries can be written as, with For reasons that will become clear later, we consider all parameters to be real. To avoid massless Goldstone bosons, we add terms to the scalar potential, which break softly the Abelian flavour symmetries. Such terms could, for instance, originate from the spontaneous breaking of a larger symmetry at very high-energies. Possible soft-breaking terms are with all parameters real. This specific form of V soft is chosen not only to avoid unwanted massless scalars, but also to open up the possibility for SCPV coming from the complex VEV of the scalar singlet S 1 , as will be discussed later. The full scalar potential is then given by V (Φ a , S a ) = V U(1) + V soft . As argued in Section IV, one of the motivations for adding two complex singlet scalars is to account for the mass hierarchy in the inverse-seesaw approximation through u 1 , v 1 , v 2 u 2 . In order to simplify the analysis, we will use that VEV hierarchy to consider the case in which S 2 is decoupled from the remaining scalars and, thus, the quartic term ∼ |S 2 | 4 will dominate over the terms ∼ Φ † a Φ a |S 2 | 2 and ∼ |S 1 | 2 |S 2 | 2 . The analysis is then simplified by taking In order to ensure vacuum stability, the scalar potential has to be bounded from below in any direction of the field space as the fields become large. This can be guaranteed by requiring the Hessian matrix of the quartic couplings in the potential to be copositive [146]. In the case under analysis, this is translated into the following conditions among λ 1−8 [147]: Since we considered all Yukawa and scalar parameters to be real, the full Lagrangian is CP-conserving. Yet, CP can be spontaneously broken if some scalar field acquires a complex VEV. The main motivation for studying this possibility is to provide a dynamical explanation of CPV in the lepton sector manifested through non-trivial Dirac and Majorana phases δ and α, respectively (see Section V). To show that this is indeed possible, we start by minimising with respect to v a , u a , θ, ξ 1 and ξ 2 . In particular, the extrema conditions for θ and ξ 2 lead to θ, ξ 2 = 0, π. Two possible solutions for µ 2 11 , µ 2 22 and µ 2 2 can be then obtained from the minimisation equations of v 1,2 and u 2 , respectively. Namely, Instead, minimising with respect to ξ 1 leads to three solutions: Notice that the last solution in (A9) is the only one which provides a non-trivial phase ξ 1 , leading to the possibility for SCPV. Setting θ, ξ 2 = 0, π in V 0 , the value of the potential at each ξ 1 minimum (V min ) is with Therefore, the CPV solution with ξ 1 = 0, π corresponds to the deepest minimum if µ 4 u 1 < ∓4 √ 2 µ 2 3 .

Scalar mass spectrum
We now briefly discuss the main features of the scalar mass spectrum of our model. In total we have two charged complex scalar fields φ + 1,2 and four neutral ones, φ 0 1,2 and S 1,2 . The mass matrix for the charged scalars is where the upper (lower) sign corresponds to θ = 0 (π). This matrix is diagonalised via the basis transformation where GeV. This leads to the SM massless Goldstone boson G ± and massive charged states H ± with mass The above rotation brings Φ 1,2 into the Higgs basis with H 1,2 [59] given by where Here, H 0 coincides with the 125 GeV SM Higgs in the alignment limit, G ± and G 0 are the charged and neutral Goldstone bosons, and H ± is the physical charged Higgs field. The neutral scalar mass matrix M 2 0 is diagonalised through a unitary transformation T which relates weak and mass eigenstates through As mentioned before, we assume that the mixing of the complex singlet S 2 with the remaining fields is negligible and, thus, η 4 and ρ 4 are decoupled. Furthermore, the CP-odd scalars η 1,2 from φ 0 1,2 also do not mix with the other scalars. However, ρ 1,2,3 and η 3 do mix among themselves. For the SCPV solution in Eq. (A9), the mixing among η 3 , ρ 3 and ρ 1,2 will depend on the soft-breaking parameters µ 3,4 and on VEV products u 2 1 and u 1 v a (a = 1, 2). Moreover, since the naturally small soft-breaking parameters must fulfil the SCPV condition µ 4 u 1 < ∓4 √ 2 µ 2 3 , and given that u 1 v 1 , v 2 , it is reasonable to consider that the mixing among the φ 0 1,2 and S 1 will be small. In this limit, the neutral mass matrix can be recast into a block-diagonal form composed of four 2 × 2 matrices: the CP-even M 2 CP-even for ρ 1,2 , the CP-odd M 2 CP-odd for η 1,2 , and the S 1,2 mass matrices M 2 S1 and M 2 S2 . The former is given by being diagonalised by where, as before, the upper (lower) sign corresponds to θ = 0 (π). The angle α 1 parameterises the mixing in the (ρ 1 , ρ 2 ) sector. Throughout this work we set β = α 1 + π/2, which is known as the alignment limit [148]. In this case, there is no mixing between H 0 and R and, thus, S 0 1 = H 0 and S 0 2 = R. As already mentioned, this allows us to identify H 0 with the 125 GeV Higgs boson discovered at the LHC [149]. For the CP-odd scalars we have which is diagonalised through the rotation matrix defined in Eq. (A13), leading to the massless Goldstone boson G 0 and to the massive scalar S 0 5 = I with mass The S 1 mass matrix is which is diagonalised by a rotation with an angle α 2 given by, . (A23) After diagonalizing the mass matrix M 2 S1 there remains a very small mixing in the (ρ 3 , η 3 ) sector, which depends on the soft-breaking couplings. The mass terms for the resulting physical fields, S 0 3 ρ 3 and S 0 6 η 3 , are given by which are approximately degenerate for the SCPV solution. Notice that in the absence of soft-breaking terms proportional to µ 3 and µ 4 , unwanted massless Goldstone bosons appear. SCPV originated from the complex VEV of the singlet S 1 is possible if the masses above are positive and the condition for the global minimum is satisfied. Lastly, the matrix M 2 S2 is diagonal and the corresponding scalar masses for S 0 4 = ρ 4 and S 0 7 = η 4 are for ξ 2 equal 0 or π, respectively. Once again, if the soft-breaking term proportional to µ 5 vanishes, η 4 would be a massless Goldstone boson. Since u 2 can be naturally large, the scalar fields η 4 and ρ 4 can have a large mass, which further justifies the decoupling behaviour of the singlet S 2 .

Appendix B: Interactions in the mass-eigenstate basis
In this appendix we collect the interactions relevant for our work in the mass-eigenstate basis of fermions, scalars and gauge bosons. We consider an arbitrary number of Majorana neutrinos ν i (i = 1, . . . , n f ) with masses m i , so that the results can be applied for scenarios with any number of sterile neutrinos. In the ISS(2,2) considered in this work n f = 7, being ν 1,2,3 the three light active neutrinos, and ν 4−7 the heavy sterile ones.

Charged-current and neutral-current interactions
In the Feynman-'t Hooft gauge and mass eigenstate basis, the charged-current, weak neutral-current and Goldstone boson interactions read: where we have followed closely the notation of Refs. [63,150]. P L,R = (1 ∓ γ 5 )/2 are the chirality projectors, g is the weak coupling constant, and c W ≡ cos θ W , with θ W being the Weinberg angle. The B and C matrices have been defined in Eq. (20). The last equalities in Eqs. (B2) and (B4) result from the Majorana character of neutrinos (ν = ν c ). Therefore, for Majorana neutrinos, the coupling Z ν i ν j is non-diagonal and involves both chiralities.

Scalar-fermion interactions
In this section we present the scalar-fermion interactions extracted from the Yukawa Lagrangian in Eq. (29), using the notation for the scalar fields introduced in Section A 1. For charged and neutral scalars S ± a and S 0 a we have: respectively, where Γ and ∆ are general Yukawa matrices. In this work, the scalar sector contains two Higgs doublets and two neutral scalar singlets, which we will consider to obey the alignment and decoupling limits discussed in Appendix A. In such case, the Yukawa coupling matrices between the fermions and the charged Higgs S ± 1 = H ± are given by where B is defined in Eq. (20) and The unitary matrices V L,R and U are given in Eqs. (3) and (6), respectively. Additionally, Y 1,2 d,u are the Yukawa matrices appearing in the interaction terms among the Higgs doublets and the quarks fields, and V d,u L,R are the unitary transformations that diagonalise the quark mass matrices. The Cabbibo-Kobayashi-Maskawa (CKM) quark mixing matrix is denoted by V. We explicitly present the interaction Lagrangian involving the H ± charged scalar since it has been extensively used, for instance, in the computation of the cLFV amplitudes (see Appendices C and D). Namely, The Yukawa matrices entering the interaction terms (B6) among leptons and the the neutral scalars S 0 1 = H 0 , S 0 2 = R and S 0 5 = I (those stemming from the Higgs doublets) are where the matrix C is defined in Eq. (20). We explicitly present the Lagrangians involving the above neutral scalars and ν i since they were useful in our calculations of the radiative corrections to light neutrino masses performed in Section VI. We have, Lastly, the coupling matrices appearing in interaction terms between the neutrinos and the neutral scalars S 0 3 ρ 3 , S 0 6 η 3 , S 0 4 = ρ 4 and S 0 7 = η 4 (coming from the complex scalar singlets S 1,2 ) are, where Note that, in order to obtain the Feynman rules using the interactions between ν i and the Z-boson or S 0 neutral scalars, one must multiply by a factor of 2 since Majorana neutrinos are self-conjugate fields.

Appendix C: Charged-lepton flavour violation
In this appendix we present the amplitudes and decay rates for the cLFV processes α → β γ, δ and coherent µ − e conversion in heavy nuclei. The corresponding current experimental bounds and future sensitivities are given in Table VIII, while the various topologies of one-loop Feynman diagrams are summarised in Figs. 10, 11 and 12. We performed the computations in the Feynman-'t Hooft gauge following the interaction Lagrangians in the mass-eigenstate basis given in Appendix B. We provide the results for the amplitudes and branching ratios in terms of the form factors collected in Appendix D.
In Fig. 10 we show the diagrams contributing to the effective vertex β α γ (β = α) at one-loop level. The transition amplitude can be written in the form where α W = g 2 /(4π), ε µ γ is the photon polarisation four-vector and q = p α − p β is the photon momentum. F γ is the local monopole contribution for an off-shell photon (q 2 = 0), while G γ stands for the non-local dipole contribution for an on-shell photon (q 2 = 0). The expressions for F αβ γ,L(R) and G αβ γ,L(R) are given in Section D 1. The former contributes 10. γ and Z-penguin diagrams contributing to radiative decays α → β γ and Z → ± α ∓ β , respectively. Although only the H ± loop contributions are presented, similar diagrams with W ± /G ± were also considered in the calculations. On the right, we show the penguin-diagram topology for − and µ − e conversion in nuclei (ψ1 = ψ2 = u, d).
Although we only display the W ± contributions, similar diagrams with G ± and H ± were also included in the calculations. We have also taken into account LNV diagrams (on the right) due to the Majorana character of neutrinos, and crossed diagrams corresponding to the flavour index interchange β ↔ δ.
In the limit where there is no charged scalar contribution, the branching ratio given in Eq. (C2) is consistent with the results of Refs. [1,[151][152][153]. The one-loop diagrams for the effective vertex β α Z (β = α) are also shown in Fig. 10. In this case, the transition amplitude is where ε µ Z is the Z-boson polarisation four-vector. The branching ratio for the LFV decay Z → − where F αβ Z,L and F αβ Z,R , are given in Section D 2 and the Z-boson total decay width is Γ Z = 2.4952 GeV [64]. Note that, in the limit where the scalar content coincides with the SM one, the above branching ratio is consistent with the results of Refs. [154,155].
δ amplitude receives contributions from the γ-penguin, Z-penguin and leptonic box diagrams shown in Figs. 10 and 11, which we write as where s W = sin θ W , g L = s 2 W − 1/2, g R = s 2 W and Λ X,Y A are given by the following combinations of Dirac matrices and chiral projectors: with σ µν = i [γ µ , γ ν ] /2. From the generic leptonic box diagrams presented in Fig. 11 we obtain the transition amplitude (C8) involving the form factor B αβγδ A,XY , written in terms of the spinorial structure β Λ X A α δ Λ Y A c γ . Since we are in the presence of Majorana neutrinos, LNV diagrams must be also considered, together with cross-diagrams with interchanged lepton indices (β ↔ δ). All these contributions can be written in the form (C8) by using Fierz transformations [156,157]. In Section D 4, the form factors B αβγδ A,XY are presented. In particular, B αβγδ T,LR = B αβγδ T,RL = 0, which can be readily seen from the identity σ µν γ 5 = −iε µνρλ σ ρλ /2.
The calculation of the decay rate for each process is done, whenever possible, assuming vanishing masses for the final lepton states, since m β,γ,δ m α . However, an exception is made for terms involving the photon amplitude, where light lepton masses must be treated with care. As well known, the three-body phase space integrals for the photon contribution are singular in the limit q 2 → 0. Therefore, one must first perform the phase space integration and only after take the limit m β,γ,δ → 0 for non-divergent terms. This will lead to the logarithmic term ln(m 2 α /m 2 β,γ,δ ) appearing in the pure photonic dipole contributions. Also, a symmetry factor of 1/2 has to be taken into account in the case of two identical charged leptons in the final state. The branching ratio for the 3-body LFV decays can be written in the general form: where the values for the branching ratios BR ( α → β ν α ν β ) are given in Eq. (C3). The k i coefficients and the form factors F depend on the charged-lepton final states, for which three combinations are possible. Namely, (i) Two different flavours of leptons, where leptons with the same flavour have opposite charge (β = γ ∧ δ = γ)  Generic u-type (on the left) and d-type (on the right) semi-leptonic box diagrams contributing to the µ − e conversion process. Although for illustrative purposes we only display the W ± -boson diagrams, the G ± and H ± contributions were also included in our calculations.
Our results agree with those obtained in Ref. [158,159] (see also Refs. [160][161][162][163]). We have also checked that the BRs match the well-known results in the limit of a SM scalar content [63]. Following Ref. [164], the operators relevant for coherent µ − e conversion in nuclei have the general form where G F is the Fermi constant [64]. The form factors receive contributions from γ-penguin and Z-penguin diagrams as shown in Fig. 10, the expression for their amplitudes are similar to the ones in Eqs. (C6) and (C7), respectively. The transition amplitude for the semi-leptonic box diagrams in Fig. 12 is given by, leading to, where the values for the nuclear form factors D, S (p) , S (n) , V (p) and V (n) , as well as the muon capture rates Γ capt (Z) are reported in Table X for the nuclei relevant to this work. The scalar-operator coefficients G The above expressions are general and thus can be applied to several models in which µ − e conversion is studied. We have also verified that with only the SM Higgs doublet the above CR coincides with that given in [150].

Photon form factors
The photon form factors coming from the γ-penguin diagrams, generically presented in Fig. 10, are given by where the superscript (1) and (2) correspond to the contributions stemming from W ± and H ± , respectively. Looking at the general form of the amplitude A αβ γ in Eq. (C1), we notice that due to electromagnetic gauge invariance these photon form factors vanish in the limit of zero external lepton momenta and masses. Therefore, to obtain a non-vanishing result we must expand the loop integrals up to next order in q 2 . The W ± -boson contribution yields where B has been defined in Eq. (20). The H ± scalar contribution yields where scalar-lepton couplings N e and N ν are given in Eqs.
The charged Higgs form factors and loop functions given above are consistent with the results of Refs. [167][168][169].

Z-boson form factors
The Z-penguin diagrams displayed in Fig. 10 lead to the form factors: where the W ± -boson contributions are with the matrix C defined in Eq. (20), while the H ± scalar contribution is The loop functions entering the above form factors are given by Z (x, y) = −1 2(x − y) Z (x, y) = 1 4(x − y) The charged Higgs form factors and loop functions given above are consistent with the results of Ref. [167].

Semi-leptonic box form factors
Here we present the form factors and loop functions relevant for the amplitudes (C15) corresponding to the semileptonic u-type and d-type diagrams presented in Fig. 12. We write 7 . We now present the various B αβqq A,XY for the u and d-type diagrams shown in Fig. 12.

Leptonic box form factors
Finally, we present the form factors in A αβγδ Box of Eq. (C8), corresponding to the leptonic one-loop box diagrams generically presented in Fig. 11. As explained in Appendix C, we consider all contributions, including those stemming from LNV and cross diagrams. The form factor B αβγδ V,LL has two types of contributions, coming from the W ± boson and H ± scalar. The remaining form factors are due to box diagrams with H ± . We write where the W ± contribution is Box (λ i , λ j ) . (D20) Box (ω i , ω j ) Box (ω i , ω j ) , Box (ω i , ω j ) , Box (λ i , λ j , λ ± ) Box (ω i , ω j ) , where N e and N ν have been defined in Eqs. (B8) and (B9), respectively. The loop functions relevant for the leptonic box form factors are G Box (x, y) = − √ xy x − y (4 + xy) Box (x, y, z) = − .
The contributions to the box form factors stemming from H ± in the loop are consistent with the derivation performed in Ref. [167]. We refer the reader to Refs. [158,162] where charged-Higgs contributions were considered in the context of low-scale seesaw SUSY. This enabled us to check the types of loop functions and of matrix and chiral structures that enter the resulting form factors. However, to the best of our knowledge, the box form factors and loop functions presented in a compact form in this work within the framework of a seesaw type model and 2HDM scalar sector have not been presented elsewhere.