Phenomenology of the Generalised Scotogenic Model with Fermionic Dark Matter

We study a simple extension of the Standard Model that accounts for neutrino masses and dark matter. The Standard Model is augmented by two Higgs doublets and one Dirac singlet fermion, all charged under a new dark global symmetry. It is a generalised version of the Scotogenic Model with Dirac fermion dark matter. Masses for two neutrinos are generated radiatively at one-loop level. We study the case where the singlet fermion constitutes the dark matter of the Universe. We study in depth the phenomenology of the model, in particular the complementarity between dark matter direct detection and charged lepton flavour violation observables. Due to the strong limits from the latter, dark matter annihilations are suppressed and the relic abundance is set by coannihilations with (and annihilations of) the new scalars if the latter and the Dirac fermion are sufficiently degenerate in mass. We discuss how different ratios of charged lepton flavour violating processes can be used to test the model. We also discuss the detection prospects of the charged scalars at colliders. In some cases these leave ionising tracks and in others have prompt decays, depending on the flavour in the final state and neutrino mass orderings.


Introduction
Neutrino masses and the missing mass in the Universe are among the most important evidence for physics beyond the Standard Model (SM). Some of the most popular proposed explanations for them are radiative neutrino mass models (see Ref. [1] for a recent review on the topic) and particle dark matter (DM), respectively. A simple and elegant candidate of the latter are weakly interacting massive particles (WIMPs). In this work, we study a simple model that has the nice feature of combining both ingredients, thus offering a simultaneous solution to both puzzles. In particular, we study a generalised version of the Scotogenic Model with a dark global U(1) DM symmetry. We denote it the Generalised Scotogenic Model (GScM). A similar model with a gauged U(1) DM has been introduced in Ref. [2] and the original Scotogenic Model (ScM) was proposed by E. Ma in Ref. [3]. In the last years there have been many several studies of the phenomenology of the ScM [4][5][6][7][8][9][10][11][12][13][14]. A systematic study of one-loop neutrino mass models with a viable DM candidate which is stabilised by a Z 2 symmetry has been presented in Ref. [15]. Several variants of the ScM with a U(1) symmetry instead of a Z 2 symmetry have been proposed [16,17] soon after the original ScM model.
The GScM involves two scalar doublets and one Dirac fermion, all charged under the dark global symmetry. In contrast to the models in Refs. [16,17] (and for some variants in Ref. [2]) the U(1) DM symmetry is not broken in the GScM which leads to several changes in the phenomenology of the model. Scalar DM is disfavoured in this model, because of a generically large DM-nucleus cross section mediated by t-channel Z-boson exchange. Thus we will focus on the case of fermionic DM, which in this model is a Dirac fermion, unlike the original ScM. This makes the study of WIMP scattering off nuclei in direct detection experiments very interesting, as it is generated via the DM magnetic moment at one loop. The limits from direct detection experiments already imply the need of coannihilations in the early Universe to explain the observed DM relic abundance.
Neutrino masses are generated at the one-loop level, with a flavour structure different to that involved in processes with charged lepton flavour violation (LFV). The model is very predictive, as the flavour structure of the Yukawa couplings is completely determined by the neutrino oscillation parameters (and the Majorana phase). This allows to draw predictions for LFV processes and decays of the scalars, as we will discuss in depth. The constraints from the non-observation of LFV processes are complementary to the limits from direct detection experiments.
The paper is structured as follows. In Sec. 2 we introduce the GScM and discuss the scalar spectrum and neutrino masses. In Sec. 3 we discuss the most relevant phenomenology of the model, especially charged-lepton flavour violation, the DM abundance as well as collider searches. In Sec. 4 we show the results of a numerical scan of the parameter space of the model. We discuss some other variants of the model (gauge symmetry, Z 3 and Z 4 symmetries) in Sec. 5. A comparison to the original ScM (which involves a Z 2 symmetry) is presented in Sec. 6. Finally, we conclude in Sec. 7. Further details of the model are reported in the final appendices. We discuss in detail the stability of the potential in App. A and fix the notation of neutrino masses and lepton mixing in App. B. The Yukawa parametrisation in terms of neutrino observables is presented in App. C. Some loop functions relevant for different processes are provided in App. D. Expressions for the Electroweak Precision Tests (EWPT) are given in App. E.  Table 1: Particle content and quantum numbers of the GScM. The upper block corresponds to the SM Higgs and the SM leptons, with flavour index α = e, µ, τ . The lower part shows the dark sector of the model: one Dirac fermion ψ and two scalar electroweak doublets Φ and Φ . In the last column we provide the transformation properties under the dark global symmetry.

The Generalised Scotogenic Model
The particle content of the model and its global charges were first outlined in Ref. [2]. It is the generalisation of the ScM in the sense that it is based on a global dark U(1) DM symmetry, while the ScM has a Z 2 symmetry. The SM is augmented by two additional scalar doublets and one vectorlike Dirac fermion, all charged under the U(1) DM symmetry. The particle content and quantum numbers are given in Tab. 1. Without loss of generality we choose the U(1) DM charge of the new particles as q = +1. See Sec. 5 and Ref. [2] for variants of the model and Sec. 6 for a comparison to the ScM. The lightest neutral particle of the dark sector is stabilised by the global U(1) DM symmetry and thus is a good DM candidate. If the DM is identified with the lightest neutral scalar in the dark sector, i.e. the scalar coming from one of the inert doublets, as it carries non-zero hypercharge, neutral current interactions mediated the Z-boson give scattering cross sections off nuclei well above current DM direct detection limits and thus disfavour this possibility. 1 Thus the only viable DM candidate is the SM singlet Dirac fermion ψ whose allowed parameter space we will study in detail. This, indeed, is the most interesting scenario, as there is a connection between DM phenomenology, neutrino masses, charged-lepton flavour violation and searches at colliders. The features of the model, such as neutrino masses and LFV, select the scalar mass spectrum and the possible mechanisms to obtain the right DM abundance.
We denote the SM Higgs doublet by H, which is given in the unitary gauge after electroweak spontaneous symmetry breaking by H ≡ (0, (h + v H )/ √ 2) T , with v H = 246 GeV the vacuum expectation value (VEV) and h the Higgs boson. Without loss of generality we work in the charged lepton mass basis. The Lagrangian for the Dirac fermion ψ is 2 whereL ≡ iσ 2 CL T , with C the charge conjugation matrix, andΦ ≡ iσ 2 Φ * . The neutrino Yukawa couplings y α Φ , y α Φ are three-component vectors In Sec. 2.2 we discuss neutrino masses and estimate the neutrino Yukawa couplings for the case of a neutrino mass spectrum with normal ordering (NO) and inverted ordering (IO). The scalar potential invariant under the U(1) DM symmetry reads: In our numerical analysis we will apply the stability conditions outlined in App. A, which allow for the potential to be bounded from below. Four phases in the Yukawa vectors y Φ and y Φ can be removed by phase redefinitions of the lepton doublets L and the fermion ψ. Similarly, the coupling λ HΦΦ can be chosen real and positive by redefining the scalar doublets Φ or Φ .

Scalar mass spectrum
We assume in the following that the neutral components of Φ and Φ do not take a VEV, so that the global U(1) DM symmetry is unbroken.
The physical scalar states of the theory are given by (i) one real field h, which corresponds to the SM Higgs boson, (ii) two (complex) neutral scalar fields η 0 and η 0 , which are linear combinations of φ 0 and φ 0 (see Tab. 1), and (iii) two charged scalars η ± ≡ φ ± and η ± ≡ φ ± . The Higgs boson mass takes the usual SM value The neutral mass eigenstates are defined as with s θ ≡ sin θ and c θ ≡ cos θ. The mixing angle θ is defined in terms of 2 It is convenient to use the conjugate for the Yukawa of the second doublet, i.e., (y α Φ ) * , so that the expressions for neutrino masses and LFV to be discussed in the following are symmetric under simultaneous interchange of the Yukawa couplings y α Φ ↔ y α Φ and the physical masses (m η 0 , m η ± ) ↔ (mη 0 , m η ± ). There is always a suppression induced by the quartic coupling λ HΦΦ = 0, but it can be further enhanced by a small splitting of the neutral scalar masses difference, m 2 η 0 − m 2 η 0 , which is approximately given by (a − b) 2 for λ HΦΦ 1. The resulting neutrino Majorana mass matrix is of rank two, provided the Yukawa vectors y Φ and y Φ are not proportional to each other. Hence, the neutrino mass spectrum consists of one massless neutrino and two (non-degenerate) Majorana fermions, with masses with |y| denoting the norm of the Yukawa vector y, |y| ≡ α |y α | 2 . Notice that the loop function F (m η 0 , m η 0 , m ψ ) ≥ 0, as m η 0 ≥ m η 0 . The flavour structure is entirely determined by the product of y α Φ y β Φ . As we will see, different combinations enter in LFV processes and DM annihilations, which is a particular feature of this model.
For vanishing solar mass squared difference we can easily estimate the form of the Yukawa couplings y α Φ and y α Φ with the help of the formulae given in App. C. Indeed, from Eq. (85) and taking m 2 = 0, we find both y α Φ and y α Φ to be proportional to the complex conjugate of the third column of the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) mixing matrix u 3 , as shown in Eq. (71). In particular, for neutrino masses with NO we have since θ 23 is around π/4 and θ 13 is small. The Yukawa couplings y e Φ and y e Φ are expected to be smaller in magnitude than the other ones, since they are proportional to θ 13 . Plugged into the formulae for the two non-vanishing neutrino masses m ± ν in Eq. (14), we confirm that m − ν ≈ 0 whereas m + ν does not vanish. For IO we use Eq. (86) with m 1 ≈ m 2 and find that y α Φ and y α Φ are proportional to the sum and difference between the complex conjugate of the first two columns u 1 and u 2 of the PMNS mixing matrix, respectively, i.e.
for c 12 ≡ cos θ 12 , s 12 ≡ sin θ 12 , θ 13 ≈ 0 and θ 23 ≈ π/4. This clearly shows that as well as y e Φ and y e Φ of similar magnitude, but not the same. For the proportionality constant in Eq. (16) being real and positive, as it is assumed in our numerical analysis, we expect the real part of both y e Φ and y e Φ to be positive, since c 12 > s 12 . Furthermore, the imaginary parts of y e Φ and y µ Φ (y τ Φ ) are proportional to each other with a positive (negative) proportionality constant, determined by the ratio s 12 /c 12 . The same holds for the imaginary parts of the Yukawa couplings y e Φ and y µ Φ (y τ Φ ). When plugged into the formula for m ± ν in Eq. (14), we find m − ν ≈ m + ν , as expected for neutrino masses with IO. These expectations for Yukawa couplings y α Φ and y α Φ are confirmed in our numerical analysis to a certain extent, 4 as shown in Figs. 10 and 11 in App. C.

Phenomenology
The most general amplitude for the electromagnetic charged lepton flavour transition α (p) → β (k) γ * (q) can be parameterised as [19] where e > 0 is the proton electric charge, p (k) is the momentum of the initial (final) charged lepton α ( β ), and q = p − k is the momentum of the photon and m α is the mass of the decaying charged lepton α . As is well known, the charged lepton radiative decays are mediated by the electromagnetic dipole transitions in Eq. (18) and are thus described by the form factors A L/R 2 . The monopole which is described by the form factors A L/R 1 does not contribute to processes with an on-shell photon. Thus the corresponding branching ratio (BR) is given by with the fine-structure constant α em = e 2 /(4π) and the Fermi coupling constant G F . The effective couplings in Eq. (18) are radiatively generated at one-loop via γ-penguin diagrams and receive two independent contributions from the charged scalars running in the loop. For the transition − α → − β γ * they are given by where the loop functions f (x) and g(x) are reported in Eq. (87) in App. D. They are approximately equal to 1/6 for small x. The branching ratios BR ( α → β ν α ν β ) are BR (µ → e ν µ ν e ) ≈ 1, BR (τ → e ν τ ν e ) ≈ 0.178 and BR (τ → µ ν τ ν µ ) ≈ 0.174 [20]. In the numerical scan we observe a wide range of values, e.g. the whole region from 10 −32 up to the current bound is allowed for the branching ratio of µ → eγ, depending on the Yukawa couplings and the quartic coupling λ HΦΦ , see Fig. 7. Similar ranges apply for radiative LFV τ -decays, as shown in Fig. 8. 4 For NO the approximation m2 = 0 is oversimplifying, since we neglect in the estimate for y e Φ ( ) the contribution proportional to the second column of the PMNS mixing matrix u * 2 , which is relatively suppressed by (∆m 2 21 /∆m 2 31 ) 1/4 ≈ 0.41 compared to the contribution coming from the third column u * 3 , which is suppressed by θ13 ≈ 0.15.
Notice that in these LFV processes, which are the most stringent ones, there is a different dependence on the neutrino Yukawa couplings y α Φ , y α Φ than in neutrino masses, where the product of both enters, c.f. Eq. (12). This is quite different from the original ScM, where there is only one type of Yukawa interaction, and therefore a very similar combination enters in both neutrino masses and LFV [3,5]. In Sec. 6 we review the structure of the neutrino mass matrix in the ScM.
The branching ratios of the different radiative decays µ → eγ, τ → eγ and τ → µγ are tightly correlated. Using the estimates for the Yukawa couplings y α Φ and y α Φ for neutrino masses with NO and IO given in Sec. 2.2, respectively, we expect that since y µ Φ (y µ Φ ) and y τ Φ (y τ Φ ) are of the same size, whereas y e Φ and y e Φ are suppressed by θ 13 for the case of neutrino masses with NO. For IO we instead expect both of the radiative τ decays to be of similar size and as none of the Yukawa couplings y α Φ and y α Φ is suppressed. Since the experimental bound on BR(µ → e γ) is several orders of magnitude stronger than the one on radiative τ decays, once the constraint BR(µ → eγ) ≤ 2.55 × 10 −13 at 90% CL [21] has been imposed on the parameter space of our model, the branching ratios of the radiative τ decays are automatically below the current limits, i.e. BR(τ → eγ) ≤ 3.3 × 10 −8 and BR(τ → µγ) ≤ 4.4 × 10 −8 at 90% CL [22], and future experimental limits for radiative τ decay, which are sensitive down to BR ≈ O(3 × 10 −9 ) [23].
As the neutrino oscillation parameters are already tightly constrained, it is possible to derive a constraint on the undetermined parameter ζ which affects the relative size of the Yukawa vectors y Φ and y Φ , see Eqs. (85) and (86), as a function of the masses m η ± and m ψ . For large |ζ| the branching ratios are dominated by the diagram with η ± in the loop and for small |ζ| by the diagram with η ± and we can always neglect the other contribution. The loop function f (x) in Eq. (20) takes values between 1/12 and 1/6 in the relevant parameter range 0 < x < 1. After conservatively approximating f ≈ 1/12 and the loop function in the expression for neutrino masses F ≈ 1, we find which, using the Yukawas expressed in terms of mixing angles given in App. C, translates for µ → eγ into the following ranges for any neutrino mass ordering 3 × 10 −4 100 GeV m η ± 100 GeV m ψ sin 2θ In the estimates above we used the best-fit values for the lepton mixing parameters and neutrino masses and marginalised over the two possible solutions for the Yukawas (see App. C) and the Majorana phase γ. The lower (upper) bounds ζ stemming from τ → eγ and τ → µγ are weaker, of the order of 10 −5 (10 5 ). The stronger limits from Eq. (24) apply unless there are fine-tuned cancellations among the η ± and η ± contributions to µ → eγ.

α → β γ γ
This process receives in general three independent contributions, i.e. from (i) γ-penguin, (ii) Zpenguin and (iii) box-type diagrams. We follow the notation of Ref. [19]. The γ-penguin amplitude for the transition α (p) → β (k 1 ) γ (k 2 ) γ (k 3 ) is described by where q = p − k 1 and the form factors A L/R 1,2 are reported in Eq. (20). The leading order contribution from the Z-penguin is proportional to the square of the charged lepton masses and thus negligible compared to the γ-penguin contribution. There are also box-type diagrams whose contributions are given by where for the decay − α → − β − γ + γ the form factor B reads and β it is given by where all the external momenta and masses have been neglected. The loop functions h 1 (x) and h 2 (x, y) are given in Eq. (88) in App. D. One can express the trilepton branching ratios in terms of the Wilson coefficients. Here we follow Ref. [24]. The branching ratio for α → β β β reads γ , as there are only box-contributions, we get [24] BR In the dipole dominance approximation, we can express µ → 3e in terms of µ → eγ as: We have checked in the numerical analysis that the above estimate is fulfilled to great precision, confirming that the box contributions are not relevant. They involve two extra Yukawa couplings and for Yukawa couplings smaller than one the box contributions are suppressed with respect to the dipole A L,R 2 and the monopole A L,R 1 contributions given in Eq. (20). The current experimental limits are BR(µ → 3e) < 1.0 × 10 −12 [25], with an expected future sensitivity of ∼ 10 −16 [26]. The rest of the trilepton decays involve taus, and their upper limits on their branching ratios are in all cases O(1) × 10 −8 [27], with future expected sensitivities of

µ − e conversion in nuclei
The conversion of a muon to an electron in a nucleus also imposes stringent constraints on the model parameter space. It is dominated by coherent conversions in which initial and final states of the nucleus N are the same. In this case the matrix elements of the axial-vector N |q γ α γ 5 q|N , pseudoscalar N |q γ 5 q|N , and tensor quark currents N |q σ αβ q|N vanish identically [28]. Similar to the purely leptonic decay of τ and µ leptons, the Z-penguin contribution to µ − e conversion in nuclei is proportional to the square of the charged lepton masses and thus negligible to the γ-penguin contribution. Moreover, the SM Higgs contribution is suppressed by the small firstgeneration Yukawa couplings. Thus µ − e conversion is dominated by photon exchange and the relevant terms in the effective Lagrangian contributing to µ − e conversion can be parameterised as [28] The long-range interaction mediating the process is given by the electromagnetic dipole transitions, whose effective couplings A L/R 2 were already introduced in Eq. (20), taking into account the appropriate flavour indices in the Yukawa couplings. The short-range interaction through the γ-penguin diagrams generate the vector current operator with where Q q is the electric charge of the quark q in units of e > 0 and A L 1 is the electromagnetic form factor given in Eq. (20) for the µ − e transition. A right-handed leptonic vector-current is not induced at one-loop because all new particles exclusively couple to the left-handed lepton doublets L. Accordingly, the µ − e conversion rate is where the effective vector couplingsg for the proton and the neutron arẽ while the numerical values of the overlap integrals D and V (p,n) and the total capture rate are reported in Tab. 5 in App. D. Notice that the neutron contribution is in our case approximately zero, as we neglected the Z-penguin contribution.
We can express ω conv in terms of BR(µ → e γ) using that where r g/f parameterizes the difference due to the different loop functions with 1 r g/f 1.5 for x 1 (remember m ψ is always smaller than m η ( )± , since ψ is the DM candidate), see Eqs. (20) and (87). In this case we can derive the allowed ranges for for Al (Au) {Ti}. We have checked that the numerical results for the µ − e conversion ratio are in very good agreement with the estimates in Eq. (38). The currently best experimental limits are CR conv (Au) < 7·10 −13 at 90 % CL by the SINDRUM II Collaboration [29] and CR conv (Ti) < 4.3 · 10 −12 at 90 % CL [29]. Future experiments are expected to improve the sensitivity by several orders of magnitude: PRISM/PRIME [30,31]

Lepton dipole moments
Electric dipole moments for the leptons occur in the model only at the two-loop level. 5 However, leptonic magnetic dipole moments occur, similarly to µ → eγ transitions, at one-loop level via γpenguin diagrams. They receive two independent contributions from the charged scalars running in the loop. They are given by (20). Notice the factor of two on the right side of Eq. (39), coming from the definition of ∆a [38,39].
In the case of the muon magnetic dipole moment, the discrepancy with the SM has positive sign, and therefore the model cannot explain it, as it gives a negative contribution, see Eq. (20) together with the loop function f (x) > 0. In any case, the predictions for |∆a µ | are always very small, |∆a µ | 10 −12 , as LFV limits are much stronger. The magnetic moments of the electron and the tau are subject to very weak limits.

Dark matter relic abundance
We will study the case where the DM candidate is the Dirac fermion ψ. There are several possible production channels in the early Universe. In the next subsection we argue that DM annihilation to SM leptons ψψ → ¯ , νν is too constrained by LFV processes. Coannihilation with components of the new scalars allows to obtain the correct relic abundance. We show in Fig. 2a the different channels for η ± , η 0 . Similar diagrams exist for the other new scalars.

Dark matter annihilations
The dominant DM annihilation channels are into a pair of charged leptons or a pair of neutrinos, shown in Fig. 2a. In the non-relativistic limit the s-wave annihilation cross sections to neutrinos and charged leptons are given by The minus sign between the different contributions originates from the presence of t (u) channels for the cases of η (η ) charged or neutral mediators. The factor of 2 smaller cross section for neutrinos originates from the fact that they are Majorana fermions. For identical neutrinos in the final state δ αβ = 1 and there is another factor of 1/2. We can conservatively estimate these cross sections in the limit of large and small |ζ| using the obtained limits on ζ in Eq. (24). They only differ whether primed or unprimed Yukawa couplings are present in the limit of equal masses of ψ and the new scalars, which we denote by m. Larger scalar masses only further suppress the annihilation cross section and thus we obtain the following conservative limit This allows us to use the limit on ζ from Eq. (24). A comparison with the typical freeze-out annihilation cross section, vσ th 2.2×10 −26 cm 3 /s, yields that the DM annihilation cross section is always too small to account for the observed relic density. Using the most stringent limit from µ → eγ on the parameter ζ in Eq. (24), we obtain for any ordering, which is valid unless cancellations occur.  An interesting way to break the proportionality of BR( α → β γ) and DM annihilation is to have the relic abundance set by vσ(ψψ → νν) . This can be achieved if the charged scalars are much more massive than at least one neutral scalar (η 0 in this model), giving therefore suppressed contributions to all α → β γ processes, and also to vσ(ψψ → − + ) . In this scenario the DM ψ is in the MeV range to obtain the correct DM relic density, with slightly heavier neutral scalar η 0 . However gauge invariance relates the interactions of SM neutrinos and charged leptons, and therefore this scenario requires some tuning of the parameters to circumvent experimental constraints (Z-boson decays, T parameter, LFV, etc.) and we thus do not consider it any further. Examples of similar scenarios have been studied in Refs. [40][41][42][43].

Dark matter coannihilations
As DM annihilation into charged fermions and neutrinos is strongly constrained by LFV processes, coannihilation processes may become important. The explanation of the correct DM relic abundance requires a small mass splitting between the DM candidate ψ and the scalars η 0 and η ( )± . The relative contribution of (i) annihilations of the DM particle with the coannihilation partner into a lepton and a SM gauge boson or Higgs boson in the final state (see Fig. 2b), and (ii) annihilations of the coannihilation partner via gauge interactions (γ/Z/W ) into SM particles, direct annihilations to Higgs bosons or DM-mediated t-channel annihilations into SM leptons (shown in Fig. 2c), depends on the size of the Yukawa couplings and the mass splitting. The coannihilation channels dominate the abundance, because the corresponding cross sections only depend on the square of one of the Yukawa couplings y α Φ , y α Φ compared to the annihilation cross section, see Eq. (40), which involves four Yukawa couplings. In our numerical scan we use micrOMEGAs 4.3.5 [44] to calculate the DM relic abundance and thus take all relevant (co)annihilation channels into account. See the seminal work [45] by Griest et al. for an analytic discussion of coannihilation.

Dark matter direct detection
One of the nice features of this model is that the fermionic DM is Dirac and not Majorana, unlike in the original ScM. Therefore direct detection, absent at tree level, has sizeable long-range contributions at one loop via photon exchange. 6 It can be parameterised by the magnetic (and electric) dipole interactions, namely In this model the electric dipole moment d ψ vanishes at one-loop level, because the DM ψ only couples to the left-handed lepton doublets and not to the right-handed charged leptons simultaneously [46]. The magnetic dipole moment is given by The loop function f DD (x, y, z) is defined in Eq. (89) in App. D. We verified our result against the well-known expressions for the magnetic moment for a Yukawa interaction reported in Ref. [47]. Similar results are given in Refs. [9,46,[48][49][50][51][52]. In Sec. 4 we show how results from the latest Xenon experiments XENON1T [53], PandaX-II [54], and LUX [55] constrain the parameter space using the LikeDM package [56].

Electroweak precision tests
The dominant contribution from new physics to electroweak radiative processes is generally expected to affect the gauge boson self-energies which are parameterised by the oblique parameters S, T , and U [57,58]. The limits on the oblique parameters are obtained from a global fit to electroweak precision data. The Gfitter collaboration finds the values: S = 0.05±0.11, T = 0.09±0. 13 and U = 0.01 ± 0.11, with correlation coefficients ρ ST = 0.90, ρ SU = −0.59 and ρ TU = −0.83 [59].
The strongest constraints on the model parameter space are set by the T parameter, which is sensitive to the mass splitting between the neutral and charged scalar components of the two inert doublets. We use the expressions for the oblique parameters in Refs. [60][61][62]. Details are reported in App. E.

Production and decay of the new scalars at colliders
Searches of neutral and charged scalar particles at colliders set constraints on the scalar mass spectrum of the model. In fact, from the precise measurement of the W and Z decay widths at LEP-II, the following kinematical bounds can be derived: At the LHC the production of these states proceeds mainly via neutral and charged current Drell-Yan processes. Other production channels are via an off-shell Z/W boson. A sub-leading contribution is given by Higgs mediated gluon fusion, provided the relevant couplings in the scalar potential Eq. (3) are sizeable [4,63].
In the case one of the charged scalars is the next-to-lightest particle in the dark sector, the expected signature at the LHC consists in the pair production of η ( )± followed by the prompt decay η ( )± → ψ ± α (α = e, µ, τ ). 7 The DM ψ escapes the detector and is revealed as missing transverse energy. The decay branching ratios of η ( ) into the different leptons only depend on the neutrino Yukawa couplings, namely Using the estimates for the Yukawa couplings y α Φ and y α Φ reported in Sec. 2.2 we expect for neutrino masses with NO that both charged scalars η ± and η ± have very similar branching ratios with the one to e ± being suppressed by two powers of the reactor mixing angle θ 13 with respect to those to µ ± and τ ± . Since θ 23 ≈ π/4, the branching ratios to the two flavours µ ± and τ ± are expected to be very similar for both NO and IO, see Eqs. (15) and (17), respectively. Moreover, for neutrino masses with IO the branching ratios of both charged scalars η ± and η ± to µ ± and τ ± are expected to be very similar, whereas the ones to e ± are expected to be different, but of similar size. We note that for IO BR(η ± → ψ e ± ) ≈ 2 BR(η ± → ψ µ ± (τ ± )) and BR(η ± → ψ e ± ) ≈ 2 BR(η ± → ψ µ ± (τ ± )). A measurement of at least one of the branching ratios may allow to extract information on the neutrino mass ordering and the Majorana phase γ, while there is only a weak dependence on the Dirac phase δ.
In the coannihilation region corresponding to m η ( ) 0 > m η ( )± m ψ , the decay width of the charged scalar is kinematically suppressed. As shown in the next section, in this case the lightest charged scalar may be long-lived, leaving ionising tracks in the detector [4,66].

Higgs and Z decays
If the scalars are sufficiently light, the Higgs and Z-boson can decay into them at tree level. In specific cases some of the limits on the neutral scalar masses from Z decays can be evaded by tuning the mixing angle θ. For instance, the Z-boson decay rate into the lightest neutral scalar, Γ(Z → η 0 η * 0 ), is proportional to cos 2 (2θ) and therefore vanishes in the case of maximal mixing, θ = π/4. 8 In this case the mass of the lightest neutral scalar, m η 0 , can be smaller than m Z /2. For the Higgs the decay to the lightest neutral scalar can be suppressed for sufficiently small quartic couplings and/or a suitable choice of the mixing angle θ.
There can also be Higgs and Z-boson decays at one-loop level. The new heavy charged scalars couple to the SM Higgs boson and thus modify the decay of the SM Higgs boson to two photons. The relative change of the SM Higgs partial decay width to two photons compared to the SM prediction can be parameterised as [67][68][69] where λ HΦ , λ HΦ are the couplings of the charged scalars η + , η + to the Higgs field, see Eq. (3).
As the new particles are not coloured and thus the Higgs production cross section is unchanged the signal strength is simply given by µ γγ = R γγ in this model. ATLAS observes µ γγ = 1.14 +0.27 −0.25 [70], and CMS µ γγ = 1.11 +0.25 −0.23 [71] which can be interpreted as a constraint on the new charged scalars. The combined measurement is µ γγ = 1.14 +0. 19 −0.18 [72]. If the charged scalar masses are light, deviations in the h → γγ channel are generically expected, but their size crucially depends on parameters in the scalar potential. This constraint can be easily fulfilled in the GScM.
In principle there can be new invisible decay channels of the Higgs and Z-boson to DM as well as of the Higgs to neutrinos. In general, the scalars are constrained to be heavier than ∼ 100 GeV due to a combination of collider searches, the limits from the invisible decay width of the Z-boson, and EWPT. Consequently, also DM cannot be light in the case of coannihilations and thus the Higgs and the Z-boson cannot decay into DM, which would in principle occur at one-loop level, see Ref. [46].
Other possible processes are LFV (and lepton flavour universality violating) Higgs and Z-boson decays, like h → τ µ and Z → τ µ. These, however, are very suppressed by a loop factor and due to experimental constraints from other LFV processes (like τ → µγ, etc.). They are therefore well beyond the expected sensitivity of future experiments [73].

Numerical analysis
The Yukawa couplings of the model can be expressed in terms of the neutrino oscillation data, the Majorana phase γ and the arbitrary parameter ζ (which can be taken without loss of generality to be positive) as explained in App. C. Therefore, we scan on the rest of the parameters of the model (and ζ) as outlined in Tab. 2, as well as on the 3σ range of the neutrino parameters (mixing angles, phases and masses) using the results from a the global fit of the group NuFIT 3.1 (November 2017) [74,75] Table 2: Priors on the 12 free real parameters used in the scan. λ i includes the following 8 quartic couplings of the potential: λ Φ ( ) , λ HΦ ( ) , λ HΦ ( ) ,2 , λ ΦΦ , λ ΦΦ ,2 . The parameter ζ is defined in App. C. We scan in these parameters using logarithmic priors.   [74,75]. Here ∆m 2 3 = ∆m 2 31 > 0 for NO and ∆m 2 3 = ∆m 2 32 < 0 for IO. We scan in these using flat priors.
We impose several constraints directly in the scan: (i) Direct searches on singly-charged scalars from LEP II imply m η ( )± 100 GeV, with some dependence on the channel. We impose a lower bound of 100 GeV on the scan; (ii) we also impose the 3σ constraints from EWPT, see Sec. 3.7. These constraints restrict the splittings of the scalars, specially the T parameter; (iii) thus we also demand than that the neutral scalars are heavier than 100 GeV motivated by EWPT in order to avoid contributions to the Higgs or Z decay widths; (iv) we apply the stability conditions on the scalar couplings given in App. A; (v) we impose the limits from radiative α → β γ as hard cuts; (vi) we assume that the fermion ψ constitutes all of the DM of the Universe and thus we impose its relic abundance to lie within the 3σ range of the latest results from Planck [76], Ω DM h 2 = 0.1198 ± 0.0026. We have evaluated the relic density numerically using micrOMEGAs 4.3.5 [44]. All the observables for which we impose constraints in the numerical scan are provided in Tab. 4.
In the following subsections we show the results of a numerical scan with about 10 4 random points, using the input parameters in Tabs. 2     We show in Fig. 3 left panel the lightest and heaviest neutral scalar masses, m η 0 and m η 0 , versus the DM mass m ψ , in red and blue, respectively. The lightest neutral scalar mass is close to the DM mass. This is completely driven by the fact that coannihilations need to be efficient enough to obtain the correct relic abundance. In the right panel, we show the neutral scalar masses versus the charged scalar mass m η + . We observe that η 0 is always the heaviest state. Roughly in 50% (30%) of the points the lightest neutral (one of the charged scalars) is very degenerate (with a splitting smaller than 5%) with the DM mass and contributes to the coannihilations. In addition, there are significant regions of the parameter space (18% of the points) where both masses of lightest neutral and one of the charged scalars (η ± or η ± ) are degenerate with the DM mass. Only in around 1% of the points the three lightest scalars η 0 , η ± , η ± are very degenerate with the DM mass. There is no difference between NO and IO.

Scalar mass spectrum
In Fig. 4, where we show the normalised mass splitting δ(m x ) = (m x − m ψ )/m ψ of the charged scalars, η + versus η + , in blue, and of the neutral scalars, η 0 versus η 0 , in red. Notice that the red region of points is superimposed over part of the blue one. The normalised mass splitting needs to be below ∼ 0.5, and typically ∼ 0.05, for at least one of the scalar states η 0 , η ± and/or η ± coannihilations to be efficient. It is typically much larger for the heaviest neutral scalar η 0 and thus the red points are all below the diagonal in Fig. 4.
In Fig. 5, left panel, we plot the lifetime of the neutral scalars, η 0 (blue) and η 0 (red), versus the normalised mass splitting δ(m x ) = (m x − m ψ )/m ψ . One can observe how the lifetime of η 0 can be much larger than that of η 0 . Indeed, when the splitting with the DM mass is small, the only decays of η 0 are into charged leptons, and even those can be closed for very small splittings, and/or suppressed for small Yukawa couplings. In Fig. 5, right panel, we plot the lifetime of the charged scalar η + versus δ(m η + ), for different ranges of λ HΦΦ : 10 −8 λ HΦΦ 0.01 in red, 0.01 λ HΦΦ 0.5 in blue, and 0.5 λ HΦΦ 4π in green. The plot for the charged scalar η + is analogous to that of η + . We observe two effects: first, for very large mass splittings (δ(m η + ) 0.1 which corresponds to m η + − m ψ 80 GeV), the main decay channel is η + → W + η 0 , and the larger the quartic coupling λ HΦΦ , the larger the neutral scalars mixing cos θ, see Eq. (7), and the larger this decay; second, the larger the mass splitting with the DM mass, the smaller the lifetime. Indeed, the charged scalar can be long-lived at colliders scales (τ η + 10 −8 s, shown with a horizontal dotted black line) for mass splittings smaller than ∼ 80 GeV, when the decay channel η + → W + η 0 is closed. In that region, the decays η + → + α ψ mediated by the Yukawa couplings y α Φ dominate. Therefore, the larger the quartic coupling λ HΦΦ , the smaller the Yukawa couplings, and the larger the lifetime. This is the reason why there are mostly blue and green points (i.e., medium and large λ HΦΦ values) in this region.

Neutrino masses
On the left panel of Fig. 6 we show |y Φ | versus |y Φ | for fixed intervals of λ HΦΦ : 10 −8 λ HΦΦ 0.01 in red, 0.01 λ HΦΦ 0.5 in blue, and 0.5 λ HΦΦ 4π in green. The Yukawa couplings are inversely proportional to each other as expected from the neutrino mass scale. Also, the larger 0.5 ≤ λ HΦΦ ≤ 4π 0.01 ≤ λ HΦΦ ≤ 0.5 10 −8 ≤ λ HΦΦ ≤ 0.01 the quartic coupling λ HΦΦ the smaller are the Yukawa couplings. For λ HΦΦ 0.5 the product of the Yukawa couplings is constrained to 10 −9 |y Φ ||y Φ | 10 −7 as shown in the plot. This is a direct consequence of the dependence of the couplings in the expression for the neutrino masses, see Eqs. (12) and (14).
We show in Fig. 6 (right panel) the product of the absolute values of the neutrino Yukawa couplings |y Φ ||y Φ | versus the mass splitting of the neutral scalars m η 0 − m η 0 for different ranges of | sin 2θ|: 10 −5 | sin 2θ| 0.001 in red, 0.001 | sin 2θ| 0.05 in blue, and 0.05 | sin 2θ| 1 in green. One can see two effects. First, the larger the mixing | sin 2θ|, the smaller the Yukawa couplings. This is expected as |y Φ ||y Φ || sin 2θ| is proportional to the scale of neutrino masses. Second, neutrino masses are also proportional to the splitting of the neutral scalars, and for a given value or range of | sin 2θ|, the larger the product of the neutrino Yukawa couplings, the smaller this splitting can be.

Constraints from lepton flavour violation
In Fig. 7, we plot the branching ratio of µ → eγ versus |y Φ |, for the same ranges of the quartic coupling λ HΦΦ used in Fig. 4 (right panel). The different sets of points form V-shaped regions whose minima are larger the lower the λ HΦΦ value. For |y Φ | 10 −4 , the branching ratio scales proportionally to |y Φ | 4 , independently of λ HΦΦ . In this region the contribution from the η + is suppressed because |y Φ | |y Φ |. If, however, |y Φ | 10 −4 η + dominates the rate. The dependence on λ HΦΦ again sets the scale of |y Φ | and thus the branching ratio of µ → eγ, i.e., the larger the quartic coupling λ HΦΦ the smaller BR(µ → eγ). The minimum value of these LFV processes occurs for |y Φ | ∼ |y Φ | ∼ 10 −4.5 , when both charged scalar contributions are of similar order, such that the overall magnitude is suppressed.
In Fig. 8 we plot BR(τ → µγ) (left panel) and BR(τ → eγ) (right panel) over BR(µ → eγ), for NO (in red) and IO (in blue). The central values of these ratios agree with our analytical estimates of Eqs. (21) and (22)  magnitude in the branching ratios. In particular, we see that BR(τ → e γ) is suppressed compared to BR(τ → µ γ) for neutrino masses with NO while they are very similar for IO. Also the largest branching ratio is achieved for τ → µγ for NO, which can be larger than the one for IO. Therefore a measurement could in principle discriminate between the neutrino mass orderings.

Lepton flavour violation and dark matter direct detection
In Fig. 9 we plot the branching ratio of µ → eγ (left axis) and the µ − e conversion rate in Al (right axis) versus the DM magnetic moment µ ψ (provided in Eq. (45)) relevant for direct detection. The magnetic moment µ ψ is correlated with the LFV processes, because the structure of the loop diagrams is similar, with a charged scalar in the loop and the same Yukawa couplings. The points in red (blue), corresponding to larger (smaller) values of µ ψ , are excluded (allowed) by the combined constraint from Xenon-based DM direct detection experiments, which are implemented in the LikeDM package [56]. We can see the interesting complementarity between direct detection and lepton flavour violation in constraining the DM parameter space. This interplay is further discussed in the generic context for a fermionic singlet DM particle in Ref. [46].

U(1) DM as a gauge symmetry
The U(1) DM symmetry is anomaly-free, because the fermion ψ is vector-like and thus it can be straightforwardly gauged. In fact a similar model with a gauge symmetry has been discussed in Ref. [2]. Three scenarios can be envisaged: (i) U(1) DM is unbroken and the corresponding dark photon γ DM is massless, (ii) U(1) DM can be realised non-linearly and the dark photon γ DM obtains its mass from the Stueckelberg mechanism, or (iii) U(1) DM is spontaneously broken by a nonvanishing VEV of an additional scalar field ρ to a residual Z N symmetry which stabilises the DM candidate. In case (i) γ DM will contribute to extra radiation and lead to large self-interactions, 0.5 ≤ λ HΦΦ ≤ 4π 0.01 ≤ λ HΦΦ ≤ 0.5 10 −8 ≤ λ HΦΦ ≤ 0.01 Figure 7: Branching ratio of µ → eγ (left axis) and the µ − e conversion rate in Aluminium (right axis) versus the absolute value of the neutrino Yukawa couplings y Φ , for the same ranges of λ HΦΦ used in Fig. 4 left panel. see Ref. [77] for a discussion. If in case (ii) and (iii) the mass of the dark photon γ DM is smaller than that of the DM candidate, its relic abundance is set by DM annihilations into dark photons. Connections of DM phenomenology to neutrino and flavour physics are thus lost so that this case is not interesting to us. In addition, in case (iii) the new scalar field ρ will mix with the SM Higgs doublet H. Such mixing is experimentally constrained by invisible Higgs decays, if these are kinematically accessible, and by direct detection limits specially for sub-GeV scalar mediators, thus requiring ρ to be either much heavier than the electroweak scale or the mixing to be small. That in turn is in conflict with the fact that the mediator needs to decay before Big Bang Nucleosynthesis, which basically rules them out [78].
Another effect of gauging U(1) DM is the kinetic mixing with U(1) Y , that is the term B µν DM B µν   45)) relevant for direct detection. We plot in blue (red) the points that are allowed (excluded) by direct detection limits.
with the field strength tensors B µν (B µν DM ) of hypercharge (U(1) DM ). Even if this is tuned to vanish at some scale, it will arise at one-loop level, since the scalars Φ and Φ are charged under both U(1) Y and U(1) DM . This effect can be easily estimated. In fact, assuming Λ UV > m Φ , m Φ > v H and a vanishing kinetic mixing at a certain high scale Λ UV , (Λ UV ) = 0, the renormalisation group running can induce a sizeable kinetic mixing at lower scales. Above the mass scales of the two scalar doublets Φ and Φ , the opposite values of the hypercharge of the two scalars will lead to an exact cancellation of , but as soon as one of the scalars decouples, the kinetic mixing term is induced via renormalisation group running, giving the approximate lower bound with the dark gauge coupling α DM = g 2 DM /4π. The annihilation into light dark photons is given by vσ ann πα 2 DM /m 2 ψ , which implies that in order to reproduce the relic abundance α DM 10 −4 (m ψ /GeV). Using the experimental value of α Y (m Z ) [20], and taking the logarithm to be O(1), we can estimate a lower bound on the kinetic mixing: | | 10 −4 (m ψ /GeV). As in the case of scalar mediators, there are very strong upper limits from direct detection, specially if the dark photons are below the GeV, with only much smaller mixings still allowed, see the recent analysis by the PandaX-II collaboration [78]. Also, the mediators should decay before Big Bang Nucleosynthesis. Therefore, certain amount of fine-tuning is required in this case.

Replacing U(1) DM with Z N
Instead of a global U(1) DM symmetry we can consider a Z N symmetry, being a subgroup of U(1) DM , i.e. we regard the charge of Φ, Φ and ψ as given modulo N .
For N > 4 the model effectively possesses a global U(1) DM symmetry, as long as we only consider renormalizable terms in the Lagrangian. Conversely, for N = 2, 3, and 4, additional terms arise at the renormalizable level. For N = 2, i.e. the smallest possible symmetry, we can identify ψ ↔ ψ c and Φ ↔Φ. The model thus contains one Majorana fermion ψ and one scalar doublet Φ which are the only particles odd under the Z 2 symmetry. The Lagrangian for the Majorana fermion Ψ = ψ + ψ c reads The scalar potential V for Φ and the SM Higgs doublet H becomes The model with a Z 2 symmetry, a Majorana fermion and one additional Higgs doublet has the same symmetries and types of particles as the ScM, which has been extensively discussed in the literature [3]. We comment on similarities and differences in phenomenology between the latter and our model with a global U(1) DM symmetry in the next section. For N = 3 and N = 4 the Lagrangian L ψ remains the same as in Eq. (1), but additional quartic terms appear in the scalar potential, see also Ref. [79]. For N = 3 there are two new quartic couplings In the case in which one of the neutral scalars is the lightest particle with non-trivial Z 3 charge, these terms give rise to DM semi-annihilations [80][81][82][83]. In principle, these new couplings are directly testable at colliders. Furthermore, for N = 4 the following term can be added to the scalar potential Again, if one of the neutral scalars is the lightest particle with non-trivial Z 4 charge, this term gives rise to DM self-interactions.

Comparison with the Scotogenic Model
As already mentioned in the preceding section, if the global U(1) DM symmetry is replaced by a Z 2 symmetry, our model coincides in symmetries and types of particles with the ScM [3]. In the following, we highlight similarities and differences between the GScM discussed here and the ScM. In order to generate at least two neutrino masses, in addition to the SM Higgs doublet, two Majorana fermions ψ 1,2 (with masses m ψ 1,2 ) and one additional scalar doublet are needed in the ScM, whereas in our model one Dirac fermion and two scalar doublets are needed. We have thus one more charged and one more neutral complex scalar in our model compared to the ScM. In the latter model the scalar mass spectrum, derived from the potential in Eq. (50), reads with φ R 0 and φ I 0 being the real and imaginary components of φ 0 , φ 0 ≡ (φ R 0 + i φ I 0 )/ √ 2. The physical scalars φ R 0 and φ I 0 acquire a mass splitting proportional to the quartic coupling λ HΦ, 3 . In contrast, in our model the mass spectrum, given in Eq. (8), clearly shows that real and imaginary parts of the neutral scalars have the same mass and combine to form complex neutral scalars, denoted η 0 and η 0 .
In the ScM the neutrino masses are generated by the real and imaginary parts of the neutral scalar running in the loop. The mass matrix is given by where the loop function is defined in Eq. (13). We can see how the difference in mass between the two complex neutral scalars η 0 and η 0 in our model, that appears in the neutrino masses, see Eq. (13), is traded for the difference in mass between the neutral scalars (real and imaginary parts) φ R 0 and φ I 0 in the ScM. As is well-known, in the ScM lepton number is broken by the simultaneous presence of the Yukawa couplings, the masses of the Majorana fermions ψ 1,2 and the quartic coupling λ HΦ,3 of the scalar potential in Eq. (50). Thus neutrino masses crucially depend on all three of them. While the dependence on the first two ones is obvious from Eq. (56), the one on λ HΦ,3 is easiest revealed in the limit m 2 )/2 where the expression for the neutrino mass matrix takes the form This is similar to what happens in the GScM, where the simultaneous presence of both Yukawa couplings y α Φ and y α Φ , the fermion mass m ψ and the quartic coupling λ HΦΦ is required in order to break lepton number. Consequently, neutrino masses are proportional to all these quantities, as can be read off from Eq. (12) together with Eqs. (7) and (8).
In the original ScM LFV processes have been studied in detail in Refs. [3,5]. It turns out that the individual penguin diagram contributions of a charged scalar and a fermion to LFV processes in the original ScM are the same as the ones in the GScM, see Sec. 3. However, the number of charged scalars and fermions differs in the two models. We therefore obtain a different number of contributions to LFV processes in the two models. Moreover, there are new box diagrams for tri-lepton decays in the GScM.
The DM phenomenology is different in the two models. For scalar DM the main channels for direct detection in the ScM are the tree-level mediated processes by the Z-boson and the Higgs. For λ HΦ,3 = 0 the scalar (φ R 0 ) and pseudoscalar (φ I 0 ) have different masses and DM scattering off nuclei is an inelastic process, with the Z-boson exchange typically dominating. This imposes a lower bound on λ HΦ,3 in order to kinematically forbid such scattering. In the GScM, scalar DM (η 0 ) also naturally has a very large direct detection rate mediated by the Z-boson, unless some cancellation occurs in Z interactions (like for maximal mixing θ = π/4). Moreover, there is an elastic contribution via the Higgs portal.
For fermionic DM in the ScM, DM-nucleon scattering occurs at loop-level [9] (similar as in Sec. 3.6 for the GScM) via penguin diagrams. If the mass splitting between ψ 1 and ψ 2 is sufficiently small, there is a transitional magnetic moment interaction with the charged leptons in the loop. This leads to inelastic DM-nucleon scattering. As discussed in Sec. 3.6 the dominant DM-nucleon scattering also occurs via a magnetic moment interaction with the charged leptons in the loop.

Conclusions
We have studied a model in which neutrino masses are absent at tree level, while they are radiatively generated at one-loop level with a fermionic WIMP DM particle running in it. The model considered is a generalised version of the Scotogenic Model (GScM) in which the symmetry is not discrete, but global. The minimal model has one massless neutrino and we have found that both neutrino mass orderings (NO and IO) can be accommodated. The flavour structure of the Yukawa couplings can be expressed in terms of the neutrino oscillation parameters and the Majorana phase γ. Thus the model possesses a definite flavour structure and is very predictive. The relevant flavour structure for neutrino masses is different from that involved in LFV, which is an interesting feature of the model, in contrast to the original Scotogenic Model (ScM). We have obtained interesting correlations among the ratios of different lepton flavour violating processes, which may allow to test the model and to discriminate between the two neutrino mass orderings.
In this work we focused on fermionic DM, given the fact that scalar DM would require fine-tuned scalar mixing. The main WIMP annihilation channels are into charged leptons and neutrinos. As they are mediated by the same Yukawa couplings involved in LFV processes, they are too suppressed to explain the observed DM relic density and thus coannihilations with one of the new scalar particles are needed. 9 Roughly in half of the parameter space, the next-to-lightest particle is the lightest neutral scalar (η 0 ), and in the other half it is the lightest charged scalar (which can be either η + or η + ).
LFV processes and DM direct detection give complementary information on the parameter space of the model. Future µ − e conversion and µ → 3e experiments will provide the most sensitive probe of the remaining allowed parameter space, but also direct detection experiments will further test a complementary region of the allowed parameter space. Another interesting signature of the model is the direct production of two charged scalars via the Drell-Yan process at colliders and their subsequent decay to a charged lepton and DM. For NO the dominant channels are into taus and muons, while for IO the decay into electrons is of similar magnitude.
In comparison to the original ScM, the generalised version discussed here has more degrees of freedom in the scalar sector (two doublets versus one in the ScM), but it needs only one Dirac fermion (vector-like), unlike the ScM which needs at least two Majorana fermions. The flavour structure in the GScM is determined up to an overall normalisation, while the ScM with three fermionic singlets provides more freedom. They both are simple explanations for neutrino masses and DM with a rich and testable phenomenology.
Only the last condition involves λ HΦΦ and bounds the latter from above.
V 4 is also minimised for non-extremal values of ρ 12 and ρ 13 . This, however, does not imply conditions different from those already shown above, but only leads to an equality involving λ HΦΦ , λ HΦ,2 and λ HΦ ,2 which needs to be fulfilled in addition. We have checked that the presented conditions can also be applied to the special directions in which one or two of h 1,2,3 vanish.

B Neutrino masses and lepton mixing
We can diagonalize the mass matrix as where ν L are the left-handed neutrino flavour states, D ν is an 3 × 3 diagonal matrix with positive real eigenvalues (in our model with m 1 = 0 for NO, and m 3 = 0 for IO) and U is the PMNS mixing matrix, which relates the neutrino mass states ν i (i = 1, 2, 3) with definite masses m i and the neutrino flavour fields ν α (α = e, µ, τ ): The standard parametrisation for U for one massless neutrino is where c ij ≡ cos θ ij and s ij ≡ sin θ ij (θ 12 , θ 13 , and θ 23 being the three lepton mixing angles). γ is the Majorana and δ is the Dirac phase. Since the lightest neutrino is massless in the GScM, there is only one physical Majorana phase.

C Parametrisation of Yukawa couplings
We want to express the Yukawa couplings in terms of neutrino masses and lepton mixing parameters. We follow the discussion in Ref. [86]. On the one hand the rank-two neutrino mass matrix can be expressed in terms of the two non-vanishing mass eigenvalues and two columns of the PMNS matrix in the flavour basis The vectors v i ≡ √ m i u * i are linearly independent and span a two-dimensional vector space. On the other hand the calculation of the neutrino mass matrix results in the following form with (see Eqs. (12) and (13)) If sin 2θ < 0 the square roots will yield a complex number. The two linearly independent vectors x i can be expressed in terms of the vectors v i where (a ij ) forms an invertible 2 × 2 matrix, i.e. det a = a 11 a 22 − a 12 a 21 = 0. Using the two different parametrizations of the neutrino mass matrix, we can find possible solutions of a ij : In particular the matrix elements are non-zero, a ij = 0. There are two sets of solutions. Writing the complex matrix element a 11 = re iα in terms of two real parameters r, α we obtain is trivially satisfied for the two solutions. Thus the two vectors can be uniquely written as For Normal Ordering (NO) we obtain while for Inverted Ordering (IO), we have that Without loss of generality we choose ζ to be real. Any phase of ζ can be absorbed with phase redefinitions of the lepton doublets L and the fermion singlet ψ. In Figs. 10 and 11 we show the results for y α Φ (left) and y α Φ (right) for neutrino masses with NO and IO, respectively, separated according to real (in red) and imaginary (in blue) parts, as obtained in the numerical scans. We show different flavours: the upper panel shows the flavour τ versus µ, while the lower shows µ versus e. We clearly see that y α Φ and y α Φ of different flavour α are correlated. These correlations can be understood analytically, see Sec. 2.2.

D Loop functions
In this appendix we write the loop functions and other ingredients appearing in LFV processes, h → γγ and DM direct detection.
The following ones enter in the box diagrams of trilepton decays: The numerical values of the overlap integrals D and V (p,n) for µ − e conversions are reported in Tab. 5 for three different nuclei. We also report in the table the total capture rate for each nucleus.

E Oblique parameters
The two inert doublets contribute to the EWPT at one-loop level. The contribution in our model to the T parameter is given by [60,61] T = 1 16π 2 α em v 2 H 2 s 2 θ F(m 2 η + , m 2 η 0 ) + 2 c 2 θ F(m 2 η + , m 2 η 0 ) + 2 c 2 θ F(m 2 η + , m 2 η 0 ) + 2 s 2 θ F(m 2 η + , m 2 η 0 ) , where the loop function is defined as F(x 2 , y 2 ) = x 2 + y 2 2 − x 2 y 2 x 2 − y 2 ln The loop function is symmetric in x and y. It vanishes in the custodial symmetry limit, x → y, and diverges for x/y going to 0 or infinity. Extending the results of Ref. [61], the S parameter reads in our model with the definitions in d space-time dimensions, where γ E 0.577 is the Euler-Mascheroni constant. Note that B 22 and B 0 are symmetric in their last two arguments. We use the compact analytic expressions given in App. B of Ref. [62]. We confirmed that the expressions agree with the ones in the inert doublet model [88] when taking the appropriate limit.