Lepton Flavor Violation in the Scotogenic Model

We investigate lepton flavor violation in the scotogenic model proposed by Ma in which neutrinos acquire non-zero masses at the 1-loop level. Although some works exist in this direction, they have mainly focused on the radiative decay $\ell_\alpha \to \ell_\beta \gamma$. Motivated by the promising new projects involving other low-energy processes, we derive complete analytical expressions for $\ell_\alpha \to 3 \, \ell_\beta$ and $\mu-e$ conversion in nuclei, and numerically study their impact on the phenomenology. We will show that these processes can actually have rates larger than the one for $\ell_\alpha \to \ell_\beta \gamma$, thus providing more stringent constraints and better experimental perspectives.


I. INTRODUCTION
The search for lepton flavor violation (LFV) is going to live an unprecedented era with great experimental efforts in many different fronts. In addition to the well-known searches for the radiative decay ℓ α → ℓ β γ, new projects involving other low-energy processes, such as ℓ α → 3 ℓ β or µ − e conversion in nuclei, are going to look for a positive LFV signal.
For many years, the experiment leading to the most stringent constraints has been MEG [1]. This experiment, which searches for the radiative decay µ → eγ, recently published a new limit, Br(µ → eγ) < 5.7 × 10 −13 , obtained with an updated analysis of the 2009-2010 data sample together with the analysis of the new data collected in 2011 [2]. The expectation is that MEG can reduce the current bound by another order of magnitude, with sensitivities of about 6 × 10 −14 after 3 years of acquisition time [3].
Finally, the limits for τ observables are less stringent, although significant improvements are expected at B factories [10,11]. Table I  For this reason, one needs to fully understand the anatomy of LFV in each model in order to determine the expected hierarchies among observables, which then become indirect tests of the model.
In this paper we pursue this goal in the context of a model proposed by Ma in which neutrinos acquire masses at the 1-loop level [17]. The same symmetry that forbids the tree-level contribution to Dirac neutrino masses, a Z 2 parity, also gives rise to a dark matter candidate. This simple extension of the Standard Model (SM), usually called Scotogenic Model, constitutes a very simple framework to address the most important motivations to go beyond 1 . Although some works have been already done regarding LFV in this model [29][30][31][32], they have either focused on µ → eγ or neglected contributions beyond the photonic dipole. To the best of our knowledge, this is the first time ℓ α → 3 ℓ β and µ − e conversion in nuclei are fully considered. As we will see, these processes might actually have rates larger than the one for µ → eγ, thus providing better bounds and experimental perspectives.
The rest of the paper is organized as follows: in Sec. II we describe the model and its basic features. In Sec. III we present our analytical results, whereas Sec. IV contains a numerical discussion addressing some phenomenological issues of interest. Finally, we summarize our results and conclude in Sec. V.

II. THE MODEL
The model under consideration [17] adds three right-handed neutrinos N i (i = 1-3) and one SU (2) L doublet η to the SM particle content. In addition, a Z 2 parity is imposed, under which the new particles are odd and the SM ones are even 2 . The interaction of the right-handed neutrino sector is described by the Lagrangian Note that one can always write the right-handed neutrino mass term as a diagonal matrix without loss of generality. The scalar potential V is given by We assume that the parameters in the scalar potential are such that the doublet η does not get a vacuum expectation value. This is fundamental in order to keep the Z 2 symmetry unbroken. After electroweak symmetry breaking, the masses of the charged component η + and neutral component where the mass difference between η R and η I is m 2 neutrino mass matrix can be expressed as where m R and m I are the masses of η R and η I respectively, and the Λ matrix is defined as In particular, when m 2 R ≈ m 2 I ≡ m 2 0 (λ 5 ≪ 1), the mass matrix gets the simplified form This neutrino mass matrix is diagonalized as where is the PMNS (Pontecorvo-Maki-Nakagawa-Sakata) matrix. Here c ij = cos θ ij , s ij = sin θ ij , δ is the Dirac phase and ϕ 1 , ϕ 2 are the Majorana phases 4 .
The Yukawa matrix y iα can be written using an adapted Casas-Ibarra parametrization [33] as where R is an complex orthogonal matrix which satisfies R T R = 1.
3 Note that the tree-level contribution is actually forbidden by the Z2 discrete symmetry. 4 We will neglect Majorana phases in all our computations.

III. ANALYTICAL RESULTS
In this section we present our analytical results for the LFV processes ℓ α → ℓ β γ, ℓ α → 3 ℓ β and µ − e conversion in nuclei. Before we proceed to the analytical discussion a comment is in order. It is well-known that the rates for LFV processes get greatly enhanced in models with right-handed neutrinos at the electroweak scale [34][35][36][37][38][39][40][41][42]. This is due to the fact that the GIM suppression at work in the SM contribution is spoiled by the mixing between left-and right-handed neutrinos.
One could naively think that this is also the case in the scotogenic model. However, the unbroken Z 2 symmetry forbids this mixing, (see footnote 2), and thus the enhancement in the W − ν loops is not present. We will show that the enhancement is still possible, but with η ± − N loops instead.
The most popular searches for LFV have focused on the radiative process ℓ α → ℓ β γ. This is described by the effective Lagrangian where µ βα is a transition magnetic moment. It proves convenient to define it in terms of the dipole form factor A D as µ βα = em α A D /2, where terms proportional to m β have been neglected and e is the electromagnetic coupling, related to the electromagnetic fine structure constant as α em = e 2 /(4π). In the model under consideration, A D gets contributions at the 1-loop level from the Feynman diagrams in Fig. 1. They lead to the following expression where the ξ i parameters are defined as ξ i ≡ m 2 N i /m 2 η + and the loop function F 2 (x) is given in appendix A. Finally, the branching fraction for ℓ α → ℓ β γ is calculated as where G F is the Fermi constant.
Next we consider the process ℓ α → 3 ℓ β (more precisely denoted as ℓ α → ℓ βlβ ℓ β ). Although this has attracted less attention, important projects are going to be launched in the near future, with the Mu3e experiment as the leading one. There are four types of 1-loop diagrams that contribute to ℓ α → 3 ℓ β . These are γ-penguins, Z-penguins, Higgs-penguins and box diagrams. In our computations we did not consider Higgs-penguins, since we are mostly interested in processes involving the first two charged lepton generations, whose small Yukawa couplings suppress Higgs contributions. Notice that this assumption would not be valid for LFV processes involving τ leptons. However, the experimental limits in this case are not as stringent as those found for processes involving the first two generations, and thus their consideration would not change the phenomenological picture.
Let us consider the momentum assignment ℓ α (p) → ℓ β (k 1 )l β (k 2 )ℓ β (k 3 ). Then, the γ-penguin diagrams shown in Fig. 2 lead to the amplitude 5 where q ≡ k 1 − p is the photon momentum. Other operators turn out to be suppressed by charged lepton masses and thus they are neglected in Eq. (15). The coefficient A D was given in Eq. (13), whereas the coefficient A N D , which corresponds to the photonic non-dipole contributions, is given by where the loop function G 2 (x) is given in appendix A.
Similarly, we now consider the contributions from Z-penguin diagrams, also shown in Fig. 2.
Neglecting sub-dominant terms proportional to q 2 , q being the 4-momentum of the Z-boson, the resulting amplitude can be written as where are the tree-level Z-boson couplings to a pair of charged leptons. Here g 2 is the SU (2) L gauge coupling and θ W is the weak mixing angle. The coefficient F is given by Equation (19) shows that Z-penguins are suppressed by the charged lepton masses m α and m β .
Therefore, although we fully derived and included them in our computation, we found that they always have negligible contributions to the LFV processes considered in this paper. For this reason, the total decay width for ℓ α → 3 ℓ β will be mainly given by the γ-penguins and the box contributions, whose relative size will determine the phenomenology.
Finally, the box diagrams contributing to the process ℓ α → 3 ℓ β are shown in Fig. 3. One finds the following amplitude The coefficient B is given by 6 where the loop functions D 1 (x, y) and D 2 (x, y) are given in appendix A. The branching ratio for ℓ α → 3 ℓ β is given by where F RR and F RL are given by In Eq. (22), the mass of the charged lepton in the final state, m β , is kept only in the logarithmic term, where it plays the role of regulating the infrared divergence that would appear otherwise.
The most remarkable experimental projects in the near future will be devoted to searches for µ − e conversion in nuclei. The great sensitivities announced by the different collaborations might make this observable the most stringent one in most neutrino mass models. We will present our results using the notation and conventions of Refs. [45,46]. The conversion rate, relative to the the muon capture rate, can be expressed as Here Z and N are the number of protons and neutrons in the nucleus, Z eff is the effective atomic charge (see [47]), F p is the nuclear matrix element and Γ capt represents the total muon capture rate. The values of these parameters for the nuclei used in experiments can be found in [46] and references therein. Furthermore, p e and E e (taken to be ≃ m µ in the numerical evaluation) are the momentum and energy of the electron and m µ is the muon mass. In the above, g XK and g (with X = L, R and K = S, V ) are given by The numerical values of the G K coefficients can be found in [45,46,48].
Note, however, the absence of box contributions (besides the tiny SM contribution). This is due to the unbroken Z 2 symmetry, which forbids the coupling between the η ± scalars and the quark sector. Moreover, we neglect again the Higgs-penguin contributions due to the smallness of the involved Yukawa couplings. Therefore, the corresponding couplings are The photon and Z-boson couplings can be computed from the Feynman diagrams in Fig. 4. One The form factors A N D , A D and F are given in section III B, see equations (16), (13) and (19).
Furthermore, Q q is the electric charge of the corresponding quark and are the tree-level Z-boson couplings to a pair of quarks.

IV. PHENOMENOLOGICAL DISCUSSION
In this section we present and discuss our numerical results. We will explore the parameter space and highlight some relevant phenomenological issues which, to the best of our knowledge, have not been discussed in the existing literature.
In the numerical evaluation of our results we considered both hierarchies for the light neutrino spectrum 7 , normal hierarchy (NH) and inverted hierarchy (IH), and randomly chose the neutrino oscillation parameters in the 1σ ranges found by the global fit [49] (Free Fluxes + RSBL results).
We note that these ranges are in good agreement with the ones found by other fits, see Refs. [50,51].
For θ 23 , the atmospheric angle, we selected the local minimum in the first octant, in agreement with [51].
Unless explicitly expressed otherwise, all our numerical results were obtained for a degenerate right-handed neutrino spectrum, assuming a random real R matrix and λ 5 = 10 −9 . This value was found in [52] to be compatible with a correct right-handed neutrino DM relic density due to the resulting size of the Yukawa couplings. Moreover, note that it is natural for λ 5 to be very small since, in case it was exactly zero, a definition of a conserved lepton number would be possible [29].
Most LFV phenomenological studies focus on the radiative decay µ → eγ, ignoring other LFV observables. There are two reasons for this. First, the great performance of the MEG experiment, that recently set the quite impressive bound Br(µ → eγ) < 5.7 × 10 −13 . And second, the dipole dominance in many models of interest. When the dipole contributions originated in photon penguin diagrams dominate, the rate for µ → 3e is correlated with the rate for µ → eγ. In this case a simple relation can be derived [43] Br Since the proportionality factor is much smaller than one, µ → 3e is suppressed with respect to µ → eγ and the latter becomes the process leading to the most stringent constraints. This assumption has been present in all previous works on lepton flavor violation in the scotogenic model [29][30][31][32]. They have either assumed explicitly that photon penguin diagrams dominate or simply ignored 4-fermion observables (like µ → 3e) and concentrated on µ → eγ (an approach consistent with the assumption that the photonic dipole contributions dominate). Here we want to study under what conditions that is a bad simplification of the phenomenology. In order to do so, we consider the ratio 8 In those regions of parameter space where R µe > 1, the observable that provides the most stringent limits is Br(µ → 3e), whereas Br(µ → eγ) would be the most relevant observable in regions where R µe < 1. Since the photonic dipole operators contribute to both observables, the only way to obtain R µe > 1 is to have dominant contributions from box and/or photonic non-dipole diagrams in µ → 3e (Z-penguins are suppressed by charged leptons and thus their contribution is always negligible). Since the photonic non-dipole diagrams, given by the A N D form factor, never exceed the dipole ones as much as to compensate the large factor that multiplies |A D | 2 in the branching ratio formula (see Eq. (22)), they are never dominant. We are therefore left with a competition between photonic dipole operators and box diagrams.
Assuming box dominance in µ → 3e and a degenerate right-handed neutrino spectrum one can estimate R µe ∼ y 4 48π 2 e 2 H(ξ), where y is the average size of the Yukawa coupling and the function H(ξ) is defined as The function H(ξ) is shown in Fig. 5. Notice the cancellation for ξ = 1. This pole is caused by an exact cancellation between the contributions from the loop functions D 1 and D 2 . However, for ξ ≪ 1 and ξ ≫ 1 one always has H(ξ) > 1.
It is clear from Eq. (31) and Fig. 5 that in order to increase the value of R µe one requires large Yukawa couplings and a large mass difference between the right-handed neutrinos and the η scalars (in order to be far from ξ = 1). This is illustrated in Fig. 6, where we show Br(µ → eγ) NH, whereas the right-hand side shows our results for IH. A random Dirac phase δ has been taken.
As can be derived from the spread of the points, this parameter has a much larger influence for NH. As expected from our previous estimate, one can in principle have R µe > 1 (or equivalently, Another parameter that turns out to be very relevant in the determination of the ratio R µe is m ν 1 , the mass of the lightest neutrino. In order to illustrate this fact, we consider two scenarios: TeV. In both cases we assume a degenerate right-handed neutrino spectrum, a random Dirac phase and a random real R matrix. Our numerical results for scenario A are shown in Figs. 7 and 8. The left-hand side of these figures were obtained with NH, whereas the right-hand side shows our results for IH. We see that large values of the lightest neutrino mass may lead to large variations in the Br(µ → eγ) and Br(µ → 3e) branching ratios, and thus in the ratio R µe . We conclude that the LFV rates in the scotogenic model are very sensitive to the absolute scale for neutrino masses. Figure 7 demonstrates that in scenario A the neutrino mass hierarchy also has a clear impact on the LFV rates. While for low m ν 1 , Br(µ → 3e) is clearly below the upper bound for NH, it is largely excluded for IH since it exceeds it. Similarly, while for low m ν 1 the ratio R µe is ∼ 10 −2 in NH (as expected from dipole domination), the contributions from box diagrams already lead to a small increase for IH, where R µe ∼ 0.5.
These figures can be understood by analyzing how the Yukawa couplings depend on m ν 1 . In particular, we must study the combinations of Yukawa couplings that contribute to the LFV processes considered here. Let us suppose that box diagrams dominate ℓ α → 3 ℓ β (otherwise we would be in a dipole dominated scenario where the ratio R µe would not deviate significantly from ∼ 10 −2 ). Then we have the relations Assuming degenerate right-handed neutrinos and a real R matrix, we can use the Casas-Ibarra parametrization in Eq. (11) to obtain Analytical results for the relevant elements of the matrix combinations (or flavor structures) in the previous expressions can be found in appendix B.
Let us first focus on the NH case. Notice that in scenario A we have ξ = (1/4) 2 = 0.0625.
With such a small value for ξ, we expect the D 1 term in the box contribution to dominate over the D 2 term. Therefore, we must inspect the expressions for y † y 21 and y † y 22 or, as shown above, U PMNSmν U †

PMNS 22
. On the one hand, we see that  R µe is then trivially deduced from these considerations. Notice that this quantity can reach values as high as ∼ 50. In this case it is obvious that one cannot ignore Br(µ → 3e), but in fact this branching ratio becomes the most relevant LFV observable.
The discussion for IH would be a bit more involved. In this case we find a larger relevance of the D 2 piece. In fact, for m ν 1 ∼ 10 −2 eV this term competes with the D 1 term, leading to the feature observed on the right-hand sides of Figs. 7 and 8.
Let us now consider our results for scenario B, shown in Figs. 9 and 10. Again, we present our results for NH on the left-hand side and our results for IH on the right-hand side. Regarding NH, it is already clear at first sight that the results are qualitatively very similar to those found in scenario A. Although the LFV rates are very different (much lower in this case), the dependence on m ν 1 is very similar. Notice that all points in these figures are actually allowed by the current limits. This was expected, since it is well-known that LFV constraints are more easily satisfied in scenarios with m N > m η + [29]. On the other hand, the difference between NH and IH found in scenario A is not present in scenario B, in which both cases show the same behavior.
Finally, let us briefly discuss a scenario with non-degenerate right-handed neutrinos. The spectrum in the right-handed neutrino sector has an impact on the LFV rates, as we want to illustrate here. In order to do so, we consider a spectrum of the type m N = (m N ,m N ), with two fixed mass eigenvalues (m (1,2) N ) and one varying (m N ). Although one can imagine other scenarios, this simple family of non-degenerate spectra serves to show the qualitative behavior that we want to emphasize. Fig. 11 shows a representative example of how the LFV rates can change in a non-degenerate right-handed neutrino spectrum. On the left, we show Br(µ → eγ) (blue) and Br(µ → 3e) (red) N = 3 TeV, m η + = 1 TeV and m ν 1 = 10 −3 eV were assumed. A NH spectrum for the light neutrinos was chosen for this figure and we allowed for a random Dirac phase δ.
As naively expected, lowm N values enhance both branching ratios, with Br(µ → 3e) being the one that typically gets the larger enhancements. This is caused by the large box contributions induced by the lightest right-handed neutrino. On the other hand, whenm N ≫m (1,2) N , the contribution of the heaviest right-handed neutrino (with a massm N ) becomes sub-dominant and the LFV rates remain barely the same as in the degenerate case. This implies that the general conclusions drawn from the numerical results shown in this section are not restricted to degenerate scenarios. Besides this fact, we do not find any other remarkable feature in the LFV phenomenology for non-degenerate right-handed neutrinos.

B. Sensitivity to low-energy neutrino parameters
We have already shown the relevant role played by the lightest neutrino mass in the resulting LFV branching ratios. Let us now extend the discussion to the other undetermined low-energy parameter (besides the Majorana phases), the Dirac phase δ.
As starting point, we discuss how a non-zero Dirac phase can change the prediction for Br(ℓ α → ℓ β γ). In order to do that, we consider the ratio Br(ℓ α → ℓ β γ)/Br(ℓ α → ℓ β γ) δ=0 , where Br(ℓ α → ℓ β γ) δ=0 is the value of the branching ratio for δ = 0. This is explicitly shown in Fig. 12, where contours of these ratios are drawn in the m ν 1 − δ plane. In this figure we chose normal hierarchy for the light neutrinos, a degenerate right-handed neutrino spectrum and a real R matrix. Although these results were obtained for specific values of the remaining parameters, we emphasize that the Br(ℓ α → ℓ β γ)/Br(ℓ α → ℓ β γ) δ=0 does not depend on them when the right-handed neutrinos are degenerate and R is a real matrix, see Eqs. (33) and (35).
For the former, we find that the branching ratio can be reduced by almost an order of magnitude, depending on the value of δ. In the latter case, the branching ratio can be increased by a factor of 4 just by switching on the Dirac phase. Moreover, in both cases we find that m ν 1 is also determinant.
We do not show our results for the remaining case, ℓ α = τ and ℓ β = µ, since we found very little dependence on the Dirac phase.
These results tell us that the LFV rates are highly sensitive to the low-energy neutrino parameters. The question then arises as to whether one can get information about them by measuring LFV observables. In case of Br(ℓ α → ℓ β γ), we have already seen that, for the specific scenario of degenerate right-handed neutrinos and a real R matrix, the flavor dependence of the amplitude will be determined just by low energy parameters: neutrino masses, mixing angles and CP violating phases. Therefore, by taking ratios of branching ratios (what we call flavor ratios), the dependence on the high-energy parameters cancels out and we are left with functions of m ν 1 and δ. More precisely, we can make use of Eqs. (33) and (35) to write Note that there is no sum over α, α ′ , β and β ′ in the previous expression.
Our results for these ratios are presented in Fig. 13. We show Br(µ → eγ)/Br(τ → eγ) (to the left) and Br(τ → µγ)/Br(τ → eγ) (to the right) contours in the m ν 1 − δ plane. Under the assumptions of a degenerate right-handed neutrino spectrum and vanishing phases in the R matrix, this figure would allow one to set important constraints on m ν 1 and δ in case two branching ratios To the left for ℓ α = µ, ℓ α ′ = τ and ℓ β = ℓ β ′ = e, to the right for ℓ α = ℓ α ′ = τ , ℓ β = µ and ℓ β ′ = e. Normal hierarchy for the light neutrinos, a degenerate right-handed neutrino spectrum and specific (but generic) values for the free parameters have been assumed, see text for more details.
were measured. Furthermore, we see again the important dependence on these two low-energy parameters, since the ratios can change by more than one order of magnitude.
The same will be true for Br(ℓ α → 3 ℓ β ) when one of the two pieces, D 1 or D 2 , dominates. A particularly interesting scenario arises when the term containing the loop function D 2 ⊂ B gives the dominant contribution. As we have found numerically, this assumption is typically valid for ξ 10 or for large m ν 1 . In this case, the special dependence on the Yukawa matrices, y T y ββ y T y βα , implies that the R matrix drops from the flavor ratios even when it contains complex entries, since We have investigated this scenario and obtained the results in Fig. 14. We concentrate on Br(µ → 3e)/Br(τ → 3e) (on the left) and Br(τ → 3µ)/Br(τ → 3e) (on the right). In the derivation of these plots, we neglected the contribution from the D 1 term. Moreover, we assumed normal hierarchy for the light neutrinos and a degenerate right-handed neutrino spectrum. It is clear that, again, the parameters δ and m ν 1 may have a very strong impact on the 3-body branching ratios.
On the left-hand side of the figure we see that (for this parameter configuration) Br(µ → 3e) is typically larger than Br(τ → 3e). The ratio between these two observables is only close to 1 for δ = π, whereas in the rest of the m ν 1 − δ plane one has Br(µ → 3e) ≫ Br(τ → 3e). On the other hand, the right-hand side of the figure shows that the ratio Br(τ → 3µ)/Br(τ → 3e) is mostly FIG. 14: Br(ℓ α → 3 ℓ β )/Br(ℓ α ′ → 3 ℓ β ′ ) contours in the m ν1 − δ plane. To the left for ℓ α = µ, ℓ α ′ = τ and ℓ β = ℓ β ′ = e, to the right for ℓ α = ℓ α ′ = τ , ℓ β = µ and ℓ β ′ = e. Normal hierarchy for the light neutrinos and a degenerate right-handed neutrino spectrum haven assumed, see text for more details. determined by m ν 1 , with δ playing a secondary role. As for the previous case, the ratio could be close to 1 (for low m ν 1 ) or much larger (for high values of the lightest neutrino mass).
Our study reveals that LFV observables in the scotogenic model are highly sensitive to lowenergy parameters such as the Dirac phase or the lightest neutrino mass. However, it also reveals a large degeneracy, this is, the LFV rates are not correlated with a single parameter. Furthermore, our results regarding flavor ratios have been obtained for a special case: degenerate right-handed neutrinos and real R matrix. In a more general scenario one expects departures from the values of the flavor ratios obtained here. In conclusion, it is not possible to determine the value of a single parameter by measuring a flavor ratio. Only the combination of measurements of the low-energy parameters with the discovery of one (or several) LFV processes can really put the flavor structure of the scotogenic model under experimental test.

C. µ − e conversion in nuclei
So far we have discussed our results on ℓ α → ℓ β γ and ℓ α → 3 ℓ β . Now we move on to discuss µ−e conversion in nuclei. In this model, we have found that Z-penguins give a very little contribution to µ − e conversion in nuclei compared to that of the γ-penguins. In this situation one could naively expect dipole operators to dominate the conversion rate. When this is the case, one expects a simple relation [53] CR(µ − e, Nucleus) where f (Z, N ) is a function that depends on the nucleus and ranges from 1.1 to 2.2 for the nuclei of interest. However, in addition to the dipole contribution given by A D , γ-penguins also have the non-dipole contribution given by A N D . In fact, Eqs. (13) and (16) tell us that, for degenerate right-handed neutrinos, one has A D = 3 F 2 (ξ)/G 2 (ξ)A N D . Therefore, the relative weight of these two different photon contributions depends on the loop functions F 2 (ξ) and G 2 (ξ). These are shown in Fig. 15, where one can see that G 2 (ξ) > F 2 (ξ). For ξ ≪ 1 the difference between G 2 (ξ) and F 2 (ξ) is small and both contributions have similar weights. However, for ξ ≫ 1 (m N ≫ m η + ) one has G 2 (ξ) ≫ F 2 (ξ) and A N D becomes the most relevant contribution 9 . This is illustrated in Fig. 16, where we show our results for Br(µ → eγ) and CR(µ − e, Ti), as well as their ratio.
The same parameter configuration as in Fig. 6 has been selected: fixed values m η + = 1 TeV and m ν 1 = 10 −3 eV, random real R matrix and Dirac phase. These numerical results were obtained for NH, although very similar results are found for IH. We focused on µ − e conversion in titanium, although the same behavior is found for other nuclei.
We find that for large values of ξ, the µ−e conversion rate in titanium gets enhanced by photonic non-dipole contributions. This is a positive result, given the great experimental perspectives for µ − e conversion in nuclei in the near future.

V. SUMMARY AND CONCLUSIONS
The scotogenic model is a popular extension of the standard model that accounts for neutrino masses and dark matter. As for most neutrino mass models, lepton flavor violation is one of the most attractive phenomenological issues, as it may reveal the underlying mechanism that leads to neutrino masses and mixings. In this work we have studied the predictions obtained in the scotogenic model for the LFV processes with the best experimental perspectives in the near future: ℓ α → ℓ β γ, ℓ α → 3 ℓ β and µ − e conversion in nuclei. Full analytical expressions have been derived, going beyond the usual dipole dominance approximation. Our computation includes, besides the dipole photon penguin contribution, non-dipole photon contributions, Z-penguins as well as box diagrams.
The full consideration of all contributions to LFV processes leads to a very interesting picture.
Given the rich LFV phenomenology in the scotogenic model, we are sure that more complete studies can be performed. Here we have explored some of the phenomenological consequences of our analytical results. This may serve as a summary of our main conclusions: • Box diagrams dominate the LFV amplitudes in some parts of parameter space. This scenario leads to a deviation from the naive expectations obtained from the dipole dominance assumption and makes ℓ α → 3 ℓ β more constraining than ℓ α → ℓ β γ.
• The mass hierarchy between the right-handed neutrinos and the inert doublet scalars is of fundamental relevance for LFV observables. We have found that parameter points with large Yukawa couplings and m N ≫ m η + or m N ≪ m η + typically have enhanced box diagrams, thus leading to Br(ℓ α → 3 ℓ β ) > Br(ℓ α → ℓ β γ). This is caused by the particular behavior of the loop functions.
• In the scotogenic model, there are two dark matter candidates: the lightest right-handed neutrino N 1 and the lightest neutral η scalar (η R or η I ) [17]. When ξ > 1, the lightest neutral η constitutes the dark matter of the universe. Otherwise, N 1 is the dark matter particle [29][30][31]52]. In case of N 1 DM (ξ < 1), the only possible annihilation channel is N 1 N 1 → ℓ αlβ , via the Yukawa interaction. For this reason, Yukawa couplings of O(1) are required in order to obtain the observed dark matter relic density Ωh 2 ≈ 0.12 [54], and this may lead to incompatibility with the LFV bounds. It is thus clear that the dark matter phenomenology of N 1 and LFV are closely related. We have explicitly constructed parameter points where all the requirements for right-handed neutrino dark matter are met: m N < m η , large Yukawa couplings and m N in the appropriate range, as found in dedicated studies [52].
Our investigation reveals that although most of these points lead to violation of the LFV bounds, a small fraction of them are perfectly compatible. These valid points involve some small tuning of the parameters and could only be found due to the generality of our scans (not limited to any fixed structure of the Yukawa couplings). These results can be seen as a positive indication in favor of the validity of right-handed neutrino dark matter, although detailed studies are required to get a definitive and robust conclusion. These are, however, beyond the scope of this paper. On the other hand, we would like to point out that in case the dark matter is provided by the scalar η, one can always obtain the correct relic density since, in addition to the Yukawa interactions, this particle has gauge and scalar interactions [55,56], not correlated with LFV.
• The LFV rates are highly sensitive to the low-energy parameters m ν 1 (the mass of the lightest neutrino) and δ (the Dirac phase). In particular, large m ν 1 typically enhances box diagrams.
• In some specific scenarios (with degenerate right-handed neutrinos), the ratios of branching ratios depend only on m ν 1 and δ. Under some assumptions, this may allow us to test the flavor structure of the model.
• Interestingly, the rate for µ − e conversion in nuclei can also be enhanced beyond the dipole contribution in some regions of the parameter space. Our study reveals that non-dipole photon contributions become very relevant for m N ≫ m η + . This may lead to µ−e conversion rates in nuclei as large as the branching ratio for µ → eγ. These are good news given the promising experimental projects in µ − e conversion in nuclei.
We would like to stress that our (qualitative) conclusions are not restricted to Ma's scotogenic model, but should apply to a much wider class of radiative neutrino mass models. In particular, extended versions of the scotogenic model (like the model proposed in [57]) should have, at least in some corners of parameter space, a similar phenomenology.
The presence of TeV scale particles with sizable couplings to the SM states also leads to interesting prospects at the LHC. Although the direct production of the right-handed neutrinos is typically suppressed due to their singlet nature, they will be produced in the decays of the η scalars when this is kinematically allowed. In turn, the η scalars may have non-negligible production crosssections provided they are light. This possibility, not related to the lepton sector, has been studied in some detail. In this case one expects multilepton final states with a significant amount of missing energy [58]. Furthermore, the scotogenic states may also modify the usual Higgs boson decays, with observable implications at the LHC [59,60].
To conclude, the anatomy of lepton flavor violation in the scotogenic model has been fully determined and some interesting phenomenological aspects have been explored. Some definite predictions have been made, and these may be used to put the model under experimental test.
The connection between neutrino masses and lepton flavor violation is a powerful test for this purpose. Hopefully, a positive signal in one (or several) experiments in the next few years will provide valuable hints on the mechanism behind neutrino masses.