A Radiative Neutrino Mass Model with SIMP Dark Matter

TL;DR

This paper introduces a novel radiative neutrino mass model where SIMP dark matter is stabilized by a discrete symmetry, with its relic abundance and astrophysical constraints successfully addressed.

Contribution

It presents the first viable radiative seesaw model incorporating SIMP dark matter with a $ ext{Z}_5$ symmetry and detailed phenomenological analysis.

Findings

Neutrino masses generated via two-loop diagrams involving SIMP DM.

SIMP DM relic abundance achieved with perturbative couplings in the resonant scenario.

Astrophysical bounds from Bullet cluster and halo shapes satisfied.

Abstract

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.

Figures

Tables

Table 1. Table 1: Charge assignments of the fermions and scalars in the ν 𝜈 \nu SIMP model, where E = ( ν ℓ − ) T E=\big{(}\,\nu~{}\,\,\ell^{-}\big{)}{}^{\text{T}} is the SM lepton doublet, Φ Φ \Phi is the SM Higgs doublet, and ω = exp ⁡ ( 2 ​ π ​ i / 5 ) 𝜔 2 𝜋 𝑖 5 \omega\,=\,\exp\big{(}2\pi i/5\big{)} is the quintic root of unity.

	$E$	$\Phi$	$N_{1,2,3}$	$\eta$	$\chi$	$S$
SU(2)	$\mathrm{𝟐}$	$\mathrm{𝟐}$	$\mathrm{𝟏}$	$\mathrm{𝟐}$	$\mathrm{𝟏}$	$\mathrm{𝟏}$
U(1)${}_Y$	$-$1/2	1/2	0	1/2	0	0
${\mathbb{Z}}_5$	1	1	${\omega }^2$	${\omega }^2$	${\omega }^2$	$\omega$

Table 2. Table 2: Charge assignments of the particles in the gauged U ​ ( 1 ) B − L U subscript 1 B L \text{U}(1)_{\text{B}-\text{L}} extension of the ν 𝜈 \nu SIMP model.

	$E$	$\Phi$	$N_{1,2,3}$	$\eta$	$\chi$	$S$	$\zeta$
SU(2)	$\mathrm{𝟐}$	$\mathrm{𝟐}$	$\mathrm{𝟏}$	$\mathrm{𝟐}$	$\mathrm{𝟏}$	$\mathrm{𝟏}$	$\mathrm{𝟏}$
U(1)${}_Y$	$-$1/2	1/2	0	1/2	0	0	0
U(1)${}_{B-L}$	$-$1	0	$-$3/5	2/5	2/5	6/5	2

Equations

$$\mathcal{L}=({\mathcal{D}}^{\rho }\Phi {)}^{†}{\mathcal{D}}_{\rho }\Phi +({\mathcal{D}}^{\rho }\eta {)}^{†}{\mathcal{D}}_{\rho }\eta +{\partial }^{\rho }{\chi }^{\ast }{\partial }_{\rho }\chi +{\partial }^{\rho }{S}^{\ast }{\partial }_{\rho }S-\mathcal{V},$$

$$\mathcal{V}$$

$$\displaystyle\Phi\,=\,\begin{pmatrix}0\\[2.84544pt] \frac{1}{\sqrt{2}}\big{(}h+\upsilon\big{)}\end{pmatrix}~{},\quad\eta\,=\,\begin{pmatrix}\,\eta^{+}\,\\[2.84544pt] \,\eta^{0}\,\end{pmatrix}~{},$$

$$m_h^2=\frac{1}{2}{\lambda }_{\Phi }{\upsilon }^2,m_S^2={\mu }_S^2+\frac{1}{2}{\lambda }_{\Phi S}{\upsilon }^2,m_{{\eta }^{+}}^2={\mu }_{\eta }^2+\frac{1}{2}{\lambda }_{\Phi \eta }{\upsilon }^2.$$

$$\displaystyle M^{2}_{\eta\chi}\,\equiv\,\begin{pmatrix}m_{\eta}^{2}&m_{\eta\chi}^{2}\,\\[4.26773pt] \,\,m_{\eta\chi}^{2}&m_{\chi}^{2}\,\end{pmatrix}\,=\,\begin{pmatrix}\,\,\mu_{\eta}^{2}+\frac{1}{2}\big{(}\lambda_{\Phi\eta}+\lambda_{\Phi\eta}^{\prime}\big{)}\upsilon^{2}&\frac{1}{2}\kappa\upsilon^{2}\,\\[4.26773pt] \frac{1}{2}\kappa\upsilon^{2}&\mu_{\chi}^{2}+\frac{1}{2}\lambda_{\Phi\chi}\upsilon^{2}\,\end{pmatrix}~{}.$$

$$\left(\begin{array}{c}{\eta }^0\\ \chi \end{array}\right)$$

$$2m_{H,X}^2$$

$$\displaystyle{\cal L}\,\,\supset\,\,-\,\tfrac{1}{2}\mu_{1}c_{\xi}\Big{[}X^{\ast}S^{2}+X\big{(}S^{\ast}\big{)}^{2}\,\Big{]}-\tfrac{1}{2}\mu_{2}c_{\xi}^{2}\Big{[}X^{2}S+\big{(}X^{\ast}\big{)}^{2}S^{\ast}\Big{]}-\tfrac{1}{6}\lambda_{3}c_{\xi}^{3}\Big{[}X^{3}S^{\ast}+\big{(}X^{\ast}\big{)}^{3}S\Big{]}~{}.$$

$$\displaystyle\,{\cal L}\,\,\supset\,\,\frac{ig_{\text{w}}s_{\xi}^{2}}{2c_{\text{w}}}\scalebox{1.2}{\big{(}}X^{\ast}\partial^{\rho}X-X\partial^{\rho}X^{\ast}\scalebox{1.2}{\big{)}}Z_{\rho}\,-\scalebox{1.2}{\big{(}}\lambda_{\Phi X}|X|^{2}+\lambda_{\Phi S}|S|^{2}\scalebox{1.2}{\big{)}}\upsilon h~{},$$

$$\displaystyle\lambda_{\Phi X}\equiv\,\lambda_{\Phi\chi}c_{\xi}^{2}\,+\kappa c_{\xi}s_{\xi}+\scalebox{1.2}{\big{(}}\lambda_{\Phi\eta}\,+\lambda_{\Phi\eta}^{\prime}\scalebox{1.2}{\big{)}}s_{\xi}^{2}~{},$$

$${\mathcal{L}}_N$$

$$\displaystyle{\cal Y}\,=\,\left(\begin{array}[]{ccc}Y_{e1}&Y_{e2}&Y_{e3}\\ Y_{\mu 1}&Y_{\mu 2}&Y_{\mu 3}\\ Y_{\tau 1}&Y_{\tau 2}&Y_{\tau 3}\end{array}\right)~{},\quad{\cal Y}^{L,R}\,=\,\left(\begin{array}[]{ccc}{\cal Y}^{L,R}_{11}&{\cal Y}^{L,R}_{12}&{\cal Y}^{L,R}_{13}\\[4.26773pt] {\cal Y}^{L,R}_{21}&{\cal Y}^{L,R}_{22}&{\cal Y}^{L,R}_{23}\\[4.26773pt] {\cal Y}^{L,R}_{31}&{\cal Y}^{L,R}_{32}&{\cal Y}^{L,R}_{33}\end{array}\right)~{},$$

$$\displaystyle\big{(}{\cal M}_{\nu}\big{)}_{rs}\,=\,\frac{\mu_{2}{\cal Y}_{rj}{\cal Y}_{sk}s_{2\xi}^{2}}{4(4\pi)^{4}}\Big{(}{\cal Y}_{jk}^{L}{\cal C}_{jk}^{L}+{\cal Y}_{jk}^{R}{\cal C}_{jk}^{R}\Big{)}~{},$$

$${\mathcal{C}}_{jk}^L$$

$${\mathcal{C}}_{jk}^R$$

$${\mathcal{I}}_L(a,b)=\frac{a^2\ln a}{(1-a)(a-b)}+\frac{b^2\ln b}{(1-b)(b-a)},{\mathcal{I}}_R(a,b)=\frac{a\ln a}{(1-a)(a-b)}+\frac{b\ln b}{(1-b)(b-a)},$$

$$m_{XjS}^2=\frac{\hat{u}{m}_X^2+\hat{v}{M}_j^2+\hat{w}{m}_S^2}{\hat{w}(1-\hat{w})},m_{HjS}^2=\frac{\hat{u}{m}_H^2+\hat{v}{M}_j^2+\hat{w}{m}_S^2}{\hat{w}(1-\hat{w})}.$$

$$\displaystyle{\cal B}\big{(}\ell_{r}\rightarrow\ell_{s}\gamma\big{)}\,=\,\frac{3\alpha\,{\cal B}\big{(}\ell_{r}\rightarrow\ell_{s}\nu_{r}\nu_{s}\big{)}}{64\pi G_{\text{F}}^{2}m_{\eta^{+}}^{4}}\,\Bigg{|}\sum_{k\,=\,1}^{3}{\cal Y}_{rk}^{\ast}{\cal Y}_{sk}\,{\cal F}\bigg{(}\frac{M_{k}^{2}}{m_{\eta^{+}}^{2}}\bigg{)}\Bigg{|}^{2}~{},$$

$$\alpha =\frac{{\hat{e}}^2}{4\pi },G_{\text{F}}=\frac{1}{\sqrt{2}{\upsilon }^2},\mathcal{F}(z)=\frac{1-6z+3z^2+2z^3-6z^2\ln z}{6(1-z{)}^4},$$

$$\displaystyle\Delta a_{\ell_{j}}=\,-\frac{m_{\ell_{j}}^{2}}{16\pi^{2}m_{\eta^{+}}^{2}}\sum_{k\,=\,1}^{3}\big{|}{\cal Y}_{jk}\big{|}^{2}{\cal F}\bigg{(}\frac{M_{k}^{2}}{m_{\eta^{+}}^{2}}\bigg{)}~{}.$$

$$\Delta \mathcal{S}$$

$$\Delta \mathcal{T}$$

$$\Delta \mathcal{U}$$

$$F(a,b)$$

$$G(a,b)$$

$$\displaystyle\lambda_{X,S}>0~{},~{}~{}\lambda_{XS}>-\tfrac{1}{2}\sqrt{\lambda_{X}\lambda_{S}}~{},~{}~{}|\lambda_{3}|<\,\sqrt{\frac{18\lambda_{X}\lambda_{S}\lambda_{XS}-8\lambda_{XS}^{3}+\scalebox{1.1}{\big{(}}4\lambda_{XS}^{2}+3\lambda_{X}\lambda_{S}\scalebox{1.1}{\big{)}}^{3/2}}{3\lambda_{X}}}~{},$$

$${\lambda }_X>\frac{{\mu }_2^2}{m_S^2},{\lambda }_S>\frac{{\mu }_1^2}{m_X^2},{\lambda }_{XS}>0,$$

$$\frac{d{n}_{\text{DM}}}{dt}+3\mathcal{H}{n}_{\text{DM}}$$

$$\displaystyle\Omega_{\text{DM}}\hat{h}^{2}\,\simeq\,\frac{\,\,5.7\times 10^{8}~{}\text{GeV}^{-1}}{m_{X}g_{\star,f}^{3/4}m_{\text{pl}}^{1/2}J^{1/2}}~{},\quad J\,=\,\mathop{{\int_{x_{f}}^{\infty}}}dx\,\frac{\big{\langle}\sigma_{3\to 2}\upsilon_{\text{rel}}^{2}\big{\rangle}}{x^{5}}~{},\quad x_{f}\,\simeq\,20~{},\quad$$

$$\displaystyle\big{\langle}\sigma_{3\to 2}\upsilon_{\text{rel}}^{2}\big{\rangle}$$

$${\sigma }_{3\to 2}{\upsilon }_{\text{rel}}^2$$

$$\displaystyle\Gamma_{S}\,=\,\Gamma\big{(}S\to X\bar{X}\big{)}\,=\,\frac{\mu_{2}^{2}c_{\xi}^{2}}{32\pi m_{S}}\sqrt{1-\frac{4m_{X}^{2}}{m_{S}^{2}}}~{}.$$

Full text

††preprint: KUNS-2675††preprint: TUM-HEP/1082/17

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 Wendell:2010md ; Abdurashitov:2009tn ; Abe:2012tg ; Gando:2010aa and the identification of dark matter (DM) in the universe Begeman:1991iy ; Massey:2007wb ; Harvey:2015hha ; Adam:2015rua 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 Minkowski:1977sc ; Yanagida:1979as ; GellMann:1980vs , 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 Zee:1980ai ; Zee:1985rj ; Zee:1985id ; Babu:1988ki , 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\to 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) Hochberg:2014dra has gotten attention and has been explored in the literature Acharya:2017szw ; Bernal:2015bla ; Bernal:2015lbl ; Bernal:2015ova ; Bernal:2015xba ; Bernal:2017mqb ; Choi:2015bya ; Choi:2016tkj ; Choi:2016hid ; Choi:2017mkk ; Cline:2017tka ; Dey:2016qgf ; Farina:2016llk ; Forestell:2016qhc ; Halverson:2016nfq ; Hansen:2015yaa ; Hochberg:2014kqa ; Hochberg:2015vrg ; Kamada:2016ois ; Kamada:2017tsq ; Kuflik:2015isi ; Kuflik:2017iqs ; Lee:2015gsa ; Lee:2015uva ; Pappadopulo:2016pkp ; Tsumura:2017knk ; Yamanaka:2014pva . In comparison with WIMP, the relic abundance of SIMP is determined by the strength of the $3\to 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\to 2$ process should be larger than the $2\to 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 Elbert:2014bma and the DM halo separation in Abell 3827 cluster Massey:2015dkw ; Kahlhoefer:2015vua .

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 Ma:2006km , 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 Ma:2007gq ; Krauss:2002px ; Aoki:2008av ; Gustafsson:2012vj .

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. Choi:2016tkj and extend it by introducing more scalars and fermions for neutrino mass generation. Hereafter, we call it $\nu$SIMP model. In this model, the complex scalar is selected as a SIMP DM candidate and is stabilized by a ${\mathbb{Z}_{5}}$ symmetry. The resonant effect can reduce the size of the quartic couplings associated with the $3\to 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 $\nu$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 $\nu$SIMP model

To achieve the $\nu$SIMP scenario, we add three vector-like fermions, $N_{1,2,3}$, one scalar doublet, $\eta$, and two complex singlet scalars, $\chi$ and $S$ to the SM, all of which have charges under a conserved ${\mathbb{Z}_{5}}$ symmetry,111This discrete symmetry can be realized as a remnant of the U$(1)$ gauge symmetry as discussed in Ref. Choi:2016tkj ; Ho:2016aye . A concrete example is given in Appendix A. while all of the SM particles are ${\mathbb{Z}_{5}}$ neutral. The particle contents and the charge assignments are summarized in Tab. 1. It follows that the lightest mass eigenstate (denoted by $X$) of the linear combination of $\chi$ and the neutral component of $\eta$ is stable and can serve as a valid SIMP DM candidate.222In the simplest $\mathbb{Z}_{3}$ SIMP model Hochberg:2014dra , the quartic coupling in the scalar potential is too large to satisfy the bound from perturbativity.

The renormalizable Lagrangian for the interactions of the scalar particles in this model with one another and with the SM gauge bosons is

[TABLE]

where $\,{\cal D}_{\rho}\,$ is the SM covariant derivative, and the scalar potential ${\cal V}$ is

[TABLE]

with $\upsilon\simeq 246.22\,\,\rm{GeV}$ being the vacuum expectation value (VEV) of $\Phi$. The Hermiticity of ${\cal V}$ implies that the parameters in the scalar potential $\mu_{\Phi,\eta,\chi,S}^{2}$, $\lambda_{\Phi,\eta,\chi,S,\Phi\eta,\Phi\chi,\Phi S,\eta\chi,\eta S,\chi S}$, and $\lambda^{\prime}_{\Phi\eta}$ must be real. In the later sections, we will choose $\mu_{1,2},\lambda_{3}$, and $\kappa$ to be real and assume $\lambda_{\eta,\Phi\eta,\Phi\chi,\Phi S,\eta\chi,\eta S}$ and $\lambda_{\Phi\eta}^{\prime}$ are negligible since these quartic couplings are irrelevant to our numerical analysis.

After spontaneously symmetry breaking, the scalar bosons can be parametrized by

[TABLE]

with $h$ being the physical Higgs boson. The masses of $\,h,\,S$ and $\,\eta^{+}$ are then given by

[TABLE]

The $\kappa\upsilon$ term in the scalar potential causes the mixing between the neutral scalars $\eta^{0}$ and $\chi$. In the basis $\scalebox{1.1}{(}\,\eta^{0}~{}\,\chi\,\scalebox{1.1}{)}{}^{\text{T}}$, the corresponding mass matrix is written as

[TABLE]

Upon diagonalizing $M^{2}_{\eta\chi}$, we get the mass eigenstates $\,H$ and $\,X$ and their respective masses $m_{H}$ and $m_{X}$ given by

[TABLE]

where $c_{\xi}=\cos\xi,\,s_{\xi}=\sin\xi$, and $m_{H}>m_{X}$. Plugging $\chi=-\,s_{\xi}H+c_{\xi}X$ into Eq.(2), one can extract the relevant interactions for the 3 $\to$ 2 annihilation processes as

[TABLE]

These couplings manifest the ${\mathbb{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 \,{\cal L}\,\supset-\,\tfrac{1}{2}\mu_{2}\big{(}s_{\xi}^{2}H^{2}-2c_{\xi}s_{\xi}XH\big{)}S+\text{H.c.} are also required. The neutrino masses will be calculated in the next section. From Eqs.(1), (2) and (II), the Lagrangian describing the invisible decay channels of the $Z$ boson and the Higgs boson is

[TABLE]

[TABLE]

where $g_{\text{w}}$ is the $\text{SU}(2)$ gauge coupling constant, and $c_{\text{w}}=\cos\theta_{\text{w}}$ with the weak mixing angle $\theta_{\text{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

[TABLE]

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}=\frac{1}{2}(1\pm\gamma_{5})$, and $\ell_{1,2,3}=e,\mu,\tau$. Explicitly, the Yukawa couplings ${\cal Y}_{rk}$ and ${\cal Y}^{L,R}_{rk}$ are of the forms as

[TABLE]

where $\,Y_{\ell_{j}k}={\cal Y}_{jk}$.

III Radiative Neutrino mass

In the $\nu$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 ${\cal L}_{\nu}=-\frac{1}{2}({\cal M}_{\nu})_{rs}\overline{\nu_{r}}\nu_{s}^{\text{c}}+\text{H.c.}$ is given as Aoki:2014cja ; Ding:2016wbd

[TABLE]

where the loop functions are

[TABLE]

with

[TABLE]

[TABLE]

The mass matrix in Eq. (17) is diagonalized by the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix $\text{U}_{\mathrm{PMNS}}$ as \text{U}_{\mathrm{PMNS}}^{\dagger}\mathcal{M}_{\nu}\text{U}_{\mathrm{PMNS}}^{\ast}=\mathrm{diag}\big{(}m_{\nu_{1}},m_{\nu_{2}},m_{\nu_{3}}\big{)}. The mixing angles in the PMNS matrix and neutrino mass eigenvalues are given by the global fitting to the neutrino oscillation data Gonzalez-Garcia:2014bfa .

As we will discuss in Sec. VI, the order of magnitude of the Yukawa couplings is ${\cal Y}_{jk}\sim{\cal O}(0.01-1)$ with $0.1\,\text{GeV}\lesssim M_{k}\lesssim 1\,\text{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_{\xi}\lesssim 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 $\mu_{2}\sim{\cal O}(100\,\,\text{MeV})$. Accordingly, if one takes ${\cal Y}_{jk}\sim 0.1,\,s_{\xi}\sim 0.05,\,\mu_{2}\sim 100\,\text{MeV},\,{\cal Y}^{L,R}_{jk}\sim 0.1$, and ${\cal C}{}^{L,R}_{jk}\sim 1$, the correct neutrino mass scale $m_{\nu}\sim 0.1\,\text{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 $\nu$SIMP scenario. Experimentally, the flavor-changing radiative decay $\ell_{r}\to\ell_{s}\gamma$ constrains the Yukawa couplings ${\cal Y}_{rk}$. The Feynman diagram depicted such decay process is shown in Fig. 2. The branching fraction of the decay process is

[TABLE]

where the fine structure constant $\alpha$, the Fermi constant $G_{\text{F}}$ and the loop function ${\cal F}(z)$ are

[TABLE]

with $\hat{e}$ the electromagnetic charge. The most stringent experimental limit on the $\mu\to e\gamma$ process comes from the MEG Collaboration TheMEG:2016wtm . The up-to-date upper bound on its branching ratio is $\,{\cal B}(\mu\to e\gamma)<4.2\times 10^{-13}$. If we take $m_{\eta^{+}}\sim 300\,\text{GeV}$ and ${\cal F}(M_{k}^{2}/m_{\eta^{+}}^{2})=1/6$ with $m_{\eta^{+}}\gg M_{k}$,333If $\eta^{\pm}$ decays dominantly into electron or muon, $m_{\eta^{+}}\gtrsim 270~{}\text{GeV}$ is required in order to avoid the constraint from the left-handed slepton search Aad:2014vma . the Yukawa couplings are limited in ${\cal Y}_{rk}\lesssim 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 ${\cal Y}$. In this solution, the pattern of neutrino mixing is pinned down by the structures of the other Yukawa matrices ${\cal Y}^{L,R}$ and the mass hierarchy of the vector-like fermions.444If $m_{\tau}>m_{e,\mu}+2M_{k}$, the new decay modes $\tau\to(e,\mu)\bar{N}N^{\prime}\to(e,\mu)\nu\nu^{\prime}\bar{X}X$ open and would contribute to $\tau\to(e,\mu)+\text{missing energy}$. However, this constraint is not so stringent.

At one-loop level, the presence of $\eta^{\pm}$ and $N_{k}$ also induces a contribution to the anomalous magnetic moment $a_{\ell_{j}}$ of charged lepton $\ell_{j}$ given by

[TABLE]

In particular, the current experimental value for the muon anomalous magnetic moment has more than $3\sigma$ deviation from the SM prediction  : $a^{\text{exp}}_{\mu}-a^{\text{SM}}_{\mu}=(288\pm 80)\times 10^{-11}$ Agashe:2014kda . Since the new contribution given by Eq.(22) is negative, we then require $|\Delta a_{\mu}|<8\times 10^{-10}$, this gives an upper bound for the Yukawa coupling ${\cal Y}_{rk}$. For example, by taking $m_{\eta^{+}}\sim 300\,\text{GeV}$ and $m_{\eta^{+}}\gg M_{k}$, the Yukawa couplings are limited in ${\cal Y}_{rk}\lesssim{\cal 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 $\Gamma_{Z,\,h\,\to\,\text{inv}}^{\text{new}}=\Gamma\scalebox{1.1}{(}Z,h\to X\bar{X}\scalebox{1.1}{)}$. The present experimental bounds on these invisible decay widths are $\Gamma_{Z\,\to\,\text{inv}}^{\text{new}}<2~{}\text{MeV}$ (at the 95% C.L.) Carena:2003aj and $\Gamma_{h\,\to\,\text{inv}}^{\text{new}}\lesssim 0.78~{}\text{MeV}$. To derive the latter one, we interpret ${\cal B}^{\text{new}}_{h\,\to\,\text{inv}}={\cal B}^{\text{exp}}_{\text{BSM}}<0.16$ Khachatryan:2016vau reported by the ATLAS and CMS combined measurements and adopt the SM Higgs width $\Gamma_{h}^{\text{SM}}=4.08~{}\text{MeV}$ Heinemeyer:2013tqa at $m_{h}=125.1\,\text{GeV}$  Aad:2015zhl . From Eqs.(II) and (II), these upper limits consequently are translated into $|s_{\xi}|\lesssim 0.4$ and $|s_{\xi}|\lesssim 0.165\,(m_{H}/100\,\text{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}\sim 300\,\text{GeV}$, we then reach the upper bound $|s_{\xi}|\lesssim 0.06$.

The scalar masses are constrained by the oblique parameters due to their modifications to the SM gauge boson propagators Peskin:1991sw . From Eq.(45) and Fig. 8 in Appendix B, those parameters are calculated as

[TABLE]

where the loop functions are given by

[TABLE]

[TABLE]

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 Burgess:1993mg . However, we have checked the difference is not important because the most stringent constraint comes from the ${\cal T}$-parameter whose definition does not change even for below electroweak scale. The current constraints are given in Ref. Baak:2012kk ; Baak:2014ora as $\Delta{\cal S}=0.05\pm 0.11$, $\Delta{\cal T}=0.09\pm 0.13$, $\Delta{\cal U}=0.01\pm 0.11$ with the correlation coefficients $0.90\,$(between $\Delta{\cal S}$ and $\Delta{\cal T}\,$), $-\,0.59\,$(between $\Delta{\cal S}$ and $\Delta{\cal U}$), and $-\,0.83\,$(between $\Delta{\cal T}$ and $\Delta{\cal U}$). These limits imply that the heavier neutral component $H$ and the charged component $\eta^{+}$ should be nearly degenerate \big{(}m_{H}\approx m_{\eta^{+}}\big{)} in the case of $s_{\xi}\ll 1$ and $m_{X}\ll m_{H},m_{\eta^{+}}.$

Theoretically, the quartic parameters $\lambda_{j}$ are subject to the conditions of vacuum stability and perturbativity. To ensure the vacuum to be stabilized at large field values, we demand Choi:2016tkj

[TABLE]

where $\lambda_{X}\equiv\lambda_{\chi}c_{\xi}^{4}\,\approx\lambda_{\chi}$, and $\lambda_{XS}\equiv\lambda_{\chi S}c_{\xi}^{2}\,\approx\lambda_{\chi S}$ because of the smallness of the mixing angle $\xi$. In the previous work Choi:2016tkj , a condition of perturbativity on the quartic couplings has been taken, which corresponds to $\lambda_{X,S}<16\pi$ 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 $\lambda_{X,S}<4\pi$ in our numerical work. Furthermore, since $X$ plays the role of DM, it should not develop the VEV. The sufficient conditions to guarantee $\langle X\rangle=0$ (as well as $\langle S\rangle=0$) are given by

[TABLE]

here we have assumed $\lambda_{3}=0$ for simplicity.555One can find the necessary conditions to ensure $\langle X\rangle=0$ by using the method in the literature Kannike:2012pe . However, the analytical result is too long to read.

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_{\text{DM}}=n_{X}+n_{\bar{X}}=2n_{X}$ (we assume there is no asymmetry between particles $X$ and $\bar{X}$) as follows

[TABLE]

with ${\cal H}$ being the Hubble parameter, $n_{\text{DM}}^{\text{eq}}$ the DM number density at the chemical equilibrium, and $\langle\sigma_{3\to 2}\upsilon_{\text{rel}}^{2}\rangle\equiv\frac{1}{24}\langle\sigma_{XXX\to\bar{X}\bar{X}}\upsilon_{\text{rel}}^{2}\rangle$ the thermal averaged effective $3\to 2$ annihilation cross section.666The definition of the effective $3\to 2$ annihilation cross section depends on the model. For example, in the $\mathbb{Z}_{3}$ SIMP model Choi:2015bya , \big{\langle}\sigma_{3\to 2}\upsilon_{\text{rel}}^{2}\big{\rangle}\equiv\frac{1}{24}\big{\langle}\sigma_{XXX\to X\bar{X}}\upsilon_{\text{rel}}^{2}\big{\rangle}+\frac{1}{8}\big{\langle}\sigma_{XX\bar{X}\to\bar{X}\bar{X}}\upsilon_{\text{rel}}^{2}\big{\rangle}.

By applying the standard derivation Kolb:1990vq , the approximate solution to the Boltzmann equation for the current relic density $\Omega_{\text{DM}}$ is given by

[TABLE]

where $x=m_{X}/T$, $\hat{h}$ denotes the normalized Hubble constant, $g_{\star,f}$ is the number of relativistic degrees of freedom at the freeze-out temperature, $T_{f}\,=\,m_{X}/x_{f}$, $m_{\text{pl}}=1.22\times 10^{19}\,\text{GeV}$ is the Planck mass. To evaluate the thermal average of the $3\to 2$ annihilation cross section, we employ the formula in Ref. Choi:2017mkk

[TABLE]

where \beta=\frac{1}{2}\big{(}\upsilon_{1}^{2}+\upsilon_{2}^{2}+\upsilon_{3}^{2}\big{)} with $\upsilon_{i}$ the velocities of three initial DM particles. In the $\nu$SIMP model, the Feynman diagrams of the $3\to 2$ process $XXX\to\bar{X}\bar{X}$ are shown in Fig. 3. From Eq.(7), the effective $3\to 2$ annihilation cross section under CP invariance is calculated as

[TABLE]

where \hat{s}\,=(p_{1}+p_{2}+p_{3})^{2}\,\simeq\,9m_{X}^{2}\big{(}1+2\beta/3\big{)} and the momenta of DM are neglected except around the resonance $\hat{s}\approx m_{S}^{2}$. By taking the mass spectrum $2M_{k}>m_{S}>2m_{X}$, the decay width of the particle $S$ is computed as

[TABLE]

To enhance the $3\to 2$ annihilation cross section, we pick the resonant pole $\,m_{S}\simeq\sqrt{\hat{s}}\simeq 3m_{X}$ in Eq.(30), and it is convenient to adjust the resonant behavior by defining the following dimensionless parameters as

[TABLE]

where $\epsilon_{S}$ indicates the degeneracy between $m_{S}$ and $3m_{X}$, and $\gamma_{S}$ is the width of the resonance.777From Eqs.(26), (31) and (32) with $\xi\ll 1$ and $m_{S}\simeq 3m_{X}$, one can easily show that $\gamma_{S}\,\simeq 10^{-3}{\cal R}_{2}^{2}\,\lesssim 10^{-2}\lambda_{X}$. Thus one obtains $\gamma_{S}\lesssim 0.1\ll 1$ with the perturbative bound $\lambda_{X}<4\pi$.

With these variables, the $3\to 2$ annihilation cross section can be expressed in the Breit-Wigner resonant form similar to the one in Ref.Gondolo:1990dk

[TABLE]

where the coefficient $c_{X}$ is

[TABLE]

with ${\cal R}_{1,2}=\mu_{1,2}/m_{X}$. Utilizing Eq.(28), we present the plots of the DM relic density $\Omega_{\text{DM}}$ versus $m_{S}$ with different values of $\,m_{X}$ and ${\cal 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\to 2$ annihilation cross section. The orange region is the latest relic density data $\,\Omega_{\text{DM}}\hat{h}^{2}=0.1197\pm 0.0022\,$ given by the Planck Collaboration Ade:2015xua . In these plots, we do not vary the mixing angle $\xi$ since dependence of the mixing angle is extremely small as long as $s_{\xi}\ll 1$. Also, in order to examine the conditions of $\langle X\rangle=0$ in Eq.(26) easily, we again assume $\lambda_{3}=0$. For nonzero $\lambda_{3}$, the numerical results are similar as pointed out in Ref. Choi:2016tkj . We have checked our choices of the values of $\,{\cal R}_{1,2}$ can accommodate the requirements of perturbativity, ${\cal R}_{1}^{2}<\lambda_{S}<4\pi$ and ${\cal R}_{2}^{2}/9\lesssim\lambda_{X}<4\pi$ with $m_{S}\simeq 3m_{X}$. According to the plots, one can see that the low values of ${\cal 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 Markevitch:2003at ; Clowe:2003tk ; Randall:2007ph and spherical halo shapes Peter:2012jh , which impose the bound $\sigma_{\text{self}}/m_{X}\lesssim 1\,\text{cm}^{2}/\text{g}$ with $\sigma_{\text{self}}=\frac{1}{4}(\sigma_{XX\to XX}+\sigma_{X\bar{X}\to X\bar{X}}+\sigma_{\bar{X}\bar{X}\to\bar{X}\bar{X}})$ the effective self-interacting cross section. We depict in Fig. 5 the Feynman diagrams of the DM self-interacting processes in our SIMP model, and their cross sections are calculated as

[TABLE]

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 \lambda_{X}\big{(}{\cal R}_{2}^{2}m_{X}^{2}/m_{S}^{2}<\lambda_{X}<4\pi\big{)}, the bounds from the Bullet cluster and spherical halo shapes can be satisfied. For instance, if we take $m_{X}(m_{S})=30\,(93)\,\text{MeV},{\cal R}_{1,2}=2,5$ and $\xi=0.05$ with $\lambda_{X}=7$, we find $\sigma_{\text{self}}/m_{X}\simeq 0.26\,\text{cm}^{2}/\text{g}$. More examples and discussions can be found in Ref. Choi:2016tkj .

VI SIMP condition

In the SIMP paradigm, DM is thermally produced through the $3\to 2$ annihilation process into the particles in the dark sector rather than the $2\to 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 Hochberg:2014dra

[TABLE]

which should be held during the freeze-out temperature. In this inequality, each reaction rate is defined by $\Gamma_{2\to 2}=n_{X}\langle\sigma_{2\to 2}\upsilon_{\text{rel}}\rangle,\,\Gamma_{3\to 2}=n_{X}^{2}\langle\sigma_{3\to 2}\upsilon_{\text{rel}}^{2}\rangle$, and $\Gamma_{\text{kin}}=n_{\text{SM}}\langle\sigma_{\text{kin}}\upsilon_{\text{rel}}\rangle$,888In the WIMP paradigm, the Boltzmann equation of the DM number density is given by

\dot{n}_{\text{DM}}+3{\cal H}n_{\text{DM}}=-\,\langle\sigma_{2\to 2}\upsilon_{\text{rel}}\rangle\scalebox{1.2}{\big{[}}n_{\text{DM}}^{2}-(n_{\text{DM}}^{\text{eq}}\big{)}{}^{2}\scalebox{1.2}{\big{]}}~{},\\

where $\langle\sigma_{2\to 2}\upsilon_{\text{rel}}\rangle$ is the thermal averaged effective $2\to 2$ annihilation cross section. Due to this definition, an extra factor $1/2$ is multiplied to the DM cross sections (Eq.(39) and (41)). This comes from the fact that DM and anti-DM particles are not identical in our case Gondolo:1990dk .

where the number densities of DM and the SM particles are given as Lee:2015gsa

[TABLE]

with $g$ counts the internal degrees of freedom, $m_{\text{SM}}$ being the mass of the SM particle, $(+)$ applies to fermions, and $(-)$ pertains to bosons. It is also pointed out in Ref. Kuflik:2015isi 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 $|\dot{E}_{3\to 2}|<|\dot{E}_{\text{kin}}|$, where $\dot{E}_{3\to 2}$ is the rate of the DM mass transferring to the kinetic energy in the DM bath, and $\dot{E}_{\text{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. Hochberg:2015vrg , we then impose $\Gamma_{3\to 2}<\Gamma_{3\to 2}<10^{-2}\Gamma_{\text{kin}}$ in our numerical study.

From Eq.(9), the particle $X$ can interact with the active neutrinos via the Yukawa couplings ${\cal Y}_{rk}$.999The particle $S$ can also have the $3\to 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. 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

[TABLE]

where the neutrino number density $n_{\nu}$ and the thermally averaged effective scattering cross section are given by

[TABLE]

with $\hat{\zeta}(3)\simeq 1.202$ the Riemann zeta function of 3.

By the crossing symmetry, the Feynman diagrams for the $2\to 2$ annihilation process of a DM pair into a pair of the SM neutrino are shown in Fig. 6(b),101010Here we have neglected the scattering processes $X\ell^{\pm}\to X\ell^{\pm}$ and the annihilation channels $X\bar{X}\to\ell^{+}\ell^{-}$ due to the mass suppression of the $Z$ boson and the Higgs boson. and the reaction rate is calculated as

[TABLE]

where the thermally averaged effective $2\to 2$ annihilation cross section is given by

[TABLE]

with $n_{X}{}$ given by Eq.(37). 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\to 2$ annihilation process and the elastic scattering are then reduced to the form as

[TABLE]

where \,\mathbb{Y}\equiv\scalebox{1.1}{\big{[}}\sum_{r,s,k,l}\text{Re}\big{(}{\cal Y}_{rk}^{\ast}{\cal Y}_{rl}{\cal Y}_{sk}{\cal Y}_{sl}^{\ast}\big{)}\scalebox{1.1}{\big{]}}^{1/4}. Using the SIMP condition at $T_{f}$, we illustrate the plots of the magnitude of $|\mathbb{Y}|$ as a function of $M$ in Fig. 7 with different numerical inputs based on Fig. 4. As indicated in the plots, the order of the Yukawa coupling is $|\mathbb{Y}|\sim{\cal O}(0.01-1)$ with $0.1\,\text{GeV}\lesssim M\lesssim 1\,\text{GeV}$.

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 ${\mathbb{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\to 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.

Acknowledgments

T. T. acknowledges support from JSPS Fellowships for Research Abroad. The work of K. T is supported by JSPS Grant-in-Aid for Young Scientists (B) (Grants No. 16K17697), by the MEXT Grant-in-Aid for Scientific Research on Innovation Areas (Grants No. 16H00868), and by Kyoto University: Supporting Program for Interaction-based Initiative Team Studies (SPIRITS).

Appendix A Gauged $\text{U}(1)_{\text{B}-\text{L}}$ extension of the $\nu$SIMP model

It is believed that there is no global symmetry can exist in a theory of quantum gravity Ibanez:1991pr ; Rai:1992xw . Under this context, the discrete symmetry we introduced in our $\nu$SIMP model may originate from a gauge symmetry (gauge redundancy). At certain energy scale, this gauge symmetry is broken down to the ${\mathbb{Z}_{5}}$ discrete symmetry by a nonzero VEV of a scalar field. In the following, we demonstrate an extension of the $\nu$SIMP model to the gauged $\text{U}(1)_{\text{B}-\text{L}}$ version by adding one more SM singlet complex scalar $\zeta$. The particle contents and the charge assignments are summarized in Tab. 2.

In this extended model, the Lagrangian associated with the $3\to 2$ processes is given by

[TABLE]

After spontaneous symmetry breaking, the complex scalar $\zeta$ can be expanded around its VEV as \zeta=\tfrac{1}{\sqrt{2}}\big{(}\varsigma+\upsilon^{\prime}\big{)}, where $\upsilon^{\prime}\equiv\sqrt{2}\langle\zeta\rangle$. The scalar interactions between $\chi$ and $S$ are then extracted as

[TABLE]

which corresponding to the first three terms in the last line of Eq.(2), respectively, with $\mu_{1,2}=\lambda_{1,2}\upsilon^{\prime}$. 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^{\prime}$-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.

Appendix B Gauge interactions

The kinetic part of the Lagrangian in Eq.(1) contains the interactions of the new scalars with the photon and the weak bosons,

[TABLE]

where

[TABLE]

With these gauge interactions, we draw the Feynman diagrams of the contributions to the SM gauge boson propagators in Fig. 8.

Appendix C Benchmark points

Assuming ${\cal Y}^{L}\gg{\cal Y}^{R}$ and the other parameter set

[TABLE]

two benchmark Yukawa couplings are given as

[TABLE]

for $m_{X}=30~{}\text{MeV}$, $m_{S}=93~{}\text{MeV}$, $\mu_{2}=150~{}\text{MeV}$, and

[TABLE]

for $m_{X}=40~{}\text{MeV}$, $m_{S}=128~{}\text{MeV}$, $\mu_{2}=320~{}\text{MeV}$. One can check that Eq. (53) and (60) satisfy the SIMP condition $(\Gamma_{3\to 2}/\Gamma_{2\to 2}\,\simeq 10^{4}$ and  $\Gamma_{\text{kin}}/\Gamma_{3\to 2}\,\simeq 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 ${\cal Y}^{L}\ll{\cal Y}^{R}$ and inverted hierarchy, though these are not shown here.