Dark Matter candidates in a Type-II radiative neutrino mass model

We explore the connection between Dark Matter and neutrinos in a model inspired by radiative Type-II seessaw and scotogenic scenarios. In our model, we introduce new electroweakly charged states (scalars and a vector-like fermion) and impose a discrete $\mathbb{Z}_2$ symmetry. Neutrino masses are generated at the loop level and the lightest $\mathbb{Z}_2$-odd neutral particle is stable and it can play the role of a Dark Matter candidate. We perform a numerical analysis of the model showing that neutrino masses and flavour structure can be reproduced in addition to the correct dark matter density, with viable DM masses from 700 GeV to 30 TeV. We explore direct and indirect detection signatures and show interesting detection prospects by CTA, Darwin and KM3Net and highlight the complementarity between these observables.


Introduction
The presence of Dark Matter (DM) and its role in the formation of large scale structures of the universe [1] altogether with the observation of neutrino oscillations [2][3][4][5][6][7] are some of the strongest indications that the Standard Model (SM) of particle physics lacks essential ingredients. Indeed, the SM provides a successful description of the microscopic interactions with a remarkable accuracy. However, neutrinos, which are precisely massless within the SM, should possess non-vanishing masses in order to explain successfully the observed oscillation patterns. In addition, the SM do not include a viable dark matter candidate among its large particle content. Both topics have became the ground base motivation for a large variety of the Beyond-the-Standard-Model (BSM) constructions. As stated before, neutrino oscillations require massive neutrinos and therefore the SM must be extended in order to accommodate non-vanishing masses. One of most famous and minimal mechanisms accounting for non-vanishing neutrino masses is known as the "see-saw" mechanism [8][9][10][11][12][13]. In seesaw-based models, extra fields such as SU (2) L fermion singlets, triplets, scalar triplet, with or without hypercharge, are added to the particle content of the SM, allowing to generate neutrino Majorana mass terms. On the effective point of view, the lowest dimensional operator responsible for Majorana mass terms is the Weinberg operator which is a dimension-5 operator constructed with SM fields that preserves every SM gauge symmetries but violate lepton number conservation by 2 units: where L i are the lepton doublets, H is the SM higgs doublet, and Λ is a mass scale. The Majorana mass term arises after the SU (2) L × U (1) Y spontaneous symmetry breaking (SSB): Within the BSM models landscape, the Weinberg operator is typically generated by integrating out of the spectrum heavy mediators. Models reproducing light neutrino Majorana masses at tree-level by adding a minimal number of extra fields can be generally sorted in three main categories. These categories are commonly referred to as Type-I, Type-II, and Type-III seesaws (see [14] for an overview). Among them, diagrams responsible for the Weinberg operator in T-I and T-III seesaw models are similar. Within the T-I seesaw framework, a right-handed neutrino (Majorana fermion SU (2) L singlet) is introduced instead of the Majorana fermion SU (2) L triplet used in the T-III. T-II seesaw is different as a SU (2) L scalar triplet with hypercharge is added in this case. This field carries 2 units of lepton number, which is broken spontaneously once this scalar acquires a non-vanishing vacuum-expectation-value (vev), generating subsequently neutrino mass terms proportional to this vev. Beyond this framework, non-vanishing neutrino masses can also be generated at the loop level, providing an additional argument to justify the large hierarchy between the electroweak scale and the light neutrino masses. For these cases, an interplay between neutrino masses and DM presents interesting features. First of all, tree-level processes could become negligible or forbidden by invoking an additional discrete or gauge symmetry [15,16]. For instance the simplest realization is typically obtained by considering a Z 2 discrete symmetry. This symmetry also protects one-loop processes, making them the only contributors to the neutrino masses. As a side effect of this symmetry, particles running inside the loop are protected, among these, the lightest state is stable and could be a viable DM candidate. There are several scenarios considering DM candidates in this framework (see for instance [17]). However, we highlight the models known as "scotogenic" , in which, the Weinberg operator is generated at one-loop by using as basis the Type-I [18] or Type-III [19] seesaw BSM fields which are, in this case, charged under the new symmetry. Constructions combining both models present interesting DM and neutrino phenomenology (for instance [20][21][22][23]). Nevertheless, models using the BSM fields of the Type-II seesaw are less common due to the ambiguity in their definition [24][25][26]. In addition, some of the phenomenological aspects are not explored. In this paper we investigate the phenomenology of a model accounting both for the observed dark matter abundance in our universe in addition to light neutrino masses, based on a Type-II seesaw mechanism.
Our manuscript is constructed as follows: in Section 2, we provide a description of various aspects of the model, such as the particle content and corresponding charge assignments. In Section 3, we explore the neutrino physics, DM phenomenology and lepton flavour violating observables implied by the model. The results of an exhaustive scan of the parameter space is presented in Section 4 which is subsequently analyzed and commented. We consider further constraints and discuss additional relevant observables in Section 5 before concluding in Section 6.

The Model
In addition to the SM particle content, we introduce four new fields charged under a discrete Z 2 symmetry with non-trivial charge assignment under the SM gauge symmetry group: a hyperchargeless Y = 0 real SU (2) L triplet scalar Ω, a hypercharged Y = 1 complex SU (2) L triplet scalar ∆ and two fermionic SU (2) L doublets f L,R with opposites chiralities but identical quantum numbers. The SM particle content is uncharged under the discrete Z 2 symmetry or equivalently Z 2 -even. The field content of this model is summarized in Table 1. The Z 2 discrete symmetry will play a two-folded role here. On the first hand, as in the original scotogenic model [18], neutrino masses are generated at the loop-level, by receiving contributions of fields charged under a discrete Z 2 symmetry. On the other hand, the neutral fields charged under the Z 2 will render a suitable stable Weakly-Interacting-Massive-Particle (WIMP) DM candidate. There are two meaningful differences arising from the original proposal [18]. The inclusion of one triplet field with Y = 1 and also charged under a Z 2 symmetry comes with doubly charged scalars which must decay in pairs of particles that are also charged under the Z 2 , thus forbidding the well-known smoking gun prospect of ∆ ++ −→ + α + β existing in models extending the TII-seesaw [27]. The inclusion of a second triplet also charged under the Z 2 , but with Y = 0 renders two features. Firstly, it gives rise to a loop made-up of charged particles contributing to neutrino masses. This comes in a similar way as it comes in Bilinear R-parity violation, where charged loops are fundamental in order to give rise to a full description of neutrino masses. Secondly, this second triplet also renders a richer scalar Z 2 -odd sector, with different features as other scotogenic models, for instance the model of Ref. [21]. The Lagrangian allowed by gauge invariance can be parametrized as where L i with i = 1, 2, 3 denote the i−th flavour SM lepton SU (2) L doublet. y i ∆,Ω are dimensionless free parameters carrying a flavour index i and m f is a free mass parameter. In this work we assume that the new terms are CP-conserving and therefore we consider the Yukawa couplings as real parameters. c denotes the Lorentz charge conjugation operation and V scalar is the scalar potential is explicited in Sec. 2.2. Notice that the absence of righthanded neutrinos implies the absence of couplings between lepton doublets L i and the SM Higgs. On the other hand, the new triplets are coupled, via Yukawa terms, to the neutrinos. Their components can be parametrized as where the +, − exponents on the fields represent their electric charges. As discussed further on, such Z 2 -odd charged and neutral scalar fields induce radiative contributions to neutrino masses. In the following we provide a more detailed description of the fermionic and scalar sector of the model.

Fermion sector
The new fermionic SU (2) L doublets can be parametrized as f T L,R = (f + L,R , f 0 L,R ). At treelevel since only Dirac mass terms are present in the new fermionic sector, the two doublets with opposite chiralities combine to form Dirac f + = f + L + f + R and f 0 = f 0 L + f 0 R which are the mass eigenstates with respective charges Q = 1 and Q = 0. In term of these fermionic mass eigenstates the Lagrangian at tree-level can simply be written as Notice here that the Q = 1 and Q = 0 fermionic states are degenerate at tree-level with masses m f + = m f 0 = m f . However, loop corrections give rise to a mass splitting ∆m f = m f + − m f 0 between these two components such that the mass terms in Eq. (2.3) can be written at the one-loop level as The relative mass splitting is typically small ∆m f /m f 1% for the relevant part of the parameter space, as shown further on. For instance ∆m f ∼ 10 −4 m f for m f = 1 TeV. The analytical expression for the mass splitting at the one-loop level ∆m f can be found in Appendix C.

Scalar Sector
The scalar sector of this model is composed of the Higgs field H, a Y = 1 triplet ∆ and a Y = 0 triplet Ω. The only Z 2 -even scalar is the Higgs doublet H and it is responsible of the electroweak spontaneous symmetry breaking. In the following, we explicit all the terms present in the scalar-sector Lagrangian as well as the various mass eigenstates and eigenvalues, obtained after diagonalization of the scalar mass matrix.

Kinetic terms
The kinetic terms for the scalar triplets can be written as where the covariant derivatives are defined as and using standard notations from the literature. We parametrize the neutral scalar fields as a function of real scalars where ∆ 0 R , Ω 0 R and ∆ 0 I denote the real and imaginary parts of ∆ 0 and Ω 0 .

Scalar Potential
In this section we explore the scalar potential and the mass eigenstates arising from this sector. First, a comment regarding the Higgs field. Since it is the only SU (2) and Z 2 -even field, it constitutes the only source of electroweak symmetry breaking, and thus the tadpole equations will not be trivially satisfied which implies a non-zero vev. The resulting neutral pseudoscalar mode obtained after expanding around the minimum will be the Goldstone mode responsible of the longitudinal modes of the Z-boson. In the same footing, the charged component of the Higgs field will constitute the longitudinal mode of the W -bosons. Thus, all the physical scalar fields within this model but the Higgs field are Z 2 -odd. The full scalar potential of the theory, including all renormalizable terms allowed by gauge invariance, can be parametrized as Tr Ω † Ω + λ Ω 16 Tr where s κ = ±1 is the sign of the dimensionless coupling κ chosen to be positive. All λ

Masses and mixings
The electroweak spontaneous symmetry breaking induces extra terms, proportional to the Higgs vev v h and induced by the κ-coupling, in the scalar mass-matrix which result in a mixing between scalars with identical quantum numbers, namely 2 CP even neutral scalars and 2 charged scalars with Q = 1. In the following, we explicit the mass eigenstates and eigenvalues for these scalars as well as the masses of the neutral CP-odd and Q = 2 scalars.
Neutral CP-even scalars: After electroweak symmetry breaking, the mass matrix for the CP-even neutral scalars in the basis (Ω 0 R , ∆ 0 R ) reads that can be diagonalized by the rotation where R 0 is a rotation matrix, and S 0 1 and S 0 2 are the mass eigenstates with corresponding masses where the − sign applies to m S 0 1 and the + sign to m S 0 2 . Further details on the mixing angle can be found in Appendix A.
Neutral CP-odd scalar: Only one CP-odd scalar is present in the spectrum, which is denoted byS 0 in the following. This particle corresponds to the imaginary part of ∆ 0 , S 0 ≡ ∆ 0 I . The mass of this scalar is given by

14)
Charged scalars with Q = 1: The mass matrix for charged scalars in the basis (Ω ± , ∆ ± ) reads that can be diagonalized via a rotation matrix connecting the gauge and mass basis by the transformation: where S ± 1 and S ± 2 are the mass eigenstates with masses where the − sign applies to m S ± 1 and the + sign to m S ±

2
. Further details on the mixing angle can be found in Appendix A.
Charged scalar with Q = 2 One charged scalar with Q = 2 is present in the spectrum. This particle corresponds to the state ∆ ++ and will be denoted in the following by S ±± for notation consistency, with a corresponding mass given by

Boundedness from below conditions for the scalar potential
In order to ensure stability in the scalar sector, the scalar potential in Eq. (2.9) must satisfy the so-called "boundedness from below" (BFB) conditions. Such conditions can be formalized as relations among the couplings in such a way that for a scalar potential with N complex scalar fields φ i with i = 1, .., N : where V min is the potential minimum (see [29][30][31] for further details). For a renormalizable theory, the scalar potential is a polynomial function of the fields. Therefore, the condition is translated to The shape of the potential is controlled by the couplings, this sets relations between various couplings in order to have a correct vacuum and to avoid the transitions to "unbounded" vaccua. The scalar potential in Eq. (2.9) can be written in terms of the mass eigenstates, in such a way that it has the functional form: where Φ i is an array containing the scalar fields of the model both neutral and charged. We perform a partial analysis in the gauge basis of the fields. Many of the terms within M 2 ij , ω ijk , and λ ijkl actually vanish due to gauge invariance and symmetry conditions. In our case, we focus on the quartic-terms as they govern the behaviour of the potential when |Φ i | → ∞. To obtain the BFB conditions, we compute the fourth derivatives for every combination of the fields to obtain the couplings: where ξ is the normalization factor: ξ(i, j, l, k) =12(δ ij δ ik δ jl δ kl ) + 2(δ ij + δ ik + δ il + δ jk + δ jl + δ kl ) needed to take into account the numerical factors from the derivatives. We classify the couplings depending of their parity under the sign change of any of the fields. Using this criterion, the couplings λ iiii and λ iijj (i = j) are even, and the rest of the couplings: λ iijk , λ iiij , and λ ijkl are odd for i = j = k = l. In general terms, all couplings are constrained by perturbatibility: 0 ≤ |λ ijkl | < √ 4π. However, the λ iiii must be always positive. The remaining couplings might be either positive or negative but restricted by the value of λ iiii .
In our case, many of the conditions are either trivial or redundant with other relations. Nevertheless, the analysis of the fourth derivative presents a minimal approach to the BFB conditions. Due to the complexity of the potential, there might be non-trivial correlations between fields going to infinity that can lead to an unbounded vaccua. The relations coming from the λ iiii are the following: (2.24) and the minimal copositivity criteria λ iijj + λ iiii λ jjjj ≥ 0 leads to: (2.25)

Dark matter candidates
As the Z 2 discrete symmetry is unbroken in this model, any interaction process must feature an even number of Z 2 -odd particle. As a result, the lightest neutral Z 2 -odd particle of this model is stable, electrically neutral and features weak interactions, implying that it is a viable WIMP-like dark matter candidate. In our model, there are four neutral particles: 3 scalars S 0 1 , S 0 2 ,S 0 and one fermionic state f 0 in the mass eigenstate basis. However, the mass of the state S 0 2 is larger than S 0 1 by definition of the mass eigenstates and eigenvalues in Eq. (2.13). In addition, one can show that the pseudoscalarS 0 mass is larger than the mass of S 0 1 . This statement is not very obvious by comparing respective expressions for the masses of these states. However, using the alternative parametrization described in Appendix B, Eq. (B.7) makes this statement more explicit. As a result, we are left with two viable dark matter candidates, the S 0 1 scalar and the neutral f 0 fermion.

Phenomenological implications: neutrinos, dark matter, and lepton flavour violation
In this section we provide details regarding the neutrino masses and dark matter phenomenology. We provide analytical expression of the neutrino masses in Sec. 3.1. In the following subsections, we discuss the dark matter relic abundance as well as dark matter direct and indirect detection. The complete numerical analysis of the parameter space is postponed to Sec. 4.

Neutrino Masses
As advanced, neutrino masses are generated in the scotogenic model at the loop level by fields that are odd under the Z 2 symmetry. Such one loop contributions are generated by self-energy diagrams involving loops of neutral fields f 0 , S 0 1,2 , as depicted in Fig. 1 in the flavour-eigenstates basis. In our particular case, additional charged-scalar loops exist [32]  generated by f ± , S ± 1,2 . Using the DR scheme, the light neutrino masses are finite and can be expressed as where the loop function is given by Therefore, neutrino masses are essentially controlled by the Yukawa couplings y i Ω , y i ∆ , the coupling κ, via the mixing angles, which plays the role of a lepton number violating parameter, in addition to masses of the scalar states and new fermions. In the following, for simplicity, we assume a normal ordering (NO) for the neutrino masses. In this case the latest determination of the neutrino-sector parameters are given by Ref. [33], whose best-fit values and ±1σ range are In general terms, the Yuwaka sector (y i Ω , y j ∆ ) can contain 6 complex phases. Among them, 5 phases can be absorbed by field redefinitions and the remaining phase can be related to δ. The impact of the CP-phase in our context is reduced due to the insensitivity of DM observables to neutrino flavours. For simplicity, as mentioned previously, we assume CP conservation therefore we set the usual Dirac phase δ to 0. CP-violation could give rise to a contribution to the electron dipole moment at the two-loop level which is constrained at the level of |d e |/e < 1.1 × 10 −29 cm [34]. Such effects have been studied in a similar context [35,36] and have been shown not to reduce significantly the available parameter space compatible with constraints from lepton-flavour violating observables. These effects are not expected to differ significantly in our setup.
In the following, we impose the neutrino-sector parameters of our model to reproduce the best fit values of Eq. (3.3) from Ref. [33]. In addition, we consider the lightest neutrino mass-eigenstate ν 1 as massless. Assuming this hierarchy reduces significantly the number of the free Yukawa couplings of the model to just one. Details regarding the parametrization used to reproduced parameters from Eq. (3.3) can be found in Appendix B.2.

Dark matter direct detection
Both fermionic and scalar dark matter candidates possess electroweak charges and significant interaction with light quarks therefore could potentially be observed by scattering on some detector material nuclei. In the following we evaluate the DM-nuclei scattering cross section for our potential fermionic and scalar dark matter candidates and discuss the expected values in light of current bounds from Direct Detection (DD) experiments.

Fermionic dark matter
In the case where the lightest neutral stable particle is the fermionic state f 0 , a vector coupling exists between this state and the Z−boson of the SM allowing for a Spin-Independent (SI) cross section at tree-level. The amplitude for the scattering of our DM state and firstgeneration quark is determined by the charge assignment of these particles with respect to the electroweak symmetry group, the only free parameter being the DM mass. The corresponding expression for the DD cross section with a proton can be straightforwardly expressed as [37] where G F is the Fermi constant, s W is the sine of the Weinberg angle, in the limit where the DM mass is much larger than the proton mass m f m p . As long as this hierarchy is satisfied, this expression is independent of any parameter of the model and fixed. The most stringent constraints on the SI cross section are currently achieved by the Xenon1T experiment [38] which excludes σ p 5 × 10 −46 cm 2 for a 50 GeV DM mass and up to σ p 10 −44 cm 2 for 10 TeV DM masses. As a result, by requiring the DM mass to remain below the perturbative unitarity limit O(10 − 100) TeV, the fermionic DM candidate is entirely ruled out by the current constraints on SI cross section from direct detection.

Scalar dark matter
In the case where our DM candidate is the lightest neutral scalar S 0 1 , a SI contribution to the direct detection cross section is generated at tree-level by the coupling of S 0 1 to the Higgs boson with s θ 0 ≡ sin θ 0 and c θ 0 ≡ cos θ 0 . This induces an effective coupling to nucleons (N ) of the form where the effective coupling λ eff N can be expressed as T q ≈ 0.3 being the scalar form factor. The SI DM-proton cross section can be straightforwardly derived: We checked that this expression is in agreement with Ref. [39] and Ref. [40] in the limit of vanishing mixing angles and agrees numerically with results from the micrOMEGAS code [41,42].
However, electroweakly charged DM candidates can have sizable loop-contributions to direct detection cross section induced by box diagrams involving electroweak gauge fields [37]. At the loop-level, the neutral component of a hypercharge Y scalar triplet interacts with quarks via a twist-2 operator in the large-mass expansion as [40,43,44] We use the parametrization from [40] which reduces to F (r) 4π/ √ r in the limit r 1. In our case the DM candidate is a mix of neutral scalar components of Y = 0 and Y = 1 triplets, and several quasidegenerate charged and neutral scalar states are present in the spectrum. The precise estimate of the one-loop contribution to the direct detection cross section would require a rather cumbersome and dedicated analysis which goes beyond the scope of this paper. Therefore, in the limit where the mass splitting between the various scalar components is small -assumption justified as illustrated in Appendix E -we can estimate the oneloop amplitude as the sum of the contributions from Y = 0 and Y = 1 triplets weighted respectively by cos 2 (θ 0 ) and sin 2 (θ 0 ), which should reproduce the expression from Ref. [40] in the limit of vanishing mixing angles. This gives the following estimate of the cross section by including the electroweak one-loop corrections: where f PDF p = 0.526 [45] is the second moment of proton parton distribution function (PDF) evaluated at the scale µ = m Z . In this case, the SI cross section depends on numerous parameters of the model and the numerical evaluation of this quantity is investigated in details further on in Sec. 4.3.

Dark matter relic abundance
In scotogenic models, the dark matter relic abundance can be generated by the usual freezeout mechanism. Indeed, the electroweak charge of our DM candidate ensures thermalization with the SM bath until the DM particles become non-relativistic and annihilate. In the limit where our DM candidate is Ω 0 R -like (the neutral component of the Y = 0 triplet), a large contribution to the DM annihilation cross section, induced by gauge interactions, features weak bosons in the final state such as W + W − . Such contribution can be expressed as which essentially depends only on the DM mass. This sets the typical mass scale for our DM candidate. For DM particles with masses lighter than ∼ 1.2 TeV, processes of the kind DM + DM → SM + SM are too efficient resulting in a density to small to account for all the observed DM relic abundance and must rely on additional processes. For masses above ∼ 1.2 TeV, the observed DM relic abundance can be achieved by DM annihilations to electroweak SM bosons via gauge interactions, i.e. W + W − , ZZ. In addition, annihilations to Higgs bosons pair hh via a combination of couplings κ, λ HΩ , λ H∆ and λ H∆ depend on the specific mixing between the neutral components of the Y = 0 and Y = 1 new scalar multiplets. Annihilation cross sections to SM fermions, mediated by s−channel Higgs-boson, are typically helicity suppressed and do not contribute significantly, except for the top-quark case. However, annihilations into SM leptons via the Yukawa couplings y i Ω and y i ∆ , mediated by t−channel exchange of heavy fermionic states f , can significantly contribute. As mentioned in the previous section, only one of these Yukawa coupling is a free parameter as we imposed to reproduce observed parameter of the neutrino sector.
However, as the DM candidate is the lightest and only stable Z 2 -odd particle, any heavier Z 2 -odd state is unstable and will eventually decay to a DM particle and additional states. By conservation of number density in the Z 2 -odd sector, the various new scalars (charged and neutrals) and new fermionic states will contribute to the co-annihilation cross section via 2 → 2 processes where the initial states are a pair of Z 2 -odd particles and the final states DM particles. In addition, as shown further on, the mass splitting between these states and the DM candidate is typically small implying that co-annihilations are sizable and have to be taken into account to estimate the DM relic density. The numerous diagrams involved in (co-)annihilation processes are depicted in Fig. 12 and Fig. 13 for illustration. As the number of involved states and diagrams is rather large, the relic abundance has to be computed numerically. For this purpose, we rely on the code micrOMEGAS [41,42] after implementing the model both in LanHEP [46] and Feynrules [47] as cross check.

Indirect searches
Dark matter annihilations within the halo of the galactic center or in galactic subhalos might produce γ-ray signals potentially observable by the current and future generation of telescopes, depending on the annihilation channels. In order to illustrate the relevance of various annihilation channels, represented in Fig. 2, we provide analytical expressions for the most relevant velocity-averaged annihilation cross sections assuming our DM candidate to be S 0 1 Ω 0 R , in the limit of vanishing mixing angles and assuming m S 0 We performed a velocity expansion at leading order in the mean DM velocity v DM 1. In this limit, the DM annihilation cross section into a pair of gauge bosons is given by the s-wave terms with e ≡ √ 4πα EM . Additional s−wave terms induced by the couplings λ H∆ , λ H∆ and κ are also present but suppressed by mixing angles. Annihilation cross section into a pair of Higgs bosons is given by Annihilations to a pair of charged leptons + i − i correspond to the following cross section where i = 1, 2, 3 is a flavour index. The second term being helicity suppressed, annihilations are mostly efficient for heavy leptons, i.e. τ + τ − , while annihilations into a pair of neutrinos ν i ν i are only triggered by the Yukawa coupling y i Ω as For low masses m S 0 1 3 − 4 TeV, W + W − is the most promising channel as annihilations to this final state are not suppressed by any mixing angle for S 0 1 Ω 0 R and depends only on the DM mass. In addition, for large values of the couplings of the scalar potential, the ZZ and hh final state would become equally important. The HESS experiment is currently the most sensitive to annihilations for masses m DM 100 GeV. As no signal from DM annihilations has been detected so far, a constraint of the order of σv 10 −26 − 10 −25 cm 3 s −1 has been set by the HESS collaboration [48] for τ + τ − and W + W − final states. Moreover, CTA should be able to probe the vanilla value of the velocity averaged annihilation cross section σv = 3 × 10 −26 [49][50][51][52] for DM annihilations into specific channels such asbb, τ + τ − and W + W − . For larger masses m S 0 1 3 − 4 TeV, leptonic final state could dominate for sizable values of the Yukawa couplings. CTA sensitivity for DM annihilations to neutrinos should reach σv νν ∼ 10 −24 cm 3 s −1 for m DM 1 TeV [53]. However, the KM3Net experiment should reach the best sensitivity for such DM masses, by probing values as low as σv νν 10 −25 cm 3 s −1 [54,55]. More details regarding our numerical results and the effect of such bounds are discussed in the following section.

Lepton flavour violation
One of the typical signatures of scotogenic-like models is rare lepton-flavour-violating leptonic decays [56,57] such as α → β γ or α → α β β . One of the strongest bounds on such processes impose BR(µ + → e + γ) < 4.2×10 −13 [58] and BR(µ + → e + e + e − ) < 1×10 −12 [59]. In our setup, α → β γ is generated at the loop level by the following effective operator The branching fraction for µ → eγ can be expressed as [56] BR The diagrams involved in the calculation of µ µe are depicted in figure 3. The complete expression for the coefficient µ µe , computed for our model, can be found in Eq. (D.1) of the appendices and depends on many various parameters of the model. In order to provide a simple numerical estimate, we consider the limit where the mass splittings between the DM candidate and the various scalars and vector-like fermion is small. In addition, by assuming universal Yukawa couplingsȳ = y i ∆ = y j Ω and vanishing mixing angles, the branching fraction is approximately given by with e = √ 4πα em . The constraint BR µ → eγ < 4.2 × 10 −13 translates into a bound One-loop contributions to the process µ → 3e arise from box diagrams as well as γ, Z and Higgs-penguin diagrams. The Higgs diagrams, involving suppressed yukawa couplings, are negligible for the first and second generation of leptons. The Z−penguins have been shown [56] to be suppressed by the charged lepton masses and therefore are subdominant. In case where the dipole-like contribution to the photon penguin diagram is dominant, the decay rate for µ → 3e becomes proportional to the µ → eγ rate and is suppressed by an additional fine structure constant and phase space volume. [56] have shown that non-dipole photonic diagrams never exceed the dipole contribution. As a result, only box diagrams can lead to a decay rate µ → 3e larger than µ → eγ. The ratio of these two contributions would be suppressed by factor ∼ȳ 4 /(48π 2 e 2 ) times a ratio of the corresponding loop functions. Since from Eq. (3.23), the typical Yukawa couplingȳ is already constrained by µ → eγ to be small, we expect the rate µ → 3e to be subdominant therefore less constraining than µ → eγ. This ratio has been computed explicitly in [56], for the normal hierarchy case, and is typically of order ∼ 10 −2 for m ν 1 = 0, assumption that we made throughout this work.
4 Numerical analysis of the parameter space

Scan of the parameter space
We performed a scan in the parameter space and selected the points satisfying the relic density condition Ω DM h 2 ∈ [0.11933 ± 3 × 0.00091] as a 3σ interval around the Planck best fit value Ω DM h 2 = 0.11933 for TT, TE, EE+lowE+lensing+BAO [60], computed numerically using the micrOMEGAS [41,42] code. In addition, we imposed conditions allowing to reproduce the "normal-ordering" neutrino mass hierarchy as described in Sec. 3.1 and lepton flavour violation constraints from µ → eγ as described in Sec. 3.5. In order to perform the scan more efficiently, we defined convenient variables whose definitions can be found in Appendix B. In order to perform the scan, we generated random numbers in the following range of the various couplings sin(θ 0 ) , sin(θ 1 ) ∈ [−1, 1] , λ , κ ,ŷ ∆ ∈ [10 −10 , √ 4π] , (4.1) whereŷ ∆ is an effective Yukawa coupling relevant for the neutrino sector, as detailed in Appendix B. λ denotes all the couplings labeled as λ ( ) i with i = {∆, Ω, ∆Ω, H∆, HΩ}. We take √ 4π as perturbative limit (upper bound) for the dimensionless couplings. In addition, we imposed a relative mass splitting (m − S 0 1 )/S 0 1 > 10 −5 where m denotes the mass of any Z 2 -odd state, to ensure efficient numerical convergence. We imposed the boundedfrom-below conditions for the scalar potential, described in Sec. 2.2.4. We considered the loop-induced mass splitting between the charged f + and neutral f 0 fermions as described in Sec. 2.1. The scan is performed in log-space except for the couplings sin θ 0,1 where the scan is done in linear space on the variables ψ 0,1 defined in Appendix B. We split our scan in the parameter space in 2 regions, corresponding to m S 0 1 > 1.2 TeV and m S 0 1 < 1.2 TeV, where we run a longer scan in the former case as the relic density is much more sensitive to the various couplings of the models as in the later case.

Neutrino sector
The numerical results of our scan are represented in the plane (y i Ω , y i ∆ ) in Fig. 4. As discussed in Sec. 3.1, the neutrino masses are proportional to the product of couplings ∝ y Ω y ∆ therefore the parameter space shown in Fig. 4 corresponds to a "broad line" in log-log space roughly defined by y Ω y ∆ ∼ 10 −7 − 10 −10 depending on the flavour index. The broadness of the "line" -for a fixed value of y i ∆ -corresponds to a variation of these Yukawa couplings compensated by other relevant couplings, controlling neutrino masses such as mixing angles and masses. The Yukawa couplings are strongly constrained by the neutrino masses but could be as low as ∼ 10 −10 . However they cannot be much lower as this would imply having Yukawa couplings that reach the perturbative limit that we imposed on our scan in the parameter space. The small values y ∆ , y Ω ∼ 10 −10 are theoretically not the most appealing part of our parameter space as it is not much more "natural" than adding right-handed neutrino singlets coupled to the Higgs and left-handed neutrinos. However a large y ∆ implied by a feeble value of y Ω leads to interesting detection prospects as discussed further on. In addition, such large hierarchy y ∆ ∼ 10 ±10 y Ω is present typically for large DM masses m S 0 1 3 − 4 TeV, as illustrated in Fig. 11. Indeed in this part of the parameter space, the gauge contribution to the annihilation cross section is no longer efficient enough and in order to achieve the correct DM relic density, annihilations have to occur via the Yukawa couplings in addition to the various couplings of the scalar potential such as κ for instance which is also constrained by the neutrino masses. Combining this effect with the perturbative limit on our parameters sets an upper bound on the DM mass of the order of m S 0 1 30 − 40 TeV.

Direct detection
Values of the DM spin independent nucleon cross section, predicted by our scan in the parameter space, are depicted in Fig. 5 as a function of the mass of the DM candidate m S 0 1 . In addition, we represented constraints from the Xenon1T experiment [38] in blue in addition to the future sensitivity achievement for the upcoming Darwin experiment [61] in green and the so-called neutrino floor in red [62].
The left panel of Fig. 5 shows the tree-level contribution from Eq. 750 GeV, co-annihilations are too efficient to yield the correct relic abundance. The upper diagonal limit is set by the perturbative limit reached by dimensionless couplings of the scalar potential and the Yukawa couplings as discussed in the previous subsection. The requirement of reproducing the neutrino masses and mixings as described in Sec. 3.1 imposes a lower bound on the coupling κ around κ 10 −3 for m S 0 1 ∼ 1 TeV and κ 10 −1 for m S 0 1 ∼ 10 TeV which translates into the lower diagonal limit in Fig. 5. Almost the entire part of the parameter space lies in this triangular contour, apart from very few funnel points.
The right panel of Fig. 5 shows the sum of the tree-level and one-loop electroweak contributions from Eq. (3.13). The loop contribution generate values for the SI cross section of the order of 10 −47 − 10 −46 cm 2 which can be seen as a more clustered region around these values. However, interferences between the tree-level and electroweak contributions tend to scatter the point to lower values of the cross section while leaving unmodified points corresponding to cross section larger than 10 −47 − 10 −46 cm 2 at tree-level.
Interestingly, a sizable part of the parameter space is already excluded by the Xenon1T experiment. In addition, almost the entire parameter space lies above the neutrino-floor including a majority of points that should be probed by the Darwin experiment, implying that the model could be almost entirely tested within the next decades.

Indirect detection
We represented in Fig. 6 the points from our scan in the parameter space, allowing to reproduce the observed DM relic abundance and neutrino masses, projected in the plane (m S 0 1 , σv ) as well as the constraint from HESS [63]. Values of σv are computed numerically using micrOMEGAS [41,42]. In addition we depicted the CTA sensitivity projection [49], assuming 500h of observations towards the Galactic Center and DM annihilating to W + W − final state. Similar analyses have reached the same conclusion [50][51][52], i.e. for DM annihilations intobb, τ + τ − and W + W − , CTA should probe the vanilla value of the velocity averaged annihilation cross section σv = 3 × 10 −26 . We used the same color code than for Fig. 5, as described in Sec. 3.2. The cluster of points in Fig. 6 is roughly shaped as a "broad line" decreasing towards high DM masses. This line follows essentially the behavior of the W + W − cross section of Eq. (3.15) as a function of the DM mass. The broadness of this line is determined by coannihilations that allow for several values for σv W + W − , for a given DM mass while still satisfying the relic density condition. The spread of points Tree-level + one-loop Figure 5: Dark matter spin independent direct detection cross section as a function of the lightest neutral scalar mass. Constraints from the Xenon1T experiment [38] are shown in blue in addition to the future sensitivity achievement of the upcoming Darwin experiment [61] in green and the so-called neutrino floor [62] in red. Blue dots are excluded by Xenon1T. Green dots will be probed by Darwin. Red (black) dots correspond to points above (below) the neutrino floor.
Darwin above ν-floor below ν-floor Figure 6: Dark matter perturbative velocity averaged annihilation cross section into W + W − and constraints from HESS [48] as well as sensitivity prediction for CTA [49]. The various points represent the results from our scan in the parameter space reproducing both the correct dark matter relic density and neutrino masses. The color code is the same as in Fig. 5, as described in Sec. 3.2. in particular for DM masses 1 < m S 0 1 < 4 TeV are to additional scalar-potential couplings contributing to σv W + W − such as κ or λ HΩ , as can be deduced from Eq. (3.15), being relatively large as for this part of the parameter space, as discussed previously. Fig. 6 shows that HESS is already constraining a sizable part of the parameter space which should be substantially improved by CTA. It is worth notice here that a large proportion of the points lying below the neutrino-floor in Fig. 5 seems to correspond to a relatively high value for σv W + W − that are already excluded by HESS or should be in the near future by CTA. As those points cannot be probed with the future generation of xenon-based direct detection experiments, indirect searches offer a very interesting complementary discovery prospects for this model.
In addition, as commented in Sec. 4.3, DM annihilations to leptons such as τ + τ − and νν are significant for large DM masses m S 0 1 3 − 4 TeV as the gauge annihilation channel is suppressed, as can be seen in Fig. 6. We represented the annihilation cross section to τ + τ − and νν in Fig. 7 using the same color code than Fig. 5. In addition, we represented as well the CTA sensitivity estimate of [64], assuming 500h of observation towards the galactic center on the left panel and the sensitivity prospects for KM3Net [54,55] on the right panel. Interestingly the points located below the Darwin sensitivity projection or below the neutrino floor, at large DM masses, tends to be accompanied by sizable leptonic annihilation cross section that should be probed by CTA and KM3Net.
Long range forces mediated by electroweak gauge bosons between DM states with masses m S 0 1 m W , m Z can alter the perturbative prediction for the DM annihilation cross section, the well-know effect known as Sommerfeld enhancement [65,66]. Formation of bound states, in addition to introducing new annihilation channels [67], can also manifest in resonances in the Sommerfeld enhancement, in particular for small velocities v DM 1. It has been shown in Ref. [66] that the Sommerfeld enhancement is not numerically significant at the time of freeze-out, typically occuring at m S 0 1 /T 20, therefore should not affect strongly our prediction for the relic density. However this effect could increase the cross section relevant for indirect detection by a factor of 10 2 − 10 3 [66,68,69] for a typical pure-triplet DM candidate with a mass above the TeV scale and away from resonances, but could be even larger for specific masses. This effect should be taken into account to accurately assess the viability of the model. However, given the complexity of the model and the fact that our DM candidate is a mixed state, a dedicated analysis is required to numerically estimate the Sommerfeld enhancement in a reliable way, which is beyond the scope of this paper.
As a result, the complementarity of future sensitivity achievement on the spin independent from Darwin and annihilations within the galactic halo to leptons by CTA and KM3Net offers a very interesting discovery prospects from this model addressing simultaneously the problem of the dark matter abundance and neutrino masses. In addition to searches with CTA, the future Southern Wide-field Gamma-ray Observatory (SWGO) [70] in the sourthern hemisphere could have a better sensitivity to the galactic center and other regions of the gamma-ray sky [71]. Taking into account the Sommerfeld enhancement would essentially boost the indirect detection signal and reinforce the optimistic discovery prospect message that we are addressing here.  [64], assuming 500h of observation towards the galactic center, is depicted in dash-dotted line. On the right, sensitivity prospects for KM3Net [54,55] are shown in dashdotted line. The color code is the same as in Fig. 5, as described in Sec. 3.2.

Additional constraints
Radiative breaking of the Z 2 symmetry Large fermionic masses could be responsible for breaking the discrete Z 2 symmetry by loop effects. Indeed, as shown in [72], the β-functions of the mass parameters of the scalar potential in the original scotogenic model could receive negative contributions from the fermionic states and drive the mass parameters towards negative values at some high energy scale, resulting in a breaking of the Z 2 symmetry. Depending on the Z 2 breaking scale, consistency of the low energy theory could be spoiled and affect the DM stability and density production. As shown in [73], additional Z 2 -odd scalars could help to stabilize the behavior of the β functions. The complete numerical RGE analysis lies beyond the scope of this paper. However, as in our case the DM candidate is the lightest Z 2 -odd state, such effects should occur at scales higher than the DM freeze-out temperature and would not affect the DM abundance.
Neutrinoless double beta decay Assuming normal ordering and the lightest neutrino massless, as described in Sec. 3.1, then there is only one physical Majorana phase φ in the neutrino sector. In this cases, the effective mass parameter characterizing the amplitude for neutrinoless double beta decay can be expressed as [23,28,74] where φ is the only free parameter after imposing the model to reproduce neutrino parameters described in Sec. 3.1. Numerical values for this parameter are within the range m ββ ∈ [0.0015, 0.0038] eV which is below the sensitivity of the KamLAND-Zen experiment [75] and below the sensitivity that could be achieved by the nEXO experiment in the future m ββ 6 meV [76].

Oblique parameters
The new scalar and fermionic multiplet contribute to the vacuum polarization of electroweak gauge bosons and shift their SM expected values, which translates into a shift of the oblique parameters. In particular in scotogenic-like models, deviations from the T −parameter can be sizable for large mass scales and mass splittings. Assuming vanishing mixing angles and small relative mass splitting between the charged and neutral components δ ≡ (m S +(+) − m S 0 )/m S 0 1, the contribution from the Y = 1 scalar triplet to the T −parameter is of order [26,77] T δ 2 m 2 which remains below the current experimental uncertainty T = 0.03±0.12 [28] for our viable parameter space. Numerical values of the mass splitting can be found in Appendix E.
In addition, as detailed in Appendix C, the mass splitting for the neutral and charged components of the vector-like lepton doublet is typically small ∆m f ∼ 10 −4 m f for m f = 1 TeV, therefore does not contribute significantly to the T −parameter.
Collider constraints Charged scalars can affect the effective coupling involved in the h → γγ process [78] which is constrained to be around the SM-expected value. Since the viable part of our parameter space features many scalars with similar masses typically between ∼ 800 GeV and ∼ 30 TeV, this model offers interesting complementary signatures at high luminosity colliders. For instance in proton-proton collisions, a DM pair could be produced in association with two or four charged leptons via pp → γ * , Z * →f ± * f ± * → + − S 0 1 S 0 1 or pp → γ * , Z * → S ±± * S ±± * → + + − − S 0 1 S 0 1 where the DM pairs would be interpreted as missing transverse energy at collider. However, notice here that the typical S ++ → + + signal expected from models extending the TII-seesaw [27] should not be present in our case as the Z 2 symmetry is unbroken and forbid processes with a even number of Z 2 -odd states. On the other hand, the following processes are present in our model S ++ → + + S 0 i and S ++ → f + + which somehow could be interesting to analyze. Analyses of such processes have been performed in similar models [23,[79][80][81][82]. However, their results cannot be applied directly to our case as the quantum numbers of the particles present in the theory are different. Such processes in our model would require a dedicated investigation which goes beyond the scope of this paper and is left for future work.

Conclusion
In this work, we investigated the phenomenology of a model that provides a Dark Matter candidate and a mechanism to generate neutrino masses using as inspiration the fields in the Type-II seesaw. In this model, we introduced, in addition to the SM fields, a pair of SU (2) L scalar triplets with hypercharge Y = 1 and Y = 0 and a SU (2) L doublet vectorlike fermion. All the extra fields are considered as charged under a discrete Z 2 symmetry. Three of these particles (two scalars S 0 1 , S 0 and one fermion f 0 ) are electrically neutral and depending on the specific parameters, could be the lightest state charged under the Z 2 symmetry, i.e. a WIMP-like dark matter candidate.
Various phenomenological aspects of this model were investigated. We showed that light neutrino masses, generated by loop contributions from Z 2 −odd states, could be accommodated in this model, in addition to the flavour structure reproducing the observed neutrino oscillations patterns. In this setup, two neutrinos acquire non-zero masses while one remains massless. We focused on reproducing the normal hierarchy for neutrino masses but inverted hierarchy is not excluded by our analysis and could still be viable in this framework.
We identified the points in the parameter space allowing to reproduce the dark matter relic abundance as observed by the Planck collaboration. In addition to studying the theoretical consistency of the model, we explored dark matter direct detection and indirect searches signatures. In addition, we considered constraints from the lepton flavour violating process µ → eγ. Part of the parameter space, satisfying the correct relic density condition and µ → eγ constraints, is in tension with current bounds from gamma-ray observations from HESS (W + W − and τ + τ − annihilation channels) and DM-nuclei scattering bounds from the Xenon1T experiment. We showed that among the three possible dark matter candidates, only one of them, the CP-even scalar S 0 1 , satisfies these constraints. The phenomenological viable mass for this DM candidate ranges from 700 GeV to 30 TeV. We identified points in the parameter space within the reach of the next generation of direct detection and indirect detection experiments. In particular, we showed that a large proportion of the parameter space should be probed by the Darwin experiment while a sizable part should evade these bounds but would still remain above the neutrino floor, i.e. accessible with xenon-based experiments in the future. Nevertheless, a minority of these points lie below the neutrino floor. Interestingly, an important proportion of these points lying below the neutrino floor would give rise to W + W − annihilation signals that should be observed by CTA and perhaps the Southern Wide-field Gamma-ray Observatory in the future.
In addition, we showed the complementarity between direct detection and indirect signals from DM annihilations to τ + τ − by CTA and to νν by KM3Net. Such various and complementary signatures would allow to verify or refute the model in the future, but also to discriminate between this model and others, in case some signal is observed. The Sommerfeld enhancement was not considered in this model but is necessary to predict more accurately the expected indirect detection signal. Nevertheless, It would only reinforce the optimistic detection prospect that we are addressing here and this computation is left for future work. Moreover, the model presents features that would give rise to additional interesting signatures at present and future colliders such as new charged scalars and vector-like fermions. In addition, the new states considered in this setup could impact additional lepton-flavourviolating observables and offer a complementary way of probing the model considered in this work. The detailed analysis of these signatures is left for future work.
Fonseca for discussions about BFB conditions, and Avelino Vicente for useful discussions about lepton-flavour-violation. The work of MP was supported by the Spanish Agencia Estatal de Investigación through the grants FPA2015-65929-P (MINECO/FEDER, UE), PGC2018-095161-B-I00, IFT Centro de Excelencia Severo Ochoa SEV-2016-0597, and Red Consolider MultiDark FPA2017-90566-REDC. MP would also like to thank the Paris-Saclay Particle Symposium 2019 with the support of the P2I and SPU research departments and the P2IO Laboratory of Excellence (program "Investissements d'avenir" ANR-11-IDEX-0003-01 Paris-Saclay and ANR-10-LABX-0038), as well as the IPhT for hospitality during part of the realization of this work. RL thanks the Instituto de Física Teórica (UAM/CSIC) for the hospitality. RL was supported by Universidad Católica del Norte through the Publication Incentive program No. CPIP20180343 and CPIP20200063.

A Rotation matrices
In this appendix we present the parameterization used for the diagonalization of the CPeven neutral scalar and charged scalars.

A.1 CP-even neutral scalars
In the gauge basis (Ω 0 R , ∆ 0 R ), the mass matrix is diagonalized through a basis rotation in such a way that where the mass matrix eigenvalues are ordered following m 2 where the most common parameterization for rotation matrix is in base on trigonometric functions For the case of the CP-even neutral scalar, the mixing angle θ 0 is obtained from the relation: However the mixing angle in the latter expression could be shifted in π to ensure on the ordering of the mass matrix eigenvalues. Using this parameterization, the limit when θ 0 mod (π) → 0, implies that S 0 1 → Ω 0 R and S 0 2 → ∆ 0 R ; and the opposite case occurs when θ 0 mod (π) → π/2 implying S 0 1 → ∆ 0 R and S 0 2 → Ω 0 R .
Complementary, another parameterization for the rotation matrix R 0 can produce sorted eigenvalues if we assume: and if the CP-even neutral scalar mass matrix in the gauge basis Eq. (2.10) is written as follows: where s κ = ±1 is the sign of κ present in the mass matrix. Notices that A 0 , B 0 , and C 0 are positive. Under these assumptions, we get that always produces sorted eigenvalues: Using this parameterization, ψ 0 = 1 implies the alignment S 0 1 → Ω 0 R and S 0 2 → ∆ 0 R and ψ 0 = −1 implies S 0 1 → ∆ 0 R and S 0 2 → Ω 0 R . The connection between both parameterizations is straightforward: A.2 Charged scalars with Q = 1 Similarly to the neutral scalar, the mass matrix in the gauge basis (Ω ± , ∆ ± ) is diagonalized by a rotation: where the transformation between the gauge basis and the mass basis is: (A.13) The R 1 parameterization in terms of trigonometric functions is where the mixing angle can be expressed as: Notice that the mixing angle θ 1 might have a π shift due to the ordering of the eigenvalues.
Complementary to the previous parameterization, the mass matrix (Eq. 2.15) can be written as: where A 1 , B 1 , and C 1 are positive; and s κ is the sign of κ. The corresponding rotation matrix must be: Notice the difference by a relative sign with respect to the neutral-scalar case. Here, it is easy to see that: produces ordered eigenvalues such that: The connection between both parameterizations is straightforward:

B Scan strategy of the parameter space
In this appendix we detail how the various relevant parameters of the Lagrangian can be expressed in terms of a set of convenient variables introduced to perform efficiently a scan over the parameter space.

B.1 Scalar sector
Using parametrization of the ration in terms of ψ 0,1 defined in the previous section, we define the convenient dimensionless quantities which allows to express the coupling and relates neutral scalar mass eigenvalues by the relations and Mass parameters of the scalar potential can be expressed in terms of these parameters as As a result, the parameters of the scalar potential of Eq. (2.9) can be expressed in terms of the set of variables {κ, s κ , λ Ω , λ ∆ , λ ∆ , λ HΩ , λ H∆ , λ ∆Ω , ψ 0 , ψ 1 }.

B.2 Neutrino sector
In order to reproduce neutrino masses and mixing angles as described in Sec. 3.1, we define the quantity φ N ≡ arctan ∆m 2 21 /∆m 2 32 1/4 in order to express the Yukawa couplings y i Ω and y i ∆ with i = 1, 2, 3 in term of the quantitiesŷ Ω andŷ ∆ wherê by the relations where the rotation matrix R ν is written in terms of the observed neutrino mixing angles: With this parametrization, neutrino masses can be expressed as The loop function F loop can be expressed as (B.14) Notice here that beside the scalar and vector-like state masses m f ,ŷ ∆ is the only free parameter of the neutrino sector.
To summarize, the set of independent variables used explicitly for the scan is {κ, s κ , λ Ω , λ ∆ , λ ∆ , λ HΩ , λ H∆ , λ ∆Ω , φ 0 , φ 1 ,ŷ ∆ , m S 0 1 , m f } . (B.15) C Vector-like fermions: mass splitting The neutral and charged fermions belonging to the doublet f are degenerated with a mass m f at tree-level. However a small mass splitting is generated by radiative corrections [83]. The mass splitting corresponds to where α EM is the fine structure constant, m Z is the Z−boson mass, and m f is the tree-level mass of vector-like fermions. The Π function can be expressed as Π (r) ≡ In our model, contributions to µ → eγ at the loop level are given by diagrams where a photon is emitted by a charged particle running in the loop as shown in Fig. 3. There are 8 diagrams where the charged particles in the loop can be S ± 1,2 , S ±± , f + , (f + ) c . The functional dependence of the diagrams for which the photon is emitted by a scalar or fermionic states are encoded by the functions F and G respectively and were derived using the PackageX code [84,85]. The coefficient of the operator responsible for µ → eγ, appearing in Eq. where the function F is defined as F (x) ≡ 2x 3 + 3x 2 − 6x 2 log(x) − 6x + 1 6(x − 1) 4 , (D. 2) and the function G is related to the function F by

E Additional plots
In this part of the appendix, we provide useful additional plots corresponding to our numerical scan in the parameter space as described in Sec. 4, to illustrate the interplay between the various couplings and observable quantities of this model. In Fig. 8 we depicted the perturbative velocity averaged annihilation cross section into a pair of Higgs or Z-bosons of the SM. In Fig. 9 and Fig. 10, we represented the mass splitting between various scalar or fermionic mass eigenstates and the DM mass. In Fig. 11, we show predicted values of the Yukawa couplings y 1 Ω and y 1 ∆ .     F Diagrams for dark matter production In Fig. 12 and Fig. 13 we represented the numerous diagrams involved in annihilations and co-annihilations of Z 2 -odd states, relevant for the DM density production. [10] R. N. Mohapatra and G. Senjanovic, Neutrino Mass and Spontaneous Parity