Dark matter in a charged variant of the Scotogenic model

Scotogenic models are among the most popular possibilities to link dark matter and neutrino masses. In this work we discuss a variant of the Scotogenic model that includes charged fermions and a doublet with hypercharge $3/2$. Neutrino masses are induced at the one-loop level thanks to the states belonging to the dark sector. However, in contrast to the standard Scotogenic model, only the scalar dark matter candidate is viable in this version. After presenting the model and explaining some particularities about neutrino mass generation, we concentrate on its dark matter phenomenology. We show that the observed dark matter relic density can be correctly reproduced in the usual parameter space regions found for the standard Scotogenic model or the Inert Doublet model. In addition, the presence of the charged fermions may open up new regions, provided some tuning of the parameters is allowed.


Introduction
The Standard Model (SM) of particle physics fails to address two of the most important open questions in current fundamental physics: the origin of neutrino masses and the nature of dark matter (DM). While these two issues might be completely independent, and the explanation to the DM puzzle might even come from a completely different branch of physics, it is tempting to explore extensions of the SM in which they are simultaneously addressed.
Many radiative neutrino mass models are good examples of such extensions. In this class of models, pioneered in [1][2][3][4], neutrino masses vanish at tree-level but become nonzero once loop corrections are included. This naturally explains the smallness of neutrino masses, which get suppressed by the usual loop factors. In many cases, they are induced  respectively, and two pairs of singlet vector-like fermions ψ a L,R (a = 1, 2) with hypercharge −1. The scalar doublets of the model can be decomposed into SU(2) L components as Here H is the usual SM Higgs doublet. As in the Scotogenic model, we impose an exact Z 2 parity. All the new particles are odd under this symmetry, while the SM particles are assumed to be even. The particle content of the model is summarized in Table 1. We note that H and η are only distinguished by their Z 2 charges. In fact, the scalar doublet η has the same quantum numbers as the usual inert doublet present in the Inert Scalar Doublet model [15] and the Scotogenic model [7]. The most general Yukawa Lagrangian involving the new particles can be written as where M ψ is a 2 × 2 vector-like (Dirac) mass matrix, which we take to be diagonal without loss of generality, while Y L and Y R are dimensionless 3 × 2 complex matrices. The most general scalar potential is given by where µ 1 , µ 2 and µ Φ are parameters with dimension of mass and λ j (j = 1, 2, 3, 4, 5), λ Φ , ρ 1 , ρ 2 , σ 1 , σ 2 and κ are dimensionless. We note that in the presence of a non-zero κ term it is not possible to define a conserved lepton number. As shown below, this plays a crucial role for the generation of neutrino masses.

Scalar sector
Let us now discuss the resulting scalar particle content of the model. We will assume a minimum of the potential characterized by the vacuum expectation values (VEVs) with v 246 GeV the usual SM VEV. This vacuum preserves the Z 2 parity, which remains unbroken after electroweak symmetry breaking. This forbids the mixing between H, even under Z 2 , and the η and Φ doublets, odd under Z 2 . The neutral H 0 and η 0 states can be split into their CP-even and CP-odd components as The CP-even state h can be identified with the SM Higgs boson, with m h ≈ 125 GeV, while A is the Goldstone boson that becomes the longitudinal component of the Z boson.
Assuming that CP is conserved in the scalar sector, η R and η I do not mix. Their masses are given by The with The gauge eigenstates η + and Φ + are related to the mass eigenstates H + 1 and H + 2 as where c χ = cos χ, s χ = sin χ and χ is a mixing angle. The Lagrangian mass term can be written in terms of the mass eigenstates H + 1 and H + 2 as This allows us to obtain the mixing angle χ. Combining Eqs. (8) and (13) with Eq. (12) one finds Therefore, the contribution to (m ν ) αβ of the diagram displayed in Fig. 1 is proportional to Y L αb Y R βb , with b the fermion index. It is easy to realize that there is a mirror contribution proportional to Y L βb Y R αb , obtained after exchanging the scalars. In summary, the neutrino mass matrix can be written as with (m ν ) The addition of both contributions cancels out the divergence (∆ ), as expected. We can finally use Eq. (14) to replace the mixing angle χ in Eqs. (19) and (20). By doing this, and also by assuming the hierarchy m 2 This expression matches the result in [25]. Neutrino masses are proportional to the κ parameter, as expected from symmetry arguments. In the limit κ → 0 lepton number is restored, and this explains that Majorana neutrino masses can only be induced with κ = 0. Furthermore, κ 1 is natural, in the sense of 't Hooft [28]. As in the Scotogenic model, small neutrino masses can be naturally generated in this model. One can generate m ν ∼ 1 eV with Y R,L ∼ 1, m H ± 1 ∼ 300 GeV, m H ± 2 ∼ 400 GeV and m ψ a ∼ 1 TeV if κ ∼ 10 −12 .

Dark matter candidate
As a consequence of the conservation of the Z 2 parity, the lightest Z 2 -odd state is stable and cannot decay. In contrast to the Scotogenic model, the only potentially viable candidate is one of the neutral η scalars, either η R or η I , since the other Z 2 -odd states are electrically charged. Their mass difference is controlled by the λ 5 quartic coupling, as shown in Eq. (7), Another difference with respect to the Scotogenic model is that this mass difference can in principle be large. This is because the λ 5 coupling does not enter the neutrino mass formula, see Eq. (21), and can be ∼ O (1). Also, we note that the sign of λ 5 determines the DM candidate. η I is the DM candidate if λ 5 > 0, while λ 5 < 0 selects η R as the DM candidate.
The η R and η I fields couple to the electroweak gauge bosons, since they originate from an SU(2) L doublet with hypercharge 1/2. Their production in the early Universe is thus expected to be generally dominated by gauge interactions. The exception to this rule will be found when the mass of the DM candidate is close to ∼ m h /2, when the so-called Higgs portal will be the most important channel. In general, the DM phenomenology of the model is expected to be similar to that of the Inert Scalar Doublet model [15], a popular and economical model for DM [16,17]. Finally, we point out that the DM candidates, η R and η I , carry the quantum numbers of a supersymmetric sneutrino, except for lepton number.

Analysis and experimental bounds
We now proceed to discuss the approach followed in our numerical analysis and the experimental bounds considered. The first step has been the implementation of the model in SARAH (version 4.11.0) [29], a Mathematica package for the analytical evaluation of the model. 3 With the help of this tool, we have created a SPheno (version 4.0.2) [31,32] module with the required numerical routines to obtain the particle spectrum and compute several observables of interest in our model. This includes the calculation of flavor violation observables with FlavorKit [33]. Finally, we have used micrOmegas (version 5.0.9) [34] for the evaluation of DM observables, such as the DM relic density and direct and indirect detection predictions. We performed a numerical scan with ∼ 10000 points. Our choice of parameters is summarized in Tab. 2. In particular, we choose the negative sign for λ 5 so that throughout our analysis η R plays the role of the dark matter. Let us recall that this is a choice just for definiteness, equivalent results would be obtained by assuming η I as the lightest Z 2 -odd state. Furthermore, we take µ 2 2 = µ 2 Φ in most of the parameter space to guarantee that η R (and not one of the charged states from Φ) is the lightest Z 2 -odd particle. An analogous reason lies behind the range chosen for µ 2 Φ .
The new particles introduced in this variant of the Scotogenic model may lead to different experimental signatures and affect the SM prediction of several observables. For this reason and throughout our analysis we have considered a list of experimental constraints.
Neutrino oscillation data The generation and smallness of neutrino masses is one of the main motivations behind the idea of Scotogenic models. In our analysis, we demand compatibility of the neutrino oscillation parameters with the most recent neutrino oscillation global fit [35]. This is achieved by adjusting the entries of the Yukawa matrices Y L and Y R by means of the master parametrization [36,37], which allows one to write where Σ is a diagonal 2 × 2 matrix containing the positive singular values of the matrix M , with M ij = ω i δ ij and Moreover, is a global factor, U is the leptonic mixing matrix, a unitary 3 × 3 matrix that brings the neutrino mass matrix to a diagonal form as while the matrixD √ m is defined as Here v can actually be replaced by any arbitrary scale and P depends on the neutrino mass hierarchy. In case of normal ordering, P is just the identity matrix, while for inverted ordering P is a permutation matrix that exchanges the first and third elements. Finally, the matrices W , A and B are defined in Appendix B, where explicit analytical results for the elements of the Y L and Y R matrices are also given. In our numerical analysis we assume a normal mass ordering for light active neutrino masses and take vanishing CP violating phases for simplicity. We also assume that the lightest neutrino is massless.
Lepton flavor violation While this model could in principle produce signals at facilities looking for LFV processes, these have not been observed yet. Hence we can use current null searches for LFV to constrain the parameters of the model, in particular κ which determines the magnitude of the Yukawa matrices Y L and Y R . We consider the following most stringent bounds on rare LFV processes: BR(µ → eγ) < 4.2×10 −13 [38], BR(µ → eee) < 1.×10 −12 [39] and CR(µ − , Au → e − , Au) < 7 × 10 −13 [40].
Electroweak precision observables The experimental accuracy of electroweak observables can be sensitive to the presence of new particles, like those introduced in this model. The main contribution to the higher-order calculation of electroweak precision observables is parameterized via the δρ parameter. We require an adequately small deviation of the ρ parameter from one through the following prescription: −0.00022 δρ 0.00098 [41] (3σ range).
Dark matter searches We assume our DM candidate η R to be in thermal equilibrium with the SM particles in the early Universe. We also assume a standard cosmological scenario. If no other dark matter candidates are present, then the relic abundance of η R must be in agreement with the latest observations by the Planck satellite [42], which set a limit on the cosmological content of cold dark matter: 0.1164 ≤ Ω η R h 2 ≤ 0.1236 (3σ range). The relic abundance of η R can also be subdominant, i.e. Ω η R h 2 < 0.1164, however in this case another DM candidate is required to explain the totality of the cosmological dark matter. Moreover, η R can be probed at dark matter experiments like direct detection facilities. We apply the most stringent bound on the WIMP-nucleon spin-independent (SI) elastic scatter cross-section from the XENON1T experiment [43]. We compute the direct detection cross section at tree-level. Since a more exhaustive study of this observable is out of the scope of this paper, we have considered the constraint obtained for the inert Higgs doublet [44] in order to avoid sizable loop corrections. If η R annihilates into SM particles with a cross section typical of WIMPs, it may also be detected indirectly. We consider both γ rays and antiprotons as annihilation products and compare with the respective current bounds on the WIMP annihilation cross section set by the Fermi Large Area Telescope (LAT) satellite [45], the ground-based arrays of Cherenkov telescopes H.E.S.S. [46] and the Alpha Magnetic Spectrometer (AMS-02) [47][48][49] onboard the International Space Station.

LHC searches
The new charged particles in the model can be copiously produced and detected at the LHC, see [25] for a discussion. It is however beyond the scope of this work to perform a detailed collider study of the model. In order to guarantee compatibility with the current bounds, we have chosen m ψ 1,2 > 1 TeV in our numerical scan. The singly-charged scalars H + 1,2 have masses in a wide range of values, always in the hundreds of GeV. Finally, the doubly-charged scalar Φ ++ is chosen to be always heavier than ∼ 200 GeV by imposing This may seem as too little restrictive, since the currently most stringent bounds on a doubly-charged scalar range from 220 GeV [50] (for a Φ ++ that decays into W + W + ) to 846 GeV [51] (for a Φ ++ that decays into µ + µ + ). However, notice that in our model Φ ++ is Z 2 -odd and its decays always include DM particles. In fact, we have found that in many parameter points Φ ++ decays as Φ ++ → H + 1 W + , followed by H + 1 → η R + ν, thus leading to large amounts of missing energy in the final state.
Invisible decay width of the Higgs boson If the new neutral scalars η R and η I are light enough, new invisible decay channels of the Higgs boson will be kinematically accessible.
We require that their contribution to the invisible decay width of the Higgs boson is not larger than the currently most stringent experimental bound: BR(h → inv) 19% [52]. Our computation of BR(h → inv) includes the decay BR(h → η R η R ), when kinematically possible, and also BR(h → η I η I ) when m η I − m η R < 10 GeV. We estimate that for larger mass differences between η I and η R , η I would decay visibly inside the detector. Moreover, the presence of charged scalars may also modify the Higgs photon coupling to two photons. For this reason we further impose that 0.84 η I lifetime Searches for long-lived particles at the LHC usually look for displaced vertices. In order to ensure compatibility with the non-observation of such signatures, we require τ (η I ) 1 s. This translates into the constraint |λ 5 | ≥ 5 × 10 −4 , which guarantees that η I decays fast enough into η R .
Theoretical considerations We require that the expansion of the scalar potential (see Eq. 3) around its minimum must be perturbatively valid. At this scope, we take all scalar quartic couplings to be 1. Similarly, the elements of the Y L and Y R Yukawa matrices are also bounded to be 1.
We summarize in table 3 the list of experimental constraints applied in the numerical scan.

Numerical results
In this section we summarize our results of the analysis of η R as a dark matter candidate in the model. As we previously commented, this is the only potential DM candidate in this model, given the choice of negative sign of λ 5 (see Sec. 2.3). Fig. 2 shows the expected η R relic abundance as a function of the mass of η R . Points colored in magenta denote solutions which can reproduce the observed cold dark matter relic density, as they fall within the 3σ range derived by the Planck satellite data [42], Ω η R h 2 = 0.120 ± 0.0036 (green band). From this analysis we can identify that the preferred DM mass range lies around 500 GeV m η R 800 GeV. Blue points also refer to allowed solutions but where η R would be a subdominant component of dark matter, and another candidate would be required. Finally, gray points denote solutions excluded by any of the constraints previously discussed in Sec. 3. In particular, the Planck constraint itself rules out most of the solutions leading to excessively large relic density. Another large set of solutions with 70 GeV m η R 700 GeV is excluded by the current bound on WIMP-nucleon SI elastic scattering cross section set by the experiment XENON1T [43].
Looking at the figure from left to right, the first dip refers to the Z-pole, where m η R ∼ M Z /2 and annihilations via s-channel Z exchange are relevant. A similar feature appears at m η R ∼ 60 GeV, where a narrower dip indicates efficient annihilations via s-channel Higgs exchange. However, most solutions with m η R 60 GeV appear to be in conflict with current collider limits on BR(h → inv) and with LHC searches for doubly-charged scalars and are therefore excluded. Notice that the s-channel Higgs annihilations are in general more efficient than the Z-mediated ones, not being momentum suppressed, and in principle could lead to solutions with Ω η R h 2 < 10 −5 . Only very few solutions eventually survive in the Higgs-pole region at m η R 100 GeV, due to a variety of constraints, among which BR(h → γγ) and BR(h → inv) are the most important ones. These parameter points have very small relic density. After the Higgs pole, η R quartic interactions with gauge bosons and heavy quarks become effective, when kinematically allowed. In particular, η R annihilations into W + W − via quartic couplings are efficient at m η R 80 GeV, translating into a third drop in the relic abundance. As soon as kinematically allowed, η R can then annihilate also into two Higgs bosons. In general, for m η R 100 GeV, the η R relic density increases as its annihilation cross section drops proportionally as ∼ 1 A very interesting feature appears at m η R ∼ 1.3 TeV. The relic density suddenly drops to very small values due to a very efficient co-annihilation channel H − 1 H + 2 → νν, mediated by the vector-like fermions ψ 1,2 in t-channel, as shown in Fig. 3. This contribution is made very efficient by the magnitude of the relevant entries of the Yukawa matrices entering the diagram. The same Yukawas also play a relevant role in the rare LFV process µ → eγ, which turns out to exceed the current limit set by MEG [38] and thus excludes all solutions falling in this very narrow pole. Nonetheless, given the large freedom due to the vast array of parameters that can be varied in the parameterization of the Yukawa matrices (see Appendix B), one can always find a fine-tuned combination which allows to keep the co-annihilation channel H − 1 H + 2 → νν efficient while at the same time not leading to a BR(µ → eγ) in conflict with current observations. This can be seen in Fig 4, where we show the dependence of BR(µ → eγ) and Ω η R h 2 as a function of the parameters K 12 (left panel) and T 12 (right panel). To obtain these figures we extracted one solution at m η R ∼ 1.3 TeV from our general scan (see Fig. 2) -which we recall was made fixing K 12 = 0 and T 12 = 0 -and we scanned around those central values. Eventually, few allowed solutions (both with a relic abundance matching current observations and under-abundant) will appear at m η R ∼ 1.3 TeV. Let us notice that similar features with fine-tuned solutions may be present also in other parts of the parameter space.
Next we proceed to discuss the results for η R direct detection. In full generality, the tree-level SI η R -nucleon interaction cross section receives two contri-  Figure 3: Dominant diagram contributing to the co-annihilation channel H − 1 H + 2 → νν. Notice that both Yukawa matrices, Y L and Y R , must have sizable entries to make this process efficient.
butions, from scattering through the Higgs or Z bosons. However, since the λ 5 quartic coupling induces a mass splitting between η R and η I , eventually it turns out that the interaction through the Z boson is kinematically forbidden, or leads to inelastic scattering, in a large part of the parameter space. Should this not be the case, the η R -nucleon interaction  [38]. Therefore, all blue points above this line are excluded.
through the Z boson would dominate (since the η doublet has hypercharge different from zero) and very likely exceed the current bounds from DM direct detection experiments. We show in Fig. 5 the SI η R -nucleon elastic scattering cross section, weighted by ξ = Ωη R Ω CDM,Planck versus the η R mass. We use the same color code as in Fig. 2. The plain green line and dashed area indicate the current most stringent limit from XENON1T [43]. Even if not shown here, other stringent constraints on the SI η R -nucleon elastic scattering cross section apply from the liquid xenon experiments LUX [53,54] and PandaX-II [55]. Liquid argon experiments like DarkSide-50 [56] and DEAP-3600 [57] are instead presently limited and hence less constraining, due to lower exposures and the currently low acceptance in DEAP-3600. We also illustrate for comparison the expected discovery limit corresponding to the "ν-floor" from coherent elastic neutrino-nucleus scattering (CEνNS) for a Ge target [58] (dashed orange line), as well as the sensitivity projection (90% CL) for the future experiment LUX-ZEPLIN (LZ, red dashed) [59]. Other future experiments like XENONnT [60], DarkSide-20k [61], ARGO [61] and DARWIN [62,63] (see [64] for a recent review), not shown here not to overcrowd the figure, will in principle be able to explore the parameter space of this model. Given the current experimental constraints, most of the solutions with a relic abundance falling within the 3σ C.L. cold dark matter range obtained by the Planck collaboration [42] lie in a narrow region around m η R ∼ 500 − 700 GeV. Many more solutions leading to underabundant dark matter lie at ∼ 150 GeV m η R ∼ 500 GeV and could be tested at future DD experiments. On the whole, the phenomenology of the scalar DM candidate in this model, while richer, shares similar properties to that of the analogous scalar DM candidates in the simple Scotogenic Model [7], in Inert Scalar Doublet models, and in other Scotogenic variants like the Singlet-Triplet Scotogenic model [24].
Another η R search strategy is via its indirect detection. If it annihilates into SM particles or messengers (with sizable annihilation cross section, close to the thermal value) it can contribute to the flux of cosmic particles that reaches the Earth. Photons, and more specifi-  Figure 5: Spin-independent η R -nucleon elastic scattering cross section -weighted by the relative abundance -versus m η R . Same color code as in Fig. 2. The green dashed area is excluded by the most recent upper bound from the XENON1T experiment [43]. The dashed orange curve indicates the expected discovery limit corresponding to the "ν-floor" from CEνNS of solar and atmospheric neutrinos for a Ge target [58]. The dashed red curve depicts the projected sensitivity expected at LUX-ZEPLIN (LZ) [59].
cally γ rays are among the most suitable messengers to probe WIMP dark matter indirectly. They are not deflected during propagation, so they carry information about their source, and they are relatively easy to detect. When η R is lighter than M W , only annihilations into fermions lighter than m η R are allowed at tree level, with the heaviest kinematically allowed fermion final states dominating (i.e. η R η R → bb and η R η R → τ + τ − ). At higher masses, when kinematically allowed through the corresponding quartic couplings, also the following channels open: η R η R → W + W − , hh, Z 0 Z 0 . The hadronization of the gauge bosons, Higgs boson or quarks produces neutral pions, which in turn decay into photons thus giving rise to a γ-ray flux with a continuum spectrum which may be detected at indirect detection experiments. Besides this featureless γ-ray spectrum, the model also predicts a spectral feature from the internal bremsstrahlung process η R η R → W + W − γ, similarly to the Inert Doublet model [65,66]. We include this contribution in our analysis.
While more challenging due to uncertainties in the treatment of their propagation, charged particles can also be used to probe for η R annihilations. PAMELA [67,68] and, more recently, the AMS-02 [69,70] positrons data allow to place constraints on annihilating WIMPs, which are particularly stringent if they annihilate mainly to the first two generations of charged leptons. In our case, at m η R < m W , η R annihilates predominantly into τ s (and bb). Hence, bounds from cosmic positrons are less stringent than those from γ rays. On the other hand, AMS-02 data on the antiproton flux and the Boron to Carbon (B/C) ratio can be used to constrain the η R annihilation cross section [47][48][49]. Provided that astrophysical uncertainties on thep production, propagation and on Solar modulation are reliably taken into account (see e.g. [71][72][73]), these bounds result to be stronger than γ-ray limits from dwarf spheroidal satellite galaxies (dSphs) over a wide mass range.
We focus on the main η R annihilation channels into bb, τ + τ − and W + W − to compare with current limits set by γ-ray observations of Milky Way dSphs with Fermi-LAT data [45], of the Galactic Center (GC) with the H.E.S.S. array of ground-based Cherenkov telescopes [46] and a combination ofp and B/C data of AMS-02 [47,48]. Fig. 6 shows the results of our numerical scan of the η R annihilation cross section (weighted by ξ 2 and by the correspondent branching ratio) versus its mass, for the main annihilation channels: η R η R → bb (green points), η R η R → τ + τ − (blue points) and η R η R → W + W − (dark red points). As in previous plots, gray points denote solutions excluded by any of the constraints listed in Sec. 3. In the same figure we depict the 95% C.L. upper limits currently set by the Fermi-LAT with γ-ray observations of Milky Way dSphs (6 years, Pass 8 event-level analysis) [45] (plain curves and shaded areas, assuming annihilation into bb (green), τ + τ − (blue) and W + W − (dark red)). We also show the current upper limit obtained with data accumulated over 10 years by H.E.S.S. observations of the GC [46], assuming a W + W − annihilation channel and an Einasto dark matter density profile (dark red dot-dashed curve and shaded area). Current bounds from a combination ofp and B/C data of AMS-02 [47,48] are instead shown as dashed curves (green for bb and red for W + W − channels, respectively). As for comparison we further illustrate projected sensitivities for Fermi-LAT from a stacked analysis of 60 dSphs and 15 years of data, in the bb channel [74] (dot-dashed green line) and for the Cherenkov Telescope Array (CTA), for the Milky Way galactic halo target, assuming the W + W − annihilation channel and an Einasto dark matter density profile [75]. Our model predictions, not excluded by other constraints, lie at least a factor of few below current limits set by γ-ray observations with Fermi-LAT and H.E.S.S.. In particular, allowed solutions in the "low" (and narrow) mass window around m η R ∼ 60 GeV (bb and τ + τ − annihilation channels) are hardly found, so only a couple of them are shown for illustration. Predictions for these channels lie well below current bounds. Antiprotons and B/C data of AMS-02 instead already allow to exclude some solutions with 200 GeV m η R 800 GeV (W + W − annihilation channel). We remark that these bounds are obtained under significant astrophysical uncertainties. Future telescopes like CTA as well as additional LAT data will allow to further explore these regions of the parameter space, thus allowing for multi-messenger indirect probes of this model. Finally let us comment that the annihilation cross section of non-relativistic η R at the current epoch can be affected by a non-perturbative correction, the Sommerfeld enhancement [76][77][78][79][80] which would consequently affect also its indirect detection signatures. However, we did not include it in our calculation as a thorough treatment of this effect lies beyond the scope of this paper.  Figure 6: η R annihilation cross section for annihilations into bb (green), τ + τ − (blue) and W + W − (red) final states. The green, blue and red plain lines refer to the corresponding 95% C.L. upper limits currently set by Fermi-LAT γ-ray data from dSphs [45]. The dark red dot-dashed curve is the current 95% C.L. upper limit obtained by H.E.S.S. using GC data [46]. The green and red dashed lines denote current 95% C.L. constraints derived from the antiproton and B/C data of AMS-02 [48]. For comparison, we also show sensitivity projections for Fermi-LAT (bb channel, assuming 60 dSphs and 15 years of data) [74] and for CTA (looking at the galactic halo, W + W − channel) [75].

Conclusions
We have considered a variant of the Scotogenic model with an extended scalar content, in which one of the new doublets carries hypercharge Y = 3/2. As a consequence, the model particle spectrum contains new states, including a doubly-charged scalar, singly-charged scalars and new charged fermions leading to a rich collider phenomenology. Neutrino masses are generated radiatively as in other Scotogenic scenarios. However, contrarily to the simple Scotogenic model, in this setup only the lightest Z 2 -odd neutral scalar η is a viable dark matter candidate. We have analyzed the phenomenology of the dark matter particle η R and shown that the observed relic density can be obtained for 500 GeV m η R 800 GeV. This result is equivalent to that obtained in the Inert Doublet model. In addition, we have found that the correct relic density can also be achieved at larger η R masses if the coannihilation channel H − 1 H + 2 → νν becomes efficient. However, this requires some fine-tuning to avoid the stringent constraints from lepton flavor violation. Moreover, we have presented a full numerical analysis of the signatures expected at dark matter detection experiments, both via direct and indirect probes. We have found that most recent direct detection data from XENON1T already rule out a region of the parameter space in the mass range m η R ∼ 100 GeV. Finally we have commented about indirect detection signatures with γ-rays and antiprotons data. While bounds derived from γ-ray data from Fermi-LAT and H.E.S.S. currently lie above the model predictions, AMS-02 antiproton data instead already allows to constrain some solutions in the 200 GeV m η R 800 GeV mass range.

A Neutrino mass matrix in Scotogenic models
In this Appendix we find a general analytical expression for one-loop contributions to the neutrino mass matrix in Scotogenic models. Let us consider a generic model with n S scalars S a with masses m Sa (a = 1, . . . , n S ) and n F fermions F b with masses m F b (b = 1, . . . , n F ) that couple to light neutrinos with interaction terms of the form where Y F S is a 3 × n S × n F complex object. This interaction term induces neutrino masses at the one-loop level via the diagram in Fig. 7. The resulting neutrino mass matrix can be expressed as where the sum extends over all mass eigenstates running in the loop. Each i (m ν ) represents the individual contribution by the scalar S a and the fermion F b and is given by where D = 4 − ε is the number of space-time dimensions, k is the momentum running in the loop and the external neutrinos have been taken at rest. The / k piece does not contribute since it is an odd function in k. The remaining term is logarithmically divergent. Introducing the usual Feynman's parameters, we can rewrite our expression as where One can now integrate in k and obtain Figure 7: Generic one-loop contribution to the neutrino mass matrix. F and S represent fermion and scalar mass eigenstates contributing to the neutrino mass matrix.
where ∆ = 2 − γ + log 4π and γ is Euler's constant. Finally, we integrate in x and take the limit → 0 (which corresponds to D → 4). This leads us to Eq. (35) is our master expression for the one-loop contributions to the neutrino mass matrix. In order to compute the neutrino mass matrix in different versions of the Scotogenic model, one just needs to identify the scalar and fermion mass eigenstates running in the loop and replace the Y F S couplings by their expressions in terms of the model parameters. Then, summing over all individual contributions one finds the total one-loop neutrino mass matrix. We note that individual contributions to m ν are divergent due to the ∆ term. However, the one-loop neutrino mass matrix is finite. Therefore, these divergences cancel out in the sum in Eq. (31).
As an example, we show this explicitly in the original Scotogenic model [7]. In this case one introduces the inert scalar doublet η, with the same definition as in the model in Sec. 2, as well as three generations of fermions N , singlets under the SM gauge group. The new N and η fields couple to the SM lepton doublet with the Yukawa term Y N L η N , with η = iτ 2 η * and τ 2 the second Pauli matrix. Moreover, the neutral component of the η doublet can been split into its CP-even and CP-odd components as Assuming that CP is conserved in the scalar sector, η R and η I do not mix and constitute mass eigenstates. Therefore, in the Scotogenic model one has n S = 2 and n F = 3, with S ≡ {η R , η I } and F ≡ {N 1 , N 2 , N 3 }. Their couplings are given by Finally, applying Eq. (35) and summing over all the states in the loop we are able to obtain a finite expression (the divergences of the diagrams with η R and η I cancel each other) for the neutrino mass matrix where m η R and m η I are the corresponding masses for each component. This expression differs from the expression in [7] by a factor of 1/2, an error that was first identified in v1 of [81].

B General parametrization for the Yukawa matrices
The master parametrization [36,37] allows one to write the 3 × 2 Yukawa matrices Y L and Y R in terms of neutrino oscillation parameters as in Eqs. (23) and (24). In these expressions, W is a 2 × 2 unitary matrix, given in terms of the complex angle β as A = T C 1 , with T a 2 × 2 upper-triangular matrix with positives real entries in the diagonal, and B = (T T ) −1 (C 1 C 2 + KC 1 ), with K a 2 × 2 antisymmetric matrix, hence containing only one degree of freedom. The matrices C 1 and C 2 are given by with z 1 and z 2 two complex numbers such that z 2 1 + z 2 2 = 0. With all these ingredients one can obtain explicit expressions for the elements of Y L and Y R , which will be determined from neutrino oscillation data as well as from a set of free parameters. We point out that these are the most general expressions for the elements of the Y L and Y R Yukawa matrices: Y R 32 = − √ m 1 s 12 s 23 − c 12 c 23 e −iδ s 13 (T 22 s β (z 2 K 12 − z 1 ) + c β (z 1 K 12 T 11 + z 2 K 12 T 12 − z 1 T 12 + z 2 T 11 )) √ 2f √ ω 2 T 11 T 22 + e −i(η 2 +η 3 ) e iη 3 √ m 2 (c β (K 12 T 11 + T 12 ) + T 22 s β ) √ 2f √ ω 2 T 11 T 22 + e −i(η 2 +η 3 ) c 23 e −iδ s 12 s 13 + c 12 s 23 + c 13 c 23 e iη 2 √ m 3 (c β (T 11 − K 12 T 12 ) − K 12 In these expressions c kl ≡ cos θ kl and s kl ≡ sin θ kl (k = 1, 2 and l = 2, 3; l = k) are the usual neutrino mixing angles measured in oscillation experiments, while m j ≡ m 2 1 + ∆m 2 j1 (j = 2, 3) are neutrino mass eigenvalues. Moreover, m 1 is the lightest neutrino mass, η i (i = 2, 3) are CP violating Majorana phases and δ is the CP violating Dirac phase. We also use the notation s β ≡ sin β and c β ≡ cos β. In our numerical scans we have fixed T 11 = T 22 = 1 and β = T 12 = K 12 = z 1 = z 2 = m 1 = η i = δ = 0 unless otherwise stated.