Multicomponent Dark Matter in Radiative Seesaw Models

We discuss radiative seesaw models, in which an exact $Z_2\times Z_2'$ symmetry is imposed. Due to the exact $Z_2\times Z_2'$ symmetry, neutrino masses are generated at a two-loop level and at least two extra stable electrically neutral particles are predicted. We consider two models: one has a multi-component dark matter system and the other one has a dark radiation in addition to a dark matter. In the multi-component dark matter system, non-standard dark matter annihilation processes exist. We find that they play important roles in determining the relic abundance and also responsible for the monochromatic neutrino lines resulting from the dark matter annihilation process. In the model with the dark radiation, the structure of the Yukawa coupling is considerably constrained and gives an interesting relationship among cosmology, lepton flavor violating decay of the charged leptons and the decay of the inert Higgs bosons.


I. INTRODUCTION
Neutrino oscillation experiments show that the neutrinos have tiny masses and mix with each other. It is a clear evidence for physics beyond the standard model (SM), since the SM has no mechanism for giving masses to the neutrinos. The global fit [1] shows that the mass-squared differences of the neutrinos are ∆m 2 21 = 7.50 +0.19 −0.17 × 10 −5 eV 2 and ∆m 2 31 = 2.524 +0.039 −0.040 (−2.514 +0.038 −0.041 )×10 −3 eV 2 for normal (inverted) mass hierarchy. The cosmological data, on the other hand, gives the upper bound of the sum of the neutrino masses as Σ j m ν j < 0.23 eV [2], a scale twelve orders of magnitudes smaller than the electroweak scale. It is one of the most important problems of particle physics to reveal the origin of the tiny masses for the neutrinos.
Type-I seesaw mechanism [3] is one of the attractive way to realize the tiny masses of the neutrinos, where the right-handed neutrinos are introduced to the SM. If the neutrino Yukawa coupling for the Dirac neutrino mass is O(1), the mass of the right-handed neutrino has to be around O (10 15 ) GeV to obtain eV-scale neutrinos. The mass scale of O(10 15 ) GeV is obviously beyond the reach of collider experiments. Even for the mass of the righthanded neutrinos around O(1) TeV, the direct search of the right-handed neutrinos would be difficult because of the tiny neutrino Yukawa couplings of O (10 −6 ).
Another attractive way to give the neutrino masses is a radiative generation (the so-called radiative seesaw model). The original idea of radiatively generating neutrino masses due to TeV-scale physics has been proposed by Zee [4], in which the neutrino masses are induced at the one-loop level because of the addition of an isospin doublet scalar field and a charged singlet field to the SM. Another possibility for generating neutrino masses via the new scalar particles is e.g. the , in which the neutrino masses arise at the two-loop level.
A further extension with a TeV-scale right-handed neutrino has been proposed in Ref. [6]. In this model the neutrino masses are induced at the three-loop level, where the Dirac neutrino mass at the tree level is forbidden due to an exact Z 2 symmetry. The right-handed neutrino is odd under the Z 2 and becomes a candidate of dark matter (DM). The idea of simultaneous explanation for the neutrino masses via the radiative seesaw mechanism and the stability of DM by introducing an exact discrete symmetry has been discussed in many models (see, e.g., Refs. [7][8][9][10][11][12][13][14] and the recent review [15] and references therein).
The model proposed by Ma in Ref. [7] is one of the simplest radiative seesaw model with DM candidates. The model has the Z 2 -odd right-handed neutrinos N k and the inert doublet scalar field η = (η + , η 0 R + iη 0 I ) T . The neutrino masses are generated at the one-loop level, in which the Yukawa interactions Y ν ik L i ηN k and the scalar interaction (λ 5 /2)(H † η) 2 contribute to the neutrino mass generation. The mass matrix is expressed as where M k is the Majorana mass of the k-th generation of right-handed neutrino, m η 0 R and m η 0 I are the mass of the η 0 R and η 0 I , respectively. In this model, we have two scenarios with respect to the DM candidate; the lightest right-handed neutrino N 1 or the lighter Z 2 -odd neutral scalar field (η 0 R or η 0 I ). The phenomenology of the model is studied in Refs. [16][17][18][19][20][21]. A DM candidate can be made stable by an unbroken symmetry. The simplest possibility of such a symmetry is Z 2 symmetry as in the above models. However, if the DM stabilizing symmetry is larger than Z 2 : Z N (N ≥ 4) or a product of two or more Z 2 s, the DM is consisting of stable multi-DM particles (multicomponent DM system). A supersymmetric extension of the radiative seesaw model of Ref. [7] is an example [14]. Other possibilities with multicomponent DM are widely discussed in [9][10][11][12][13][22][23][24].
In this paper we study two models of the two-loop extension of the model by Ma [7], we call them as "model A" and "model B", in which due to the Z 2 ×Z ′ 2 symmetry a set of stable particles can exist. Introducing two new scalar fields, the λ 5 term is generated radiatively in the model A [12,13]. In this model we discuss the three component DM system in which the two new scalar fields and a right-handed neutrino are the DM candidate. Such case has been discussed in [13], however, we reanalyze the model since the benchmark points studied in [13], where the masses of both new scalars are several hundred GeV, has been excluded by the recent results of the direct detection DM experiments. In this paper we focus on the case where the mass of one of the scalar DMs is close to the half of the Higgs boson mass to satisfy the constraints from the direct detection. In the model B the right-handed neutrinos have the mass radiatively generated through the one loop of internal new fermion and scalar fields. We identify the lightest right-handed neutrino as dark radiation.
We start in section II by writing down a set of the Boltzmann equations of the multicomponent DM system. The model A is discussed in section III by following Ref. [13]. In section IV we discuss the model B and relate dark radiation with the lepton flavor violating decay of the charged leptons and the decay of the inert Higgs bosons. Summery and discussion are given in section V.

II. MULTICOMPONENT DARK MATTER SYSTEMS
In the case of one-component DM the relic density of DM χ is determined by the Boltzmann equationṅ where n χ is the DM number density,n χ is the equilibrium number density and σ χχ→XX ′ v is the thermally averaged cross section for χχ → XX ′ . Here X and X ′ stand for the SM particles. The Hubble parameter H is given by H = 1.66×g 1/2 * T 2 /M PL , where g * is the total number of effective degrees of freedom, T and M PL are the temperature and the Planck mass, respectively. It is convenient to rewrite the equation in terms of dimensionless quantities; 1. Standard annihilation (left), DM conversion (middle) and semi-annihilation (right).
the number per comoving volume Y χ = n χ /s and the inverse temperature x = m/T . Here s is the entropy density s = (2π 2 /45)g * T 3 and m is the mass of the DM particle. Using the replacement of dx/dt = H| T =m /x, we obtain The thermally averaged cross section σ χχ→XX ′ v of O(10 −9 ) GeV with a DM mass of 100 GeV gives Y χ ≃ 10 −12 , which is consistent with the observed value of the relic abundance Ωh 2 ≃ 0.12 [25].
In the multicomponent DM system three types of processes enter in the Boltzmann equations 1 : See Fig. 1 for a depiction of three types of processes. Here we assume that none of the DM particles have the same quantum number with respect to the DM stabilizing symmetry. The Boltzmann equations for the DM particle χ i with mass m i are Here x = µ/T and 1/µ = ( i m −1 i ) is the reduced mass of the system. The contributions of non-standard annihilations have been discussed in e.g. [24] for two and three component DM system with a Z 2 × Z ′ 2 symmetry.
1 Semi-annihilation processes also exist in one-component DM systems when DM is a Z 3 charged particle [26] or a vector boson [27].

III. MODEL A
In the following, by extending the one-loop model in [7] we study two of the two-loop radiative seesaw models with Z 2 ×Z ′ 2 symmetry. We refer to them as "model A" and "model B". Owing to the Z 2 × Z ′ 2 symmetry, there exist at least two extra stable electrically neutral particles. The multicomponent DM system is realized in the model A, while one of the stable particles plays as the dark radiation in the model B.
The matter content of the model A is shown in Table I. In addition to the matter content of the SM model, we introduce the right-handed neutrino N k , an SU(2) L doublet scalar η, and two SM singlet scalars χ and φ. Note that the lepton number L ′ of N is zero. The Z 2 × Z ′ 2 × L ′ -invariant Yukawa sector and Majorana mass term for N can be described by where i, j, k (= 1, 2, 3) stand for the flavor indices. The scalar potential V is written as The Z 2 × Z ′ 2 is the unbroken discrete symmetry while the lepton number L ′ is softly broken by the last term in the potential V m , the φ mass tem. In the absence of this term, there will be no neutrino mass. Note that the "λ 5 term", (λ 5 /2)(H † η) 2 , is also forbidden by L ′ . A small λ 5 of the original model of Ma [7] is "natural" according to 't Hooft [28], because the absence of λ 5 implies an enhancement of symmetry. In fact, if λ 5 is small at some scale, it remains small for other scales as one can explicitly verify [17]. Here we attempt to derive the smallness of λ 5 dynamically, such that the λ 5 term becomes calculable.
The Higgs doublet field H, the inert doublet field η and the singlet scalar φ are respectively parameterized as where v h is the vacuum expectation value. The tree-level masses of the scalars are given by Althogh the tree-level mass of η 0 R is the same as that of η 0 I as shown in (13), the degeneracy is lifted at the one-loop level via the effective λ 5 term: This correction is embedded into the two-loop diagram to generate the neutrino mass (see Fig. 2). The 3×3 neutrino mass matrix M ν can be given by where we have assumed that Using λ eff 5 given in (17), the neutrino mass matrix can be approximated as We see from (17) and (19) that the neutrino mass matrix M ν is proportional to |Y ν κ| 2 m 2 5 . When m χ , m φ R , m η 0 , M k ∼ O(10 2 ) GeV, for instance, implies that |Y ν κ|m 5 ∼ O(10 −1 ) GeV to obtain the neutrino mass scale of O(0.1) eV. With the same set of the parameter values we find that λ eff 5 ∼ 10 −6 , where the smallness λ eff 5 is a consequence of the radiative generation of this coupling. As we will see, the product |Y ν κ| enters into the semi-annihilation of DM particles which produces monochromatic neutrinos, while the upper bound of |Y ν | follows from the µ → eγ constraint.

A. Multicomponent dark matter system
In the model A there are three type of dark matter candidates ; N 1 (the lightest among there are two candidates. In the following discussions we assume that N 1 is a DM candidate [13]. The other possibility, η 0 R -DM, is discussed in [12].
We discuss the three DM system of N 1 , φ R , χ. There are three types of DM annihilation process: Here we have assumed m φ R > m χ and m φ R + m χ < M 2,3 . Moreover, since the mass difference between φ R and φ I is controlled by the lepton-number breaking mass m 5 , which is assumed to be much smaller than m φ R . Then m φ R and m φ I are practically degenerate and the contribution of φ I to the annihilation processes during the decoupling of DMs is nonnegligible. The diagrams for annihilation processes which enter the Boltzmann equation are shown in Figs. 3-5. Since the reaction rate of the conversion between φ R and φ I can reach chemical equilibrium during the decoupling of DMs, we can sum up the number densities of φ R and φ I and compute the relic abundance of Ω φ R h 2 [13].
In the multicomponent DM scenario, the effective cross section off the nucleon is given by In our model, only χ and φ R scatter with the nucleus, and the right-handed neutrino N 1 does not interact with nucleus at tree level. So we can neglect the N 1 contribution at the lowest order in perturbation theory. The effective cross sections of φ R and χ are expressed as where σ χ and σ φ R are the spin independent cross sections and given by Heref ∼ 0.3 is the usual nucleonic matrix element [29] and m N is nucleon mass. The upper bounds on the cross section off the nucleon is obtained by LUX [30] and XENON1T [31]. In the cases of one-component DM system of a real or complex scalar boson, those experimental results give the strong constraint on the masses of those DM particles; the allowed DM mass region is ≃ m h /2 and > ∼ 1 TeV [32]. In the model A with the multicomponent DM system, the constrains on the cross sections off the nucleon for χ and φ R are also relatively severe. As a benchmark we take the mass of χ as m χ = m h /2 while vary the mass of φ R in the following analysis 2 .
In the original one-loop neutrino mass model in [7], the relic density of N 1 tends to be larger than the observational value [16]. The additional contributions coming from the semiannihilation can enhance the annihilation rate for N 1 so that the N 1 DM contribution to 2 Two singlet scalar DM scenario in Z 2 × Z ′ 2 model has been explored in detail in Ref. [22]. the total relic abundance can be suppressed. This situation is realized for M 1 > m φ R , m χ as can be seen later.
As the benchmark set we take the following values for the parameters.
The masses of heavier right-handed neutrinos are M 2 = M 3 = 1 TeV. The mass differences between m η 0 R and the sum of m χ and m φ R are so chosen that no resonance appears in the schannel of the semi-annihilation in Fig. 5 (right). The benchmark set satisfies the constraints from the perturbativeness, the stability conditions of the scalar potential [12,13], the lepton flavor violation (LFV) such as µ → eγ [33] and the electroweak precision measurements [34,35]. It is noted that κ is bounded as |κ| 0.4 by the perturbativeness and the stability conditions [12,13]. Figure 6 shows the relic abundances of Ω χ h 2 , Ω φ R h 2 and Ω N 1 h 2 and the total relic abundance Ω total h 2 (= Ω χ h 2 + Ω φ R h 2 + Ω N 1 h 2 ) as a function of m φ R for the benchmark set. The horizontal dashed line stands for the observed value Ω obs h 2 ∼ 0.12. It is shown that the relic abundance of the χ is Ω χ ≃ Ω obs /2. When φ R is lighter than N 1 , the semi-annihilation tends to decrease the relic abundance of N 1 . For the benchmark set, the total relic abundance is consistent with the observed value around m φ R ≃ 280 GeV.
The left panel in Fig. 7 shows the contour plot for the m φ R -γ 5 plane where the total relic density of DM can be made consistent with the observed value Ω obs h 2 ∼ 0.12. We take two values, 10 GeV (black line) and 1 GeV (red line), for the mass difference between m η 0 R and m χ + m φ R in (28). The other parameters are taken as the same in Eqs. (27)-(31). We can see the scalar coupling γ 5 increases drastically as m φ R increases for m φ R > ∼ 290 GeV. It is because the relic density of the N 1 DM, Ω N 1 h 2 , becomes significant for m φ R > ∼ 290 GeV, so that Ω φ R h 2 should be drastically suppressed. The scalar couplings of DM particles with the SM Higgs boson, γ 2 and γ 5 , and the DM masses are constrained by the DM direct detection experiments. For the χ DM, the effective cross section off nucleon σ eff χ in Eq. (24) is σ eff χ ∼ 10 −47 cm 2 for the benchmark set. It is an order of magnitude smaller of the spin independent cross section off the nucleon given by LUX [30] and XENON1T [31], respectively.
The hatched region is excluded by perturbativity. In both panels, we take two values, 10 GeV (black line) and 1 GeV (red line), for the mass difference between m η 0 R and m χ + m φR .
than the current experimental bound. For the φ R DM, the right panel in Fig. 7 shows the relation between m φ R and the effective cross section σ eff  hatched region is excluded by perturbativity. Although the scalar coupling γ 5 becomes large for m φ R > ∼ 290 GeV and then the cross sections off the nucleon σ φ R becomes large, the effective cross section σ eff φ R decreases for m φ R > ∼ M 1 (= 300 GeV), since the abundance of φ R decreases. For the case of (m χ + m φ R ) − m η 0 R =10 GeV, it can be seen that the mass region 288 GeV < ∼ m φ R is excluded by LUX and XENON1T data. On the other hand, there are no constraints from the direct DM search experiments on the mass of φ R for the case of (m χ + m φ R ) − m η 0 R =1 GeV. This is because the relic abundance of φ R becomes much smaller by the large contribution from the s-channel process of the semi-annihilation. Figure 8 shows the same as in Fig. 7 but for M 1 = 500 GeV and γ 7 = 0.28. From the right panel in Fig. 8, we see that 485 (490) (1) GeV. Before we go to discuss indirect detection, we summarize the parameter space, in which a correct value of the total relic DM abundance Ω total h 2 can be obtained without contradicting the constraint from the direct detection experiments. As in the case of the single SM singlet DM, the constraint is in fact very severe: The mass of χ has to be very close to m h /2, and γ 2 (the Higgs portal coupling) also has to be close to 0.004 for an adequate amount of Ω χ . However, as for m φ R and γ 5 , there exist a certain allowed region. The allowed region in the m φ R -γ 5 plane is controlled by the semi-annihilation (especially, the last diagram in Fig. 5, which is sensitive to the mass relation (28)) and the DM conversion (especially the right diagram in Fig. 4, which is sensitive to γ 7 ). If we increase the mass of the righthanded neutrino DM, the mass of φ R increases, but how the allowed range in the m φ R -γ 5 plane emerges remains the same. If we take the larger γ 7 , e.g. γ 7 = 0.28, in Fig. 7, the allowed region for m φ R becomes narrower as 295 GeV < ∼ m φ R < ∼ 300 GeV. The smaller m φ R ( < ∼ 295 GeV) is excluded by Ω total < Ω obs due to the larger DM conversion i.e. the larger annihilation process of φ R φ R → χχ → XX.

B. Indirect detection
For indirect detections of DM the SM particles produced by the annihilation of DM are searched. Because the semi-annihilation produces a SM particle, this process can serve for an indirect detection. In our model, especially, the SM particle from the semi-annihilation process as shown in Fig. 5 is neutrino which has a monochromatic energy spectrum. Therefore, we consider below the neutrino flux from the Sun [ [36][37][38][39][40] as a possibility to detect the semi-annihilation process of DMs.
The DM particles are captured in the Sun losing their kinematic energy through scattering with the nucleus. Then captured DM particles annihilate each other. The time dependence of the number of DM n i in the Sun is given by [37][38][39][40] where i, j, k = χ, φ R , N 1 and C i is the capture rates in the Sun: and C A 's are the annihilation rate: Here f (m i ) depends on the form factor of the nucleus, elemental abundance, kinematic suppression of the capture rate, etc., varying O(0.01−1) depending on the DM mass [39,40]. V ij is an effective volume of the Sun with µ ij = 2m i m j /(m i +m j ) in the non-relativistic limit. In the Eq. (32) we have neglected the DM production processes such as jj → ii and jk → iX because the kinetic energy of the produced particle i is much larger than that corresponding to the escape velocity from the Sun, i.e. ∼ 10 3 km/s [37,41]. Consequently, the number of the right-hand neutrino DM cannot increase in the Sun, and hence the semi-annihilation process, φ R χ → N 1 ν, is the only neutrino production process 3 , where its reaction rate as a function of t is given by Γ(φ R χ → Nν; t) = C A (φ R χ → N 1 ν)n φ R (t)n χ (t). Figure 9 shows the m φ R dependence of the neutrino production rate today Γ(φ R χ → Nν; t 0 ), where t 0 = 1.45×10 17 s is the age of the Sun, for the same parameter space as in Fig. 7 3 There are also neutrinos having continuous energy spectrum from the decay of SM particles produced by the standard annihilation. The upper bounds for the production rates of the SM particles are given in [41,42]. ( Fig. 9 (left)) and in Fig. 8 (Fig. 9 (right)). The hatched region is excluded by perturbativity. Arrows indicate the excluded regions by the direct detection experiments. For m φ R > ∼ M 1 where the relic abundance of N 1 dominates that of φ R , the neutrino production rate decreases since the capture rate of the φ R becomes small. As we can see from The upper limits on the DM DM → XX ′ from the Sun are given by IceCube experiment [42]. For instance, the upper limit on the annihilation rate of the DM of 250 (500) GeV into W + W − is 1.13 × 10 21 (2.04 × 10 20 ) s −1 and that into τ + τ − is 5.99 × 10 20 (7.96 × 10 19 ) s −1 , which is at least 10 2 times larger than the rate Γ(ν) shown in Fig. 9. Note however that the energy spectrum of the neutrino flux produced by the W or τ decay is different from the monochromatic neutrino. With an increasing resolution of energy and angle the chance for the observation of the semi-annihilation and hence of a multicomponent nature of DM can increase.

IV. MODEL B
Neutrinos have always played consequential roles in cosmology (see [43], and also [44] and references therein). While they play a role as hot dark matter, the mechanism of their mass generation is directly connected to cosmological problems such as baryon asymmetry of the Universe [45] and dark matter [ [6][7][8][9][12][13][14][15][16]. Resent cosmological observations with increasing accuracy [25,[46][47][48] provide useful hints on how to extend the neutrino sector. Here we propose an extension of the neutrino sector such that the tensions among resent different cosmological observations can be alleviated. The tensions have emerged since the first Planck result [25] in the Hubble constant H 0 and in the density variance σ 8 in spheres of radius 8h −1 Mpc: The Planck values of 1/H 0 and σ 8 are slightly larger than those obtained from the observations of the local Universe such as Cepheid variables [47] and the Canada-France-Hawaii Telescope Lensing Survey [49], respectively. The Planck galaxy cluster counts [50] and also the Sloan Digital Sky Survey data [46] yield a smaller σ 8 .
It has been recently suggested [50][51][52] that these tensions can be alleviated if the number N eff of the relativistic species in the young Universe is slightly larger than the standard value 3.046 and the mass of the extra relativistic specie is of O(0.1) eV [52]. Here we suggest a radiative generation mechanism of the neutrino mass, which is directly connected to the existence of a stable DM particle and also a non-zero ∆N eff = N eff − 3.046.
The matter content of the model is shown in Table II. It is a slight modification of the model A: χ in this model is a Majorana fermion. The Z 2 × Z ′ 2 × L ′ -invariant Yukawa sector (the quark sector is suppressed) is described by the Lagrangian where i, j, k = 1, 2, 3, and we have assumed without loss of generality that the χ mass term is diagonal. We also assume that Y e ij have only diagonal elements. The most general renormalizable form of the Z 2 × Z ′ 2 × L ′ -invariant scalar potential is given by and the mass terms are where the m 4 term in (39) breaks L ′ softly. The scalar fields H, η and φ are defined in (10). Since we assume that the discrete symmetry Z 2 × Z ′ 2 is unbroken, the scalar fields above do not mix with other, so that their tree-level masses can be simply expressed: The two-loop diagram for the neutrino mass is shown in Fig. 10. Because of the soft breaking of the dimension two operator φ 2 , the propagator between φ and φ can exist, generating the mass of N: The 3×3 two-loop neutrino mass matrix M ν is given by We can also use (45) to obtain an approximate formula for the neutrino mass where U N is the unitary matrix diagonalizing the mass matrix (M N ) ij with the eigenvalues M k and the mass eigenstates N ′ k , and we have used the fact that M k ≪ m η 0 R ≃ m η 0 I . In the following discussions we choose the theory parameters so as to be consistent with the global fit

A. Dark radiation
According to the discussion at the beginning of this section, we identify the lightest righthanded neutrino with dark radiation contributing to ∆N eff 4 . Without los of generality we may assume it is N ′ 1 with mass < ∼ 0.24 eV. The upper bound on the mass is obtained together with 3.10 < N eff < 3.42 in Ref. [52]. To simplify the situation, we require that the heavier right-handed neutrinos N ′ 2 and N ′ 3 decay above the decoupling temperature T dec N of N ′ 1 . Their decay widths are given by where we have used m η 0 = m η 0 R ≃ m η 0 I and neglected the mass of N ′ 1 and ν L s. Therefore, is satisfied, where H(T ) is the Hubble constant at temperature T , and g * s (T ) is the relativistic degrees of freedom at T . To obtain the effective number of the light relativistic species N eff [44,53], we have to compute the energy density of N ′ 1 at the time of the photon decoupling, where we denote the decoupling temperature of γ, ν L and N ′ 1 by T γ0 , T dec ν and T dec N , respectively. Further, T ν0 (T N 0 ) stands for the temperature of ν L (N ′ 1 ) at the decoupling of γ. The most important fact is that the entropy per comoving volume is conserved, so that sa 3 is constant, where s is the entropy density, and a is the scale factor. The effective number N eff follows from ρ r (T γ0 ) = (π 2 /15)(1 + (7/8)(4/11) 4/3 N eff ) T 4 γ0 and is given by [44,53,54] N eff = 3.046 + g * s (T dec ν ) g * s (T dec N ) 4/3 (51) for N ν = 3, where ρ r is the energy density of relativistic species. Since g * s (T dec ν ) = (11/2) + (7/4)N ν = 10.75, we need to compute the decoupling temperature T dec N to obtain g * s (T dec N ) and hence N eff . For 0.05 < ∼ ∆N eff < ∼ 0.38 [52] we find 101 > ∼ g * s (T dec N ) > ∼ 22 and also T dec N ≃ 165 MeV to obtain g * s (T dec N ) ≃ 30 (which gives ∆N eff = 0.25). To estimate T dec N , we compute the annihilation rate Γ N (T ) of N ′ 1 at T , which is given by where ζ(3) ≃ 1.202 . . . and n N (T ) is the number density of N ′ 1 . Then we calculate T dec N from Γ N (T dec N ) = H(T dec N ), which can be rewritten as 5 with (Y ν ) 2 = i |Y ′ν i1 | 2 . It turns out that M 2,3 ∼ O(10) GeV to obtain ∆N eff ∼ 0.25 while satisfying M 1 < ∼ 0.24 eV. To see this, we first find that which follows from (53) for ∆N eff ∼ 0.25. Further we can estimate a part of (49) from the neutrino mass (47) with M ν ∼ 0.05 eV: where we have used m 2 Then using (50) with T ≃ 165 MeV (which corresponds to ∆N eff ≃ 0.25), we obtain Note that this is an order of magnitude estimate, and indeed M 2 (3) can not be smaller than 10 GeV to satisfy ∆N eff < ∼ 0.38.
Since we require that M 1 < ∼ 0.24 eV, there exists a huge hierarchy in the right-handed neutrino mass. This has a consequence on the Yukawa coupling matrix Y ′ν : To obtain realistic neutrino masses with the mixing parameters given in (48), has to be satisfied. Note that only |Y ′ν i1 | enters into the thermally averaged annihilation cross section of N ′ 1 , as we can see from (52). Because of ∆N eff < ∼ 0.38, on the other hand, |Y ′ν i1 | can not be made arbitrarily large. The hierarchy (57) has effects on the LFV radiative decays of the type l i → l j γ, so that the LFV decays and ∆N eff are related, as we will see below. In the limit m j ≪ m i , where m i and m j stand for the mass of l i and l j , respectively, the ratio of the partial decay widthB(l i → l j γ) = Γ(l i → l j γ)/Γ(l i → ν i eν e ) can be written as [18] Here m η ± and Y ′ν ik are defined in (40) and (47), respectively, and the current upper bounds on the branching fraction of these processes [33,56] require From (59) we find that Y ′ν 31 is not constrained by the stringent constraint from µ → eγ, which will be crucial in obtaining a realistic N eff without having any contradiction with (59)-(61). Furthermore, if Y ′ν 31 is large compared with others and the hierarchy (57) is satisfied, the ratio R =B(τ → µγ)B(τ → eγ)/B(µ → eγ) is ∼ |Y ′ν 31 | 2 , and from the same reason ∆N eff depends mostly on Y ′ν 31 . A benchmark set of the input parameters is given by which yields sin 2 θ 12 = 0.305, sin 2 θ 23 = 0.441, sin 2 θ 13 = 0.0213, where we have assumed that Y ′ν ij are all real so that there is no CP phase. These values are consistent with (48), (59) -(61). With the same input parameters we find: The lhs of (50) =5.46×10 −21 (1.78×10 −20 ) GeV for N 2 (N 3 ), where the rhs is H = 2.10 × 10 −20 GeV with T dec N = 166.8 MeV and g * s (T dec N ) = 30.83, and ∆N eff = 0.245. In Fig. 11 we plot R 1/2 against ∆N eff with m η ± = 240 GeV and m η 0 R = 220 GeV, where we have varied m η 0 I between 221 and 227 GeV. In the black region of Fig. 11 the differences of the neutrino mass squared and the neutrino mixing angles are consistent with (48) for the normal hierarchy, and the constraints M 1 < 0.24 eV, (50) and (59)-(61) are satisfied. If ∆N eff and R 1/2 would depend on Y ′ν 31 only, we would obtain a line in the ∆N eff -R 1/2 plane. The Y ′ν 11 and Y ′ν 21 dependence in R 1/2 cancels, but this is not the case for ∆N eff . This is the reason why we have an area instead of a line in Fig. 11. We see from Fig. 11 that the predicted region for ∆N eff < ∼ 0.1 is absent. The main reason is that we have assumed that M 2 , M 3 < ∼ 16 GeV. This has also a consequence on the difference between m 2 η 0 R and m 2 η 0 I , because the mass difference changes the overall scale of the neutrino mass (47). To obtain a larger M 2,3 , we can decrease the mass difference, thereby implying an increase of the degree of fine-tuning. Further, the difference between m 2 η 0 R and m 2 η 0 I can not be made arbitrarily large, because it requires a smaller M 2,3 , which due to H(T ) ∝ T 2 in turn implies that the decoupling temperature T dec N has to decrease to satisfy the constraint (50). A smaller T dec N , on the other hand, means a larger ∆N eff which is constrained to be below 0.38. This is why m η 0 I is varied only in a small interval in Fig. 11.
Since the current upper bound on B(µ → eγ) ≃B(µ → eγ) is 4.2 × 10 −13 [33], the model B predicts which is about two orders of magnitude smaller than the current experimental bounds [56]. Another consequence of the hierarchy (57) is that the total decay width of η R depends on i,j |Y ′ ij | 2 , which is approximately i |Y ′ i1 | 2 (we assume that η R is the lightest among ηs). Therefore, ∆N eff is basically a function of the decay width. In Fig. 12 we show ∆N eff against Γ η R /m η 0 R , the decay width of η 0 R over m η 0 R , where we have used the same parameters as for Fig. 11. η 0 R decays almost 100 percent into neutrinos and dark radiation N ′ 1 , which is invisible. In contrast to this, η + can decay into a charged lepton and N ′ 1 , and the decay width over m η ± is the same as Γ η R /m η 0 R . Γ η R should be compared with the decay width for η + → W + * η 0 R,I → ff ′ N ′ 1 ν, which is ∼ 10 −8 m η ± for the same parameter space as for Fig. 12, where f and f ′ stand for the SM fermions (except the top quark). Therefore, η + decays almost 100 percent into a charged lepton and missing energy. In Ref. [20], a similar hierarchical spectrum of the right-handed neutrinos in the model of [7] has been assumed (the lightest one has been regarded as a warm dark matter) and collider physics has been discussed. How the inert Higgs bosons can be produced via s-channel exchange of a virtual photon and Z boson [57] is the same, but the decay of the inert Higgs bosons is different because of the hierarchy (57) of the Yukawa coupling constants. As it is mentioned above, the η ± decays in the present model almost only into the lightest one N ′ 1 and a charged lepton. Therefore, the cascade decay of the heavier right-handed neutrinos into charged leptons will not be seen at collider experiments, because they can be produced only as a decay product of η ± . The decay width of η ± into an individual charged lepton depends of course on the value of Y ′ν i1 . In the parameter space we have scanned we cannot make any definite conclusion on the difference.

B. Cold dark matter and its direct and indirect detection
Since the lightest N is dark radiation and the masses of the heavier ones are O(10) GeV (as we have seen in the previous subsection), η 0 R,I can not be DM candidates, because they decay into N and ν. So, DM candidates are χ and the lightest component of φ 6 . In the case that ηs are lighter than φ R and the lightest component of φ (which is assumed to be φ R ) is DM, a correct relic abundance Ω φ R h 2 = 0.1204 ± 0.0027 [25] can be easily obtained, because γ 3 for the scalar coupling (η † η)|φ| 2 is an unconstrained parameter so far. So, in the following discussion we assume that φ R is DM. Because of the Higgs portal coupling γ 2 , the direct detection of φ R is possible. The current experimental bound of XENON1T [31] of the spin-independent cross section σ SI off the nucleon requires |γ 2 | < ∼ 0.05 ∼ 0.14 for m φ R = 250 ∼ 500 GeV. Since γ 2 is allowed only below an upper bound (which depends on the DM mass m φ R ), γ 3 can vary in a certain interval for a given DM mass.
With this remark, we note that the capture rate of DM in the Sun is proportional to σ SI , while its annihilation rate in the Sun is proportional to the thermally averaged annihilation cross section, vσ( [36][37][38][39][40]. If a pair of φ R s annihilates into η 0 R η 0 R and also η 0 I η 0 I , a pair of ν L andν L will be produced, which may be observed on the Earth [41]. The signals will look very similar to those coming from W ± , which result from DM annihilation. The annihilation rate as a function of time t is given by [40] where C φ R is the capture rate in the Sun, , against σ SI for m φR = 250 (black) and 500 (red) GeV, where m η is fixed at 230 GeV (all ηs have the same mass) and 0.117 < Ω φR h 2 < 0.123. The black (red) vertical dashed line is the XENON1T [31] upper bound on σ SI for m φR = 250 (black) and 500 (red) GeV. and C A is given by We have used f (250 GeV) ≃ 0.5 and f (500 GeV) ≃ 0.2 [40], and we have assumed that all the ηs have the same mass and therefore C A (η 0 η 0 ) = C A (η + η − ). In Fig. 13 we plot the annihilation rate Γ(φ R φ R → η 0 η 0 ; t 0 ) today (t 0 = 1.45×10 17 s) against σ SI for m φ R = 250 and 500 GeV with m η fixed at 230 GeV and 0.117 < Ω φ R h 2 < 0.123. The vertical dashed lines are the XENON1T upper bound on σ SI [31]. The peak of Γ(φ R φ R → η 0 η 0 ; t 0 ) for m φ R = 250 (500) GeV appears at σ SI = 4.2 (4.7) × 10 −46 cm 2 and is ≃ 1.7 (0.7) × 10 18 s −1 , which is two to three orders of magnitude smaller than the upper bound on the DM annihilation rate into W ± in the Sun [42] .

V. CONCLUSION
We have discussed the extensions of the Ma model by imposing a larger unbroken symmetry Z 2 × Z ′ 2 . Thanks to the symmetry, at least two stable particles exit. We have studied two models, model A and model B, where the stable particles form a multicomponent DM system in the model A, while they are a DM and dark radiation in the model B.
The model A is an extension of the model of Ma such that the lepton-number violating "λ 5 coupling", which is O(10 −6 ) to obtain small neutrino masses for Y ν ∼ 0.01, is radiatively generated. Consequently, the neutrino masses are generated at the two-loop level, where the unbroken Z 2 × Z ′ 2 symmetry acts to forbid the generation of the one-loop mass. Such larger unbroken symmetry implies that the model involves a multi-component DM system. We have considered the case of the three-component DM system: two of them are SM singlet real scalars and the other one is a right-handed neutrino. The DM conversion and semiannihilation in addition to the standard annihilation are relevant to the DM annihilation processes. We have found that the non-standard processes have a considerable influence on the DM relic abundance. We also have discussed the monochromatic neutrinos from the Sun as the indirect signal of the semi-annihilation of the DM particles. In the cases of onecomponent DM system of a real scalar boson or of a Majorana fermion the monochromatic neutrino production by the DM annihilation is strongly suppressed due to the chirality of the left-handed neutrino. However, such suppression is absent when DM is a complex scalar boson or a Dirac fermion. Also in a multicomponent DM system, the neutrino production is unsuppressed if it is an allowed process. We have found that the rate for the monochromatic neutrino production in the model A is very small compared with the current IceCUBE [42] sensitivity. However, the resonant effect in the s-channel process of the semi-annihilation can be expected to enhance the rate.
In the model B, the mass of the right-handed neutrinos are produced at the one-loop level. Then the radiative seesaw mechanism works at the two-loop level. Thanks to Z 2 × Z ′ 2 there exist at least two stable DM particles; a dark radiation N ′ 1 with a mass of O(1) eV and the other one, DM, is the real part of φ. The dark radiation contributes to ∆N eff < 1 such that the tensions in cosmology that exist among the observations in the local Universe (CMB temperature fluctuations and primordial gravitational fluctuations) can be alleviated. Because of the hierarchy M 2,3 ≫ T dec N ≃ O(100) MeV ≫ M N 1 O(1) eV, we are able to relate to the ratio of the lepton flavor violating decays to ∆N eff . The indirect signal of the neutrino from the Sun has also been discussed. It is found that the predicted annihilation rate of the neutrinos is two to three orders of magnitude smaller than the current bound [42]. We have also expressed ∆N eff as a function of the decay width of η 0 R (which is assumed to be lightest among ηs). It decays 100 percent into left-and right-handed neutrinos, where the heavier right-handed neutrinos decay further into dark radiation (the lightest among them). Dark radiation appears as a missing energy in collider experiments. We also have found that η + decays 100 percent into a charged lepton and the missing energy. This is a good example in which, through the generation mechanism of the neutrino masses, cosmology and collider physics are closely related.