A radiative neutrino mass model with SIMP dark matter

We propose the first viable radiative seesaw model, in which the neutrino masses are induced radiatively via the two-loop Feynman diagram involving Strongly Interacting Massive Particles (SIMP). The stability of SIMP dark matter (DM) is ensured by a ℤ5 discrete symmetry, through which the DM annihilation rate is dominated by the 3 → 2 self-annihilating processes. The right amount of thermal relic abundance can be obtained with perturbative couplings in the resonant SIMP scenario, while the astrophysical bounds inferred from the Bullet cluster and spherical halo shapes can be satisfied. We show that SIMP DM is able to maintain kinetic equilibrium with thermal plasma until the freeze-out temperature via the Yukawa interactions associated with neutrino mass generation.


Introduction
The standard model (SM) of particle physics is an enormously successful theory describing the nature of the universe. Nevertheless, the origin of the non-zero neutrino mass [1][2][3][4] and the identification of dark matter (DM) in the universe [5][6][7][8] are the lack of explanations in the SM.
As it is well known, the easiest way to account for tiny neutrino masses is the canonical seesaw mechanism [9][10][11], in which heavy right-handed singlet neutrinos are added to the SM. However, such heavy fermions are very hard to probe by current colliders. Alternatively, people focus on radiative seesaw models [12][13][14][15], where neutrino masses are generated at loop level and the mass scales of the new particles involving in the Feynman diagram can be lighter than the canonical seesaw mechanism.
On the other hand, a number of well-motivated DM candidates have been suggested, the most popular among which is the Weakly Interacting Massive Particles (WIMP) with the mass range spanning from sub-GeV to TeV scale. The WIMP DM is thermally produced in the early universe, and its relic density is usually determined by the strength of the 2 → 2 annihilation cross section of DM into the SM particles. The experimental investigations for the WIMP DM have null results so far, this motivates physicists to come up with the new perspectives for the DM nature. Recently, a novel idea of DM, Strongly Interacting Massive Particles (SIMP) [16] has gotten attention and has been explored in JHEP07(2017)101 the literature . In comparison with WIMP, the relic abundance of SIMP is determined by the strength of the 3 → 2 annihilation cross section of DM into itself, while its mass scale spreads from MeV to sub-GeV, which may be insensitive to present direct searches. The annihilation rate for the 3 → 2 process should be larger than the 2 → 2 annihilation rate to consider SIMP DM instead WIMP. In addition, SIMP DM has to be in kinetic equilibrium with the SM sector until the freeze-out so that the temperature of the dark sector is the same with that in the SM sector, known as the SIMP condition. An advantage of SIMP DM opposite to the WIMP DM is that the SIMP candidate can address some astrophysical issues such as small-scale structure problems [44] and the DM halo separation in Abell 3827 cluster [45,46].
In the economic point of view, any realistic model beyond the SM should incorporate the above crucial ingredients. The most renowned one possessing these necessary components is Ma's scotogenic model [47], in which the WIMP DM is running in the loop diagram to produce the neutrino masses. There are a bunch of studies along this direction [48][49][50][51].
In this article, we propose a brand-new scheme of the scotogenic model, where the role of WIMP DM is replaced by the SIMP DM. To accomplish our thought, we refer to the resonant SIMP model constructed in ref. [24] and extend it by introducing more scalars and fermions for neutrino mass generation. Hereafter, we call it νSIMP model. In this model, the complex scalar is selected as a SIMP DM candidate and is stabilized by a Z 5 symmetry. The resonant effect can reduce the size of the quartic couplings associated with the 3 → 2 annihilation processes so that the perturbative bound and the constraints from the Bullet cluster and spherical halo shapes can be satisfied. The SIMP condition can also be fulfilled via the new Yukawa interactions, which connects the dark sector and the SM sector.
The plan of the paper is as follows. In the next section, we introduce the νSIMP model and give a description of the relevant interactions and masses for the new particles. In section 3, we write down the neutrino mass formula. In section 4, we take into account several experimental and theoretical constraints on the model. In section 5, we evaluate the relic density of the resonant SIMP DM and briefly mention the restrictions from the astrophysical sources. In section 6, we demonstrate the allowed parameter space to make the SIMP condition work. We conclude and summarize our study in section 7. Some lengthy formulas, diagrams, and the benchmark points of the model are put in the appendices.

νSIMP model
To achieve the νSIMP scenario, we add three vector-like fermions, N 1,2,3 , one scalar doublet, η, and two complex singlet scalars, χ and S to the SM, all of which have charges under a conserved Z 5 symmetry, 1 while all of the SM particles are Z 5 neutral. The particle contents and the charge assignments are summarized in table 1. It follows that the lightest mass eigenstate (denoted by X) of the linear combination of χ and the neutral component of η is stable and can serve as a valid SIMP DM candidate. 2 JHEP07(2017)101 Table 1. Charge assignments of the fermions and scalars in the νSIMP model, where E = ν − T is the SM lepton doublet, Φ is the SM Higgs doublet, and ω = exp 2πi/5 is the quintic root of unity.
The renormalizable Lagrangian for the interactions of the scalar particles in this model with one another and with the SM gauge bosons is where D ρ is the SM covariant derivative, and the scalar potential V is with υ 246.22 GeV being the vacuum expectation value (VEV) of Φ. The Hermiticity of V implies that the parameters in the scalar potential µ 2 Φ,η,χ,S , λ Φ,η,χ,S,Φη,Φχ,ΦS,ηχ,ηS,χS , and λ Φη must be real. In the later sections, we will choose µ 1,2 , λ 3 , and κ to be real and assume λ η,Φη,Φχ,ΦS,ηχ,ηS and λ Φη are negligible since these quartic couplings are irrelevant to our numerical analysis.
After spontaneously symmetry breaking, the scalar bosons can be parametrized by with h being the physical Higgs boson. The masses of h, S and η + are then given by The κυ term in the scalar potential causes the mixing between the neutral scalars η 0 and χ. In the basis ( η 0 χ ) T , the corresponding mass matrix is written as (2.5)

JHEP07(2017)101
Upon diagonalizing M 2 ηχ , we get the mass eigenstates H and X and their respective masses m H and m X given by where c ξ = cos ξ, s ξ = sin ξ, and m H > m X . Plugging χ = −s ξ H + c ξ X into eq. (2.2), one can extract the relevant interactions for the 3 → 2 annihilation processes as These couplings manifest the Z 5 discrete symmetry and can produce the 5-point interactions of X by integrating out the complex scalar field S. To generate the neutrino masses, the additional couplings L ⊃ − 1 2 µ 2 s 2 ξ H 2 − 2c ξ s ξ XH S + H.c. are also required. The neutrino masses will be calculated in the next section. From eqs. (2.1), (2.2) and (2.6), the Lagrangian describing the invisible decay channels of the Z boson and the Higgs boson is where g w is the SU(2) gauge coupling constant, and c w = cos θ w with the weak mixing angle θ w . There are also the gauge interactions of the exotic scalars with the photon and the weak bosons, which are related to the electroweak precision tests. We collect them in appendix B. The Lagrangian responsible for the masses and interactions of the vector-like fermions N 1,2,3 is where M k represent the Dirac masses, the summation over j, k = 1, 2, 3 is implicit, the superscript c refers to the charge conjugation, P R,L = 1 2 (1 ± γ 5 ), and 1,2,3 = e, µ, τ . Explicitly, the Yukawa couplings Y rk and Y L,R rk are of the forms as Figure 1. Feynman diagrams for neutrino mass generation at the two-loop level.

Radiative neutrino mass
In the νSIMP model, the neutrinos acquire mass radiatively through two-loop diagrams with internal H, X, S, and N k as shown in figure 1. The resulting neutrino mass matrix defined by L ν = − 1 2 (M ν ) rs ν r ν c s + H.c. is given as [53,54] where the loop functions are , The mass matrix in eq.
The mixing angles in the PMNS matrix and neutrino mass eigenvalues are given by the global fitting to the neutrino oscillation data [55].
As we will discuss in section 6, the order of magnitude of the Yukawa couplings is Y jk ∼ O(0.01-1) with 0.1 GeV M k 1 GeV when the SIMP condition is imposed. In the next section, we will also show that the size of the mixing angle should be s ξ 0.06 due to the constraints from the invisible decays of the Z boson and the Higgs boson. Moreover, in order to satisfy perturbative bounds on the quartic couplings and the observed relic JHEP07(2017)101 density of DM, we find that the cubic coupling µ 2 ∼ O(100 MeV). Accordingly, if one takes Y jk ∼ 0.1, s ξ ∼ 0.05, µ 2 ∼ 100 MeV, Y L,R jk ∼ 0.1, and C L,R jk ∼ 1, the correct neutrino mass scale m ν ∼ 0.1 eV can be arrived. To make this model more reliable, we display the benchmark points in appendix C.

Constraints
There are various experimental and theoretical restrictions on the masses and couplings of the new particles in the νSIMP scenario. Experimentally, the flavor-changing radiative decay r → s γ constrains the Yukawa couplings Y rk . The Feynman diagram depicted such decay process is shown in figure 2. The branching fraction of the decay process is where the fine structure constant α, the Fermi constant G F and the loop function F (z) are withê the electromagnetic charge. The most stringent experimental limit on the µ → eγ process comes from the MEG collaboration [56]. The up-to-date upper bound on its branch- 3 the Yukawa couplings are limited in Y rk 0.02 which is in conflict with the range mentioned in the previous section. The simplest solution to evade this severe constraint is to assume a diagonal Yukawa matrix Y. In this solution, the pattern of neutrino mixing is pinned down by the structures of the other Yukawa matrices Y L,R and the mass hierarchy of the vector-like fermions. 4 At one-loop level, the presence of η ± and N k also induces a contribution to the anomalous magnetic moment a j of charged lepton j given by In particular, the current experimental value for the muon anomalous magnetic moment has more than 3σ deviation from the SM prediction: a exp µ −a SM µ = (288±80)×10 −11 [58]. Since the new contribution given by eq. (4.3) is negative, we then require |∆a µ | < 8 × 10 −10 , this gives an upper bound for the Yukawa coupling Y rk . For example, by taking m η + ∼ 300 GeV and m η + M k , the Yukawa couplings are limited in Y rk O(1) which is less stringent compared to the constraints from the flavor-changing radiative decay. Figure 2. Feynman diagram of flavor-changing radiative decay r → s γ.

JHEP07(2017)101
In our study, we suggest that the lightest complex scalar X is the SIMP DM candidate. Since the mass scale of SIMP DM is MeV to sub-GeV, there is then a new physics contribution to the invisible decay width of the Z boson and the Higgs boson, namely  .8), these upper limits consequently are translated into |s ξ | 0.4 and |s ξ | 0.165 (m H /100 GeV) −1 , respectively. It turns out that the constraint from the Higgs invisible decay width is much stronger than the Z boson one. For instance, by choosing m H ∼ 300 GeV, we then reach the upper bound |s ξ | 0.06.
The scalar masses are constrained by the oblique parameters due to their modifications to the SM gauge boson propagators [63]. From eq. (B.1) and figure 8 in appendix B, those parameters are calculated as where the loop functions are given by (4.5) Since we are interested in the scale that the mass m X is below electroweak scale, one may think that more general definitions of the oblique parameters may be used [64]. However, we have checked the difference is not important because the most stringent constraint comes from the T -parameter whose definition does not change even for below electroweak scale.

JHEP07(2017)101
The current constraints are given in refs. [65,66] as ∆S = 0.05 ± 0.11, ∆T = 0.09 ± 0.13, ∆U = 0.01 ± 0.11 with the correlation coefficients 0.90 (between ∆S and ∆T ), −0.59 (between ∆S and ∆U ), and −0.83 (between ∆T and ∆U ). These limits imply that the heavier neutral component H and the charged component η + should be nearly degenerate m H ≈ m η + in the case of s ξ 1 and m X m H , m η + . Theoretically, the quartic parameters λ j are subject to the conditions of vacuum stability and perturbativity. To ensure the vacuum to be stabilized at large field values, we demand [24] ξ ≈ λ χ , and λ XS ≡ λ χS c 2 ξ ≈ λ χS because of the smallness of the mixing angle ξ. In the previous work [24], a condition of perturbativity on the quartic couplings has been taken, which corresponds to λ X,S < 16π in our convention in eq. (2.2). However, this upper bound seems to be overly optimistic when the RG running is considered. Instead, we force the relatively conserved conditions λ X,S < 4π in our numerical work. Furthermore, since X plays the role of DM, it should not develop the VEV. The sufficient conditions to guarantee X = 0 (as well as S = 0) are given by here we have assumed λ 3 = 0 for simplicity. 5

Resonant SIMP DM and relic abundance
In order to estimate the thermal relic abundance of SIMP DM, we have to solve the Boltzmann equation of the DM number density n DM = n X + nX = 2n X (we assume there is no asymmetry between particles X andX) as follows with H being the Hubble parameter, n eq DM the DM number density at the chemical equilibrium, and σ 3→2 υ 2 rel ≡ 1 24 σ XXX→XX υ 2 rel the thermal averaged effective 3 → 2 annihilation cross section. 6 By applying the standard derivation [68], the approximate solution to the Boltzmann equation for the current relic density Ω DM is given by

2)
5 One can find the necessary conditions to ensure X = 0 by using the method in the literature [67].
However, the analytical result is too long to read. 6 The definition of the effective 3 → 2 annihilation cross section depends on the model. For example, in the Z3 SIMP model [23], σ3→2υ 2 rel ≡ 1 24 σ XXX→XX υ 2 rel + 1 8 σ XXX→XX υ 2 rel . where x = m X /T ,ĥ denotes the normalized Hubble constant, g ,f is the number of relativistic degrees of freedom at the freeze-out temperature, T f = m X /x f , m pl = 1.22 × 10 19 GeV is the Planck mass. To evaluate the thermal average of the 3 → 2 annihilation cross section, we employ the formula in ref. [26] σ 3→2 υ 2 rel =

JHEP07(2017)101
where β = 1 2 υ 2 1 + υ 2 2 + υ 2 3 with υ i the velocities of three initial DM particles. In the νSIMP model, the Feynman diagrams of the 3 → 2 process XXX →XX are shown in figure 3. From eq. (2.7), the effective 3 → 2 annihilation cross section under CP invariance is calculated as whereŝ = (p 1 + p 2 + p 3 ) 2 9m 2 X 1 + 2β/3 and the momenta of DM are neglected except around the resonanceŝ ≈ m 2 S . By taking the mass spectrum 2M k > m S > 2m X , the decay width of the particle S is computed as To enhance the 3 → 2 annihilation cross section, we pick the resonant pole m S √ŝ 3m X in eq. (5.4), and it is convenient to adjust the resonant behavior by defining the following dimensionless parameters as  where S indicates the degeneracy between m S and 3m X , and γ S is the width of the resonance. 7 With these variables, the 3 → 2 annihilation cross section can be expressed in the Breit-Wigner resonant form similar to the one in ref. [69] where the coefficient c X is with R 1,2 = µ 1,2 /m X . Utilizing eq. (5.2), we present the plots of the DM relic density Ω DM versus m S with different values of m X and R 1,2 in figure 4, where the solid lines (light dashed lines) are the predicted values by using the thermal (non-thermal) averaged effective 3 → 2 annihilation cross section. The orange region is the latest relic density data Ω DMĥ 2 = 0.1197 ± 0.0022 given by the Planck collaboration [70]. In these plots, we do not vary the mixing angle ξ since dependence of the mixing angle is extremely small as long as s ξ 1. Also, in order to examine the conditions of X = 0 in eq. (4.7) easily, we again assume λ 3 = 0. For nonzero λ 3 , the numerical results are similar as pointed out in ref. [24]. We have checked our choices of the values of R 1,2 can accommodate the requirements of perturbativity, R 2 1 < λ S < 4π and R 2 2 /9 λ X < 4π with m S 3m X . According to the plots, one can see that the low values of R 1,2 are disfavored if the DM mass m X is heavier. Besides fitting the relic abundance of DM, there are the other astrophysical observations from the Bullet cluster [71][72][73] and spherical halo shapes [74], which impose the bound σ self /m X 1 cm 2 /g with σ self = 1 4 (σ XX→XX + σ XX→XX + σXX →XX ) the effective self-interacting cross section. We depict in figure 5 the Feynman diagrams of the DM self-interacting processes in our SIMP model, and their cross sections are calculated as here we have neglected the contributions from the h and Z-mediated diagrams due to their small couplings and mass suppression. By choosing an appropriate value of λ X R 2 2 m 2 X /m 2 S < λ X < 4π , the bounds from the Bullet cluster and spherical halo shapes can be satisfied. For instance, if we take m X (m S ) = 30 (93) MeV, R 1,2 = 2, 5 and ξ = 0.05 with λ X = 7, we find σ self /m X 0.26 cm 2 /g. More examples and discussions can be found in ref. [24].

SIMP condition
In the SIMP paradigm, DM is thermally produced through the 3 → 2 annihilation process into the particles in the dark sector rather than the 2 → 2 annihilation process into the SM particles. On the other hand, in order to keep the temperature of the dark sector as the same with the SM sector, the SIMP candidate needs to be in kinetic equilibrium with the SM sector. Thus the naive criteria that DM can be a SIMP candidate is given by [16] Γ 2→2 < Γ 3→2 < Γ kin , (6.1)

JHEP07(2017)101
which should be held during the freeze-out temperature. In this inequality, each reaction rate is defined by Γ 2→2 = n X σ 2→2 υ rel , Γ 3→2 = n 2 X σ 3→2 υ 2 rel , and Γ kin = n SM σ kin υ rel , 8 where the number densities of DM and the SM particles are given as [39] with g counts the internal degrees of freedom, m SM being the mass of the SM particle, (+) applies to fermions, and (−) pertains to bosons. It is also pointed out in ref. [37] that a slightly stronger SIMP condition may be derived by considering the rate of energy transfer (rather than the rate of reaction) between the SM and dark sectors. This rigorous SIMP condition can be written as |Ė 3→2 | < |Ė kin |, whereĖ 3→2 is the rate of the DM mass transferring to the kinetic energy in the DM bath, andĖ kin is the rate of the kinetic energy of DM transfers to the thermal plasma. Quoting the detailed calculations of the SIMP condition in ref. [34], we then impose Γ 3→2 < Γ 3→2 < 10 −2 Γ kin in our numerical study. From eq. (2.9), the particle X can interact with the active neutrinos via the Yukawa couplings Y rk . 9 The Feynman diagrams of the elastic scattering between X and the SM neutrinos are displayed in figure 6(a), and its reaction rate can be computed as where the neutrino number density n ν and the thermally averaged effective scattering cross section are given by withζ(3) 1.202 the Riemann zeta function of 3. By the crossing symmetry, the Feynman diagrams for the 2 → 2 annihilation process of a DM pair into a pair of the SM neutrino are shown in figure 6(b), 10 and the reaction rate is calculated as Γ 2→2 = n X r,s σ XX→νrνs υ rel , (6.5) 8 In the WIMP paradigm, the Boltzmann equation of the DM number density is given bẏ where σ2→2υ rel is the thermal averaged effective 2 → 2 annihilation cross section. Due to this definition, an extra factor 1/2 is multiplied to the DM cross sections (eq. (6.4) and (6.6)). This comes from the fact that DM and anti-DM particles are not identical in our case [69]. 9 The particle S can also have the 3 → 2 annihilation processes, but it can only interact with the SM sector through the Higgs portal. In this case, the reaction rate would be governed by the Higgs mass, and may be too small to keep kinetic equilibrium with the SM sector. 10 Here we have neglected the scattering processes X ± → X ± and the annihilation channels XX → + − due to the mass suppression of the Z boson and the Higgs boson. Figure 6. (a) Feynman diagrams of the elastic scattering between X and the SM neutrinos.

JHEP07(2017)101
(b) Feynman diagrams for the 2 → 2 annihilation process of a DM pair into a pair of the SM neutrino.  where the thermally averaged effective 2 → 2 annihilation cross section is given by with n X given by eq. (6.2). For simplification of numerical treatment, here we assume the masses of the vector-like fermions are degenerate (M 1 = M 2 = M 3 = M ). The reaction rates of the 2 → 2 annihilation process and the elastic scattering are then reduced to the form as

JHEP07(2017)101 7 Conclusion
In this paper, we have built a SIMP version of the scotogenic model, where the SIMP DM has the responsibility to generate the neutrino masses and its stability is guaranteed by the Z 5 discrete symmetry. We have considered the experimental and theoretical constraints on the masses and the couplings in the model including the neutrino masses and mixings, lepton flavor violating processes, anomalous magnetic moment, the invisible decay modes of the Z boson and the Higgs boson, the electroweak precision data, perturbativity of the couplings and vacuum stability. In the models of SIMP DM, a large coupling is generally required in order to reproduce the correct DM relic abundance measured by experiments through 3 → 2 annihilating processes. This may give a tension with perturbativity and potential stability. By employing the resonant mechanism in our model, the correct relic abundance of DM has been reproduced, and the bounds on the quartic couplings and the self-scattering cross section have been fulfilled at the same time. We found the parameter space of the new Yukawa interactions such that the SIMP condition is achieved. Since our model faces to the stringent constraints from the Higgs invisible decay and the direct search of new charged scalars, it will be tested in near future. A Gauged U(1) B−L extension of the νSIMP model It is believed that there is no global symmetry can exist in a theory of quantum gravity [75,76]. Under this context, the discrete symmetry we introduced in our νSIMP model may originate from a gauge symmetry (gauge redundancy). At certain energy scale, this gauge symmetry is broken down to the Z 5 discrete symmetry by a nonzero VEV of a scalar field. In the following, we demonstrate an extension of the νSIMP model to the gauged U(1) B−L version by adding one more SM singlet complex scalar ζ. The particle contents and the charge assignments are summarized in table 2. In this extended model, the Lagrangian associated with the 3 → 2 processes is given by After spontaneous symmetry breaking, the complex scalar ζ can be expanded around its VEV as ζ = 1  Table 2. Charge assignments of the particles in the gauged U(1) B−L extension of the νSIMP model.
which corresponding to the first three terms in the last line of eq. (2.2), respectively, with µ 1,2 = λ 1,2 υ . The Yukawa couplings contributed to the neutrino mass diagrams are the same with in eq. (2.9), and the lightest scalar particle involving in the diagrams can be a SIMP DM candidate. On the other hand, the SIMP condition can be achieved by Z -portal instead of the Yukawa interactions due to the new gauge boson in this model. We leave the detailed study of the model to future work.

B Gauge interactions
The kinetic part of the Lagrangian in eq. (2.  for m X = 40 MeV, m S = 128 MeV, µ 2 = 320 MeV. One can check that eq. (C.1) and (C.2) satisfy the SIMP condition (Γ 3→2 /Γ 2→2 10 4 and Γ kin /Γ 3→2 10 3 ) as shown in the left and right panels of figure 7, respectively. These benchmark points give normal ordering neutrino mass eigenvalues and mixing angles consistent with neutrino oscillation data. It is also possible to take benchmark parameter sets in the cases for Y L Y R and inverted hierarchy, though these are not shown here.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.