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 $\mathbb{Z}_5$ discrete symmetry, through which the DM annihilation rate is dominated by the $3 \to 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.


I. 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 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 Sec. III, we write down the neutrino mass formula. In Sec. IV, we take into account several experimental and theoretical constraints on the model. In Sec. V, we evaluate the relic density of the resonant SIMP DM and briefly mention the restrictions from the astrophysical sources. In Sec. VI, we demonstrate the allowed parameter space to make the SIMP condition work. We conclude and summarize our study in Sec. VII. Some lengthy formulas, diagrams, and the benchmark points of the model are put in the appendices.

II. ν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 Tab. I. 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 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. 1 This discrete symmetry can be realized as a remnant of the U(1) gauge symmetry as discussed in Ref. [24,52]. A concrete example is given in Appendix A. 2 In the simplest Z 3 SIMP model [16], the quartic coupling in the scalar potential is too large to satisfy the bound from perturbativity.
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 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), 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.(1), (2) and (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 where Y j k = Y jk .

III. 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 Fig. 1. The resulting neutrino mass matrix defined by where the loop functions are with The mass matrix in Eq. (11) is diagonalized by the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) The mixing angles in the PMNS matrix and neutrino mass eigenvalues are given by the global fitting to the neutrino oscillation data [55].
Feynman diagrams for neutrino mass generation at the two-loop level.
As we will discuss in Sec. VI, 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 density of DM, we find that the cubic coupling 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.

IV. 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 Fig. 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 branching ratio is 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 3 If η ± decays dominantly into electron or muon, m η + 270 GeV is required in order to avoid the constraint from the left-handed slepton search [57]. 4 If m τ > m e,µ + 2M k , the new decay modes τ → (e, µ)N N → (e, µ)νν X X open and would contribute to τ → (e, µ) + missing energy. However, this constraint is not so stringent.
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. (16) 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.
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 Γ new Z, h → inv = Γ(Z, h → XX). The present experimental bounds on these invisible decay widths are Γ new Z → inv < 2 MeV (at the 95% C.L.) [ (6) and (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.(B1) and Fig. 8 in Appendix B, those parameters are calculated as where the loop functions are given by 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. The current constraints are given in Ref. [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] where λ X ≡ λ χ c 4 ξ ≈ λ χ , 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). 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

V. 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 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 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 Fig. 3. From Eq.(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 6 The definition of the effective 3 → 2 annihilation cross section depends on the model. For example, in the Z 3 SIMP model [23], σ 3→2 υ 2 rel ≡ 1 24 σ XXX→XX υ 2 rel + 1 8 σ XXX→XX υ 2 rel .
particle S is computed as To enhance the 3 → 2 annihilation cross section, we pick the resonant pole m S √ŝ 3m X in Eq. (24), 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. (22), we present the plots of the DM relic density Ω DM versus m S with different values of m X and R 1,2 in Fig. 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.(20) 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 Fig. 5 the Feynman diagrams of the DM self-interacting processes in our SIMP model, 7 From Eqs. (20), (25) and (26) with ξ 1 and m S 3m X , one can easily show that γ S 10 −3 R 2 2 10 −2 λ X . Thus one obtains γ S 0.1 1 with the perturbative bound λ X < 4π.  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].

VI. 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] 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.(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 Fig. 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 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. (33) and (35)). 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.
where Y ≡ r,s,k,l Re Y * rk Y rl Y sk Y * sl 1/4 . Using the SIMP condition at T f , we illustrate the plots of the magnitude of |Y| as a function of M in Fig. 7 with different numerical inputs based on

VII. 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. In this extended model, the Lagrangian associated with the 3 → 2 processes is given by L ζ = 1 √ 2 λ 1 ζ * χ * S 2 + 1 √ 2 λ 2 ζ * χ 2 S + 1 6 λ 3 χ 3 S * + H.c. .

(A1)
After spontaneous symmetry breaking, the complex scalar ζ can be expanded around its VEV as ζ = 1 √ 2 ς + υ , where υ ≡ √ 2 ζ . The scalar interactions between χ and S are then extracted as L ζ ⊃ 1 2 λ 1 υ χ * S 2 + 1 2 λ 2 υ χ 2 S + 1 6 λ 3 χ 3 S * + H.c. , which corresponding to the first three terms in the last line of Eq.(2), respectively, with µ 1,2 = λ 1,2 υ . The Yukawa couplings contributed to the neutrino mass diagrams are the same with in Eq. (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.

�� � �
for m X = 40 MeV, m S = 128 MeV, µ 2 = 320 MeV. One can check that Eq. (C1) and (C2) 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 Fig. 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.