Radiative neutrino model with semi-annihilation dark matter

TL;DR

This paper introduces a two-loop radiative neutrino model with a hidden gauged U(1) symmetry, where semi-annihilation significantly influences dark matter relic density, aligning with neutrino and lepton flavor constraints.

Contribution

It presents a novel two-loop radiative neutrino model incorporating semi-annihilation dark matter and explores its parameter space consistent with experimental bounds.

Findings

Semi-annihilation is crucial for relic density.

Allowed DM mass range is 80-400 GeV.

Model satisfies neutrino oscillation and lepton flavor violation constraints.

Abstract

We propose a Two-Loop induced radiative neutrino model with hidden gauged $U(1)$ symmetry, in which a dark matter of Dirac fermion arises. The relic density gets contribution from annihilation and semi-annihilation due to a residual $\mathbb{Z}_3$ parity. After imposing the requirement of neutrino oscillation data and lepton flavour violation bounds, we find out that the semi-annihilation plays a crucial role in order to satisfy the relic density constraint $0.117 < \Omega h^2 < 0.123$, by proceeding near either one of two deconstructive scalar resonances. Our numerical analysis demonstrates the allowed region for the DM-Scalar coupling with the DM mass in $(80, 400)$ GeV.

Tables

Table 1. Table 1: Contents of fermion and scalar fields and their charge assignments under S ​ U ​ ( 2 ) L × U ​ ( 1 ) Y × U ​ ( 1 ) H 𝑆 𝑈 subscript 2 𝐿 𝑈 subscript 1 𝑌 𝑈 subscript 1 𝐻 SU(2)_{L}\times U(1)_{Y}\times U(1)_{H} , where all the new fields are singlets under S ​ U ​ ( 3 ) C 𝑆 𝑈 subscript 3 𝐶 SU(3)_{C} , and all the quark fields are neutral under U ​ ( 1 ) H 𝑈 subscript 1 𝐻 U(1)_{H} . Note that the H ′ superscript 𝐻 ′ H^{\prime} field can only present for the one-loop radiative seesaw.

	Fermion Fields					Scalar Fields				Inert Scalar Fields
	$L_L$	$e_R$	$L_{L/R}^{\prime }$	${\chi }_{L/R}$	$N_{L/R}$	$H$	$H^{\prime }$	$\Delta$	$\phi$	$s$	$\eta$	$s^{\prime }$	${\eta }^{\prime }$
$SU{\left(2\right)}_L$	$\mathrm{𝟐}$	$\mathrm{𝟏}$	$\mathrm{𝟐}$	$\mathrm{𝟏}$	$\mathrm{𝟏}$	$\mathrm{𝟐}$	$\mathrm{𝟐}$	$\mathrm{𝟑}$	$\mathrm{𝟏}$	$\mathrm{𝟏}$	$\mathrm{𝟐}$	$\mathrm{𝟏}$	$\mathrm{𝟐}$
$U{\left(1\right)}_Y$	$-\frac{1}{2}$	$-1$	$-\frac{1}{2}$	$0$	$0$	$\frac{1}{2}$	$\frac{1}{2}$	$1$	0	0	$\frac{1}{2}$	0	$\frac{1}{2}$
$U{\left(1\right)}_H$	$0$	$0$	$2x$	$x$	$y$	$0$	$-3x$	$-3x$	$-3x$	$-2x$	$x$	$x+y$	$-2x+y$

Equations

$$-{\mathcal{L}}_Y$$

$$+y_{{\eta }_{ab}^{\prime }}{\bar{L}}_{R_a}^{\prime }{\tilde{\eta }}^{\prime }{N}_{L_b}+y_{{\eta }_{ab}^{\prime }}^{\prime }{\bar{L}}_{L_a}^{\prime }{\tilde{\eta }}^{\prime }{N}_{R_b}+y_{s_{ab}^{\prime }}{\bar{N}}_{R_a}{\chi }_{R_b}^c{s}^{\prime }+y_{s_{ab}^{\prime }}^{\prime }{\bar{N}}_{L_a}{\chi }_{L_b}^c{s}^{\prime }$$

$$+M_{{\chi }_{aa}}{\bar{\chi }}_{L_a}{\chi }_{R_a}+M_{N_{aa}}{\bar{N}}_{L_a}{N}_{R_a}+M_{N_{aa}^{\prime }}{\bar{L}}_{L_a}^{\prime }{L}_{R_a}^{\prime }+h.c.,$$

$${\mathcal{V}}_1$$

$$+\sum _{\varphi }^{\eta ,{\eta }^{\prime }}{\lambda }_{H\varphi }^{\prime }(H^{†}\varphi )({\varphi }^{†}H)+\sum _{\varphi }^{\eta ,{\eta }^{\prime },s,s^{\prime }}\left[{\mu }_{\varphi }^2{\varphi }^{†}\varphi +{\lambda }_{\varphi }|{\varphi }^{†}\varphi {|}^2\right].$$

$${\mathcal{V}}_2$$

$$\displaystyle H=\left[\begin{array}[]{c}G^{+}\\ \frac{v_{H}+h+iG^{0}}{\sqrt{2}}\end{array}\right]\,,\quad\Delta=\left[\begin{array}[]{cc}\frac{\Delta^{+}}{\sqrt{2}}&\Delta^{++}\\ \frac{v_{\Delta}+\Delta_{R}+i\Delta_{I}}{\sqrt{2}}&-\frac{\Delta^{+}}{\sqrt{2}}\end{array}\right]\,,\quad\varphi=\frac{v_{\varphi}+\varphi_{R}+i\varphi_{I}}{\sqrt{2}}\,.$$

$$-{\mu }_H^2+{\lambda }_H{v}_H^2+\frac{{\lambda }_3+{\lambda }_4}{2}{v}_{\Delta }^2+\frac{{\lambda }_5}{2}{v}_{\phi }^2={\lambda }_{\Delta }{v}_{\Delta }{v}_{\phi }$$

$$(M^2+\frac{{\lambda }_3+{\lambda }_4}{2}{v}_H^2+({\lambda }_1+{\lambda }_2)v_{\Delta }^2+\frac{{\lambda }_6}{2}{v}_{\phi }^2)v_{\Delta }={\lambda }_{\Delta }{v}_H^2{v}_{\phi }/2$$

$$(-{\mu }_{\phi }^2+{\lambda }_{\phi }{v}_{\phi }^2+\frac{{\lambda }_5}{2}{v}_H^2++\frac{{\lambda }_6}{2}{v}_{\Delta }^2)v_{\phi }={\lambda }_{\Delta }{v}_H^2{v}_{\Delta }/2$$

$$v_H\simeq {\left(\frac{{\lambda }_{\Delta }{v}_{\Delta }{v}_{\phi }+{\mu }_H^2}{{\lambda }_H}\right)}^{1/2},v_{\Delta }\simeq \frac{{\lambda }_{\Delta }{v}_H^2{v}_{\phi }}{2(M^2+({\lambda }_3+{\lambda }_4)v_H^2/2)},v_{\phi }\simeq \frac{{\mu }_{\phi }}{{\lambda }_{\phi }^{1/2}}.$$

$$\displaystyle\eta=\left[\begin{array}[]{c}\eta^{+}\\ \frac{\eta_{R}+i\eta_{I}}{\sqrt{2}}\end{array}\right]\,,\quad s=\frac{s_{R}+is_{I}}{\sqrt{2}}\,;\quad\eta^{\prime}=\left[\begin{array}[]{c}\eta^{\prime+}\\ \frac{\eta^{\prime}_{R}+i\eta^{\prime}_{I}}{\sqrt{2}}\end{array}\right]\,,\quad s^{\prime}=\frac{s^{\prime}_{R}+is^{\prime}_{I}}{\sqrt{2}}\,.$$

$$\displaystyle V_{\alpha}^{T}\left[\begin{array}[]{cc}m_{s_{R}}^{2}&\frac{\lambda_{0}}{2}v_{H}v_{\varphi}\\ \frac{\lambda_{0}}{2}v_{H}v_{\varphi}&m_{\eta_{R}}^{2}\\ \end{array}\right]V_{\alpha}=\left[\begin{array}[]{cc}m^{2}_{H_{1}}&0\\ 0&m^{2}_{H_{2}}\end{array}\right],~{}~{}V_{\alpha^{\prime}}^{T}\left[\begin{array}[]{cc}m_{s^{\prime}_{R}}^{2}&\frac{\lambda^{\prime}_{0}}{2}v_{H}v_{\varphi}\\ \frac{\lambda^{\prime}_{0}}{2}v_{H}v_{\varphi}&m_{\eta^{\prime}_{R}}^{2}\\ \end{array}\right]V_{\alpha^{\prime}}=\left[\begin{array}[]{cc}m^{2}_{H^{\prime}_{1}}&0\\ 0&m^{2}_{H^{\prime}_{2}}\end{array}\right]$$

$$\displaystyle\left[\begin{array}[]{c}s_{R}+is_{I}\\ \eta_{R}+i\eta_{I}\end{array}\right]=\left[\begin{array}[]{cc}c_{\alpha}&-s_{\alpha}\\ s_{\alpha}&c_{\alpha}\end{array}\right]\left[\begin{array}[]{c}H_{1}+iA_{1}\\ H_{2}+iA_{2}\end{array}\right],\quad\left[\begin{array}[]{c}s^{\prime}_{R}+is^{\prime}_{I}\\ \eta^{\prime}_{R}+i\eta^{\prime}_{I}\end{array}\right]=\left[\begin{array}[]{cc}c_{\alpha^{\prime}}&-s_{\alpha^{\prime}}\\ s_{\alpha^{\prime}}&c_{\alpha^{\prime}}\end{array}\right]\left[\begin{array}[]{c}H^{\prime}_{1}+iA^{\prime}_{1}\\ H^{\prime}_{2}+iA^{\prime}_{2}\end{array}\right]$$

$$-{\mathcal{L}}_Y$$

$$+\frac{y_{{\chi }_{ab}}^{L/R}}{\sqrt{2}}{\bar{\chi }}_{L_a/R_a}^C{\chi }_{L_a/R_b}(c_{\alpha }{H}_1-s_{\alpha }{H}_2)+i\frac{y_{{\chi }_{ab}}^{L/R}}{\sqrt{2}}{\bar{\chi }}_{L_a/R_a}^C{\chi }_{L_a/R_b}(c_{\alpha }{A}_1-s_{\alpha }{A}_2)$$

$$+\frac{y_{S_{ia}}}{\sqrt{2}}{\bar{\nu }}_{L_i}{N}_{R_a}^{\prime }(c_{\alpha }{H}_1-s_{\alpha }{H}_2)+i\frac{y_{S_{ia}}}{\sqrt{2}}{\bar{\nu }}_{L_i}{N}_{R_a}^{\prime }(c_{\alpha }{A}_1-s_{\alpha }{A}_2)$$

$$+\frac{y_{s_{ab}^{\prime }}}{\sqrt{2}}{\bar{N}}_{R_a}{\chi }_{R_b}^C(c_{{\alpha }^{\prime }}{H}_1^{\prime }-s_{{\alpha }^{\prime }}{H}_2^{\prime })+i\frac{y_{s_{ab}^{\prime }}}{\sqrt{2}}{\bar{N}}_{R_a}{\chi }_{R_b}^C(c_{{\alpha }^{\prime }}{A}_1^{\prime }-s_{{\alpha }^{\prime }}{A}_2^{\prime })$$

$$+\frac{y_{s_{ab}^{\prime }}^{\prime }}{\sqrt{2}}{\bar{N}}_{L_a}{\chi }_{L_b}^C(c_{{\alpha }^{\prime }}{H}_1^{\prime }-s_{{\alpha }^{\prime }}{H}_2^{\prime })+i\frac{y_{s_{ab}^{\prime }}^{\prime }}{\sqrt{2}}{\bar{N}}_{L_a}{\chi }_{L_b}^C(c_{{\alpha }^{\prime }}{A}_1^{\prime }-s_{{\alpha }^{\prime }}{A}_2^{\prime })$$

$$+\frac{y_{{\eta }_{ab}^{\prime }}}{\sqrt{2}}{\bar{N}}_{R_a}^{\prime }{N}_{L_b}(s_{{\alpha }^{\prime }}{H}_1^{\prime }+c_{{\alpha }^{\prime }}{H}_2^{\prime })-i\frac{y_{{\eta }_{ab}^{\prime }}}{\sqrt{2}}{\bar{N}}_{R_a}^{\prime }{N}_{L_b}(s_{{\alpha }^{\prime }}{A}_1^{\prime }+c_{{\alpha }^{\prime }}{A}_2^{\prime })$$

$$+\frac{y_{{\eta }_{ab}^{\prime }}^{\prime }}{\sqrt{2}}{\bar{N}}_{L_a}^{\prime }{N}_{R_b}(s_{{\alpha }^{\prime }}{H}_1^{\prime }+c_{{\alpha }^{\prime }}{H}_2^{\prime })-i\frac{y_{{\eta }_{ab}^{\prime }}^{\prime }}{\sqrt{2}}{\bar{N}}_{L_a}^{\prime }{N}_{R_b}(s_{{\alpha }^{\prime }}{A}_1^{\prime }+c_{{\alpha }^{\prime }}{A}_2^{\prime })+h.c..$$

$$(m_{\nu }{)}_{ij}=m_{{\nu }_{ij}}^{(I)}+m_{{\nu }_{ij}}^{(II)}+[m_{{\nu }_{ij}}^{(I)}{]}^T+[m_{{\nu }_{ij}}^{(II)}{]}^T,$$

$$V_{\mathrm{MNS}}$$

$$\times \mathrm{diag}(1,{\mathrm{e}}^{\mathrm{i}\frac{{\alpha }_{21}}{2}},{\mathrm{e}}^{\mathrm{i}\frac{{\alpha }_{31}}{2}})$$

$$0.250\le s_{12}^2\le 0.354,0.381\le s_{23}^2\le 0.615,0.019\le s_{13}^2\le 0.024,$$

$$m_{{\nu }_3}^2-m_{{\nu }_1}^2=(2.45-2.69)\times 10^{-3}{\mathrm{eV}}^2,$$

$$m_{{\nu }_2}^2-m_{{\nu }_1}^2=(6.93-7.96)\times 10^{-5}{\mathrm{eV}}^2,$$

$$(m_{\nu }{)}_{ij}\equiv \frac{1}{(4\pi {)}^4}\left(y_{{\eta }_{ia}}[F_I+F_{II}{]}_{ab}{y}_{S_{bj}}^T+y_{S_{ja}}[F_I^T+F_{II}^T{]}_{ab}{y}_{{\eta }_{bj}}^T\right)$$

$$\equiv \frac{1}{(4\pi {)}^4}\left(y_{{\eta }_{ia}}{G}_{ab}{y}_{S_{bj}}^T+y_{S_{ja}}{G}_{ab}^T{y}_{{\eta }_{bj}}^T\right),$$

$$y_{\eta }=\frac{1}{2}[(V_{\mathrm{MNS}}^{\ast }{D}_{\nu }{V}_{\mathrm{MNS}}^{†}+A](y_S^T{)}^{-1}{G}^{-1},$$

$$\mathcal{L}=-y_{{\eta }_{ia}}{\bar{\ell }}_{L_i}{\eta }^{-}{\chi }_{R_a}+\frac{1}{\sqrt{2}}{y}_{S_{ia}}{\bar{\ell }}_{L_i}{E}_{R_a}^{\prime }[(c_{{\alpha }_R}{H}_1-s_{{\alpha }_R}{H}_2)+i(c_{{\alpha }_I}{A}_1-s_{{\alpha }_I}{A}_2)]+h.c.,$$

$$Br({\ell }_i\to {\ell }_j\gamma )\approx \frac{48{\pi }^3{\alpha }_{\mathrm{em}}}{G_{\mathrm{F}}^2{m}_{{\ell }_i}^2}{C}_{ij}\left(|a_{L_{ij}}{|}^2+|a_{R_{ij}}{|}^2\right),$$

$$a_{R_{ij}}\approx \frac{m_{{\ell }_i}}{(4\pi {)}^2}[y_{{\eta }_{ja}}{y}_{{\eta }_{ai}}^{†}H({\chi }_a,{\eta }^{-})$$

$$-\frac{y_{S_{ja}}{y}_{S_{ai}}^{†}}{2}[c_{{\alpha }_R}^{2}H(H_1,E_a^{\prime })+s_{{\alpha }_R}^{2}H(H_2,E_a^{\prime })+c_{{\alpha }_I}^{2}H(A_1,E_a^{\prime })+s_{{\alpha }_I}^{2}H(A_2,E_a^{\prime })]],$$

$$H(a,b)={\int }_0^{1}dx{\int }_0^{1-x}dy\frac{xy}{x{m}_a^2+(1-x)m_b^2}$$

$$\stackrel{m_a\ne m_b}{=}\frac{2m_a^6+3m_a^4{m}_b^2-6m_a^2{m}_b^4+m_b^6+6m_a^4{m}_b^2\log \left(\frac{m_b^2}{m_a^2}\right)}{12{\left(m_a^2-m_b^2\right)}^4}$$

$$Br(\mu \to e\gamma )\le 4.2\times 10^{-13}(6\times 10^{-14})$$

$$Br(\tau \to \mu \gamma )\le 4.4\times 10^{-8},Br(\tau \to e\gamma )\le 3.3\times 10^{-8}$$

$$\Delta a_{\mu }\approx -m_{\mu }(a_L+a_R{)}_{22}.$$

$$\displaystyle\frac{1}{2}\begin{pmatrix}Z_{0}\\ \tilde{Z}\end{pmatrix}^{T}\left[\begin{array}[]{cc}\frac{(g_{1}^{2}+g_{2}^{2})}{4}(v_{H}^{2}+4v_{\Delta}^{2})&3x\sqrt{g_{1}^{2}+g_{2}^{2}}g_{H}v_{\Delta}^{2}\\ 3x\sqrt{g_{1}^{2}+g_{2}^{2}}g_{H}v_{\Delta}^{2}&9x^{2}g_{H}^{2}(v_{\Delta}^{2}+v_{\varphi}^{2})\\ \end{array}\right]\begin{pmatrix}Z_{0}\\ \tilde{Z}\end{pmatrix},$$

Full text

I Introduction

Radiative seesaw neutrino models are one of the attractive scenarios to connect neutrino sector with dark matter (DM) sector in a natural manner. These two sectors certainly involve mysterious puzzles that are frequently interpreted as physics beyond the Standard Model (SM). When the neutrino masses are radiatively induced, the magnitude of relevant couplings could reach $\mathcal{O}(1)$ compared with the case where the neutrino mass is generated at the tree-level, so that the mass hierarchy among the SM sector and heavy fermion/scalar sectors is largely alleviated. Furthermore, new particles that are accommodated in the theory are at $\mathcal{O}$(TeV) energy scale and accessible by the extensive search at Large Hadron Collider (LHC). For radiative seesaw, a discrete symmetry is essentially implemented in order to forbid the neutrino mass at the tree-level and such symmetry will in turn stabilize the lightest neutral particle as a DM candidate. As a consequence, this type of theory provides interesting phenomenologies with the requirement to satisfy the observed relic density of $\Omega h^{2}\approx 0.120\pm 0.001$ planck and other experimental constraints.

The simplest discrete symmetry can be $\mathbb{Z}_{2}$, as the remnant of a broken $U(1)$ symmetry, and a typical DM-generated neutrino mass model at the one-loop level is proposed in Ma:2006km . However other enlarged discrete symmetry is also possible to stablize the DM candidate such as the $\mathbb{Z}_{N}$, $N>2$ discrete parity, which brings in semi-annihilation in addition to annihilation for the Lee-Weinberg scenario Lee:1977ua , allowing for an odd number of DM particles appearing in a $2\to 2$ process Hambye:2008bq ; Hambye:2009fg ; DEramo:2010keq ; Belanger:2012vp . Under the control of $\mathbb{Z}_{N}$ discrete symmetry, any field transforming as $f_{i}\to\omega^{a}f_{i}$, with $\omega=\exp(i2\pi/N)$ and $a=1,\cdots,N-1$, could serve as the dark matter candidate depending on the spectrum and interactions. In this paper we consider a two-loop induced neutrino mass model 2-lp-zB ; Babu:2002uu ; Ma:2007gq ; Kajiyama:2013zla ; Kajiyama:2013rla ; Aoki:2013gzs with new particles charged under a hidden $U(1)$ symmetry Langacker:2008yv ; Ma:2013yga ; Ko:2014loa ; Ma:2015mjd ; Chun:2008by ; Ko:2016ala ; Nomura:2018jkd , in which a Dirac fermion type of $\mathbb{Z}_{3}$ DM candidate arises, whose relic density is dominantly explained by the s-channel of semi-annihilation modes. Note that in this model, a complex scalar is also possible to behave as DM in the inverse mass pattern. The discrete $\mathbb{Z}_{3}$ symmetry origins from the spontaneous breaking of $U(1)$ symmetry and plays an important role to ensure the DM $\chi$ does not decay while the reaction in the form of $\chi\chi\to\chi^{\dagger}v_{i}$ exits. We present how the DM and neutrino mass are correlated by formulating each sector. In particular, we perform an analysis to obtain the allowed region which satisfies a set of necessary bounds, including neutrino oscillation data, Lepton Flavour Violations (LFVs), muon anomalous magnetic moment ($\Delta a_{\mu}$, aka muon $g-2$), and the DM relic density.

This paper is organized as follows. In Sec. II, we show the valid Lagrangian with charge assignments, and formulate the scalar and neutrino sectors, along with the LFVs, muon $g-2$, $Z-Z^{\prime}$ mixing and bound of electroweak precision test. In Sec. III, we analyze the Dirac fermionic DM to explain the relic density with an emphasis on the semi-annihilation and a brief illustration of the analytic derivation. In Sec. IV, we conduct a numerical analysis, and show the allowed region to satisfy all the phenomenologies that we discuss above. Finally we conclude and discuss in Sec. V.

II The Model

The model is built by extending the SM with additional scalars and vector-like fermions, which are charged under a hidden $U(1)$ symmetry before some of the scalars obtain VEVs. The field contents and their charge assignments are reported in Table 1. For the fermion sector, an isospin doublet $L^{\prime}\equiv[E^{\prime},N^{\prime}]^{T}_{i}$ plus two isospin singlets $\chi_{i}$ and $N_{i}$ with $i=1,2,3$, are added. The vector-like nature of these extra fermions ensures our extension to be anomaly-free. The quantum number assignment for $L^{\prime}$, $\chi$, $N$ under the two gauge groups of $(U(1)_{Y},U(1)_{H})$ are $(-1/2,2x)$, $(0,x)$ and $(0,y)$ respectively. Here we use two arbitrary integers $(x,y)$ with $\{x,y\}\neq 0$ to keep track of the heavy fermions running in the outer and inner loops of neutrino mass diagram (see Figure 1). As for new scalar fields, we introduce four inert scalar fields $s$, $\eta$, $s^{\prime}$, $\eta^{\prime}$, where $(\eta,\eta^{\prime})$ are $SU(2)_{L}$ doublets and $(s,s^{\prime})$ are singlets. As we can see that since $(s,\eta)$ are charged under the $U(1)_{H}$ as $(-2x,x)$, these two fields will only interact with new fermions of $L^{\prime}$ and $\chi$. On the other hand, the two prime fields $(s^{\prime},\eta^{\prime})$ are charged with $(x+y,-2x+y)$ for the hidden symmetry, thus they are allowed to connect with the exotic fermion $N$ under the assumption of $y\neq-x$, $y\neq 2x$ 111In fact we can think that the gauged $U(1)_{H}$ is a linear combination of two global $U(1)$s, which should be observed individually in the unbroken phase.. The two scalar fields $(H,\varphi)$ are needed in order to generate the neutrino mass at the two-loop level provided they will mix inert scalars $(\eta,s)$ and $(\eta^{\prime},s^{\prime})$ inside each set. The symmetry permits more scalar fields, like a doublet $H^{\prime}$ or a triplet $\Delta$, to induce the $Z$-$Z^{\prime}$ mixing for LHC collider signature. In case of adding $H^{\prime}$, the neutrino mass will be generated at the one-loop level, since the red dot in the Figure.1 can be substituted by an interaction of $\overline{L^{\prime c}}_{L/R}H^{\prime}\chi_{L/R}$. However in such case the VEV of $H^{\prime}$ should be very small (equivalent to loop generated), so that this one-loop radiative seesaw is possible to reconcile the tension between neutrino mass and relic density bound. Thus we will focus on exploring the impact of a triplet $\Delta$ interplaying with $(H,\varphi)$ via the scalar potential in the two-loop radiative seesaw. For that purpose, the scalars $H$, $\Delta$ and $\varphi$ are required to develop nonzero vacuum expectation values (VEVs), respectively symbolized by $\langle H\rangle\equiv v_{H}/\sqrt{2}$, $\langle\Delta\rangle\equiv v_{\Delta}/\sqrt{2}$, $\langle\varphi\rangle\equiv v_{\varphi}/\sqrt{2}$. The valid renormalizable Lagrangian for the fermion sector are given by,

[TABLE]

where $i,a,b=1,2,3$ are the flavor indices for the SM and exotic fermions, and $\tilde{\eta}\equiv i\sigma_{2}\eta^{*}$, with $\sigma_{2}$ being the second Pauli matrix. For simplicity, we assume that all coefficients are real and $M_{\chi},M_{N},M_{N^{\prime}}$ to be diagonal matrices. The first term of $\mathcal{L}_{Y}$ generates the SM charged-lepton masses $m_{\ell_{i}}\equiv y_{\ell_{i}}v_{H}/\sqrt{2}$, while the $2$nd to $4$th terms will be responsible for the (semi-)annihilations. The residual $\mathbb{Z}_{3}$ from the broken hidden symmetry makes the lightest neutral states with $U(1)_{H}$ charge of $x$ or $\pm 2x$, i.e. particles in the set of $(\chi_{i},N^{\prime}_{i},\eta,s)$, to be our DM candidate. While in this paper we are interested in the mass pattern where $\chi_{1}$ actually plays the role of DM. Referring to Table 1, we can see that the two scalar fields $(\Delta,\varphi)$ carrying a $U(1)_{H}$ charge $q_{H}=-3x$ with $x\in$ integer, so that they will transform under the Abelian $U(1)$ symmetry as $\Delta\to e^{-iq_{H}\alpha}\Delta$ and $\varphi\to e^{-iq_{H}\alpha}\varphi$, for an arbitrary value of $\alpha$ before the spontaneous symmetry breaking. However after these two scalars obtain VEVs, the phase is forced to be $\alpha=2\pi/3$ for any integer $x\geq 1$, thus the Lagrangian is still invariant under a discrete $\mathbb{Z}_{3}$ symmetry. And the particles with $x(2x)$ charge in $U(1)_{H}$ have $w=e^{i2\pi/3}(w^{2})$ parity assignment under this $\mathbb{Z}_{3}$.

II.1 The scalar potential

We explicitly write the nontrivial terms for the inert scalar potential which are invariant under the $SU(2)_{L}\times U(1)_{Y}\times U(1)_{H}$ gauge symmetry to be:

[TABLE]

where we assume that these terms like $s^{\prime 2}\varphi^{(*)}$, $s^{\prime 2}\varphi^{2}$, $\eta^{\prime{\dagger}}\eta\varphi^{2}$ vanish due to $U(1)_{H}$ charges (e.g. $x=y=1$). Thus no mass splitting occurs among the real and imaginary parts of any inert field. The general potential for the scalars $(H,\Delta)$ can be found in ref Bonilla:2015eha ; Primulando:2019evb , and we will modify it by adding interactions with a complex singlet $\varphi$.

[TABLE]

The scalar fields beside the inert ones are explicitly expressed as:

[TABLE]

so that the mass of $W$ boson is fixed to be $m_{W}=\frac{g_{2}\sqrt{v_{H}^{2}+2v_{\Delta}^{2}}}{2}$. The minimum of the potential is determined by derivatives $\partial\mathcal{V}_{2}/\partial v_{H}=0$, $\partial\mathcal{V}_{2}/\partial v_{\Delta}=0$, $\partial\mathcal{V}_{2}/\partial v_{\varphi}=0$, which read as:

[TABLE]

As we argue in the section [II.4] for $Z-Z^{\prime}$ mixing, $v_{\Delta}$ is very tiny due to the $\rho$ parameter, thus we will focus on the limit of $v_{\Delta}\ll v_{H}\lesssim v_{\varphi}$. Thus under the assumption of negligible mixing between $\varphi$ and $(H,\Delta)$, i.e. $\lambda_{5},\lambda_{6}\ll 1$, we obtain:

[TABLE]

In addition the mass matrices in terms of $(h,\Delta_{R},\varphi_{R})$, $(G^{0},\Delta_{I},\varphi_{I})$ and $(G^{+},\Delta^{+})$ can be diagonalised into CP-even or odd spectrum by respective orthogonal matrices. Analogously the inert bosons $(s,\eta)_{R/I}$ and $(s^{\prime},\eta^{\prime})_{R/I}$ are written as:

[TABLE]

They are rotated into the mass basis as follows:

[TABLE]

where we use shorthands of $s_{\alpha^{(^{\prime})}}=\sin\alpha^{(^{\prime})}$, $c_{\alpha^{(^{\prime})}}=\cos\alpha^{(^{\prime})}$ and the complex fields $H_{i}+iA_{i}$, $H^{\prime}_{i}+iA^{\prime}_{i}$, $i=1,2$ are mass eigenstates. Note that the semi-annihilation exists for the theory with a $\mathbb{Z}_{3}$ parity, indicating that we need to keep the degeneracy between $H_{1,(2)}$ and $A_{1,(2)}$. The reason is that a $\mathbb{Z}_{3}$ parity assignment $w=e^{i2\pi/3}$ is valid for a Dirac fermion or a complex scalar field, like $\tilde{H}_{i}=H_{i}+iA_{i}$ with $i=1,2$. Under this specific potential we obtain that: $m^{2}_{H_{1,2}}=m^{2}_{A_{1,2}}$, $m^{2}_{H^{\prime}_{1,2}}=m^{2}_{A^{\prime}_{1,2}}$. Without loss of generality, we can assume $m_{H_{1}}<m_{H_{2}}$ and $m_{H^{\prime}_{1}}<m_{H^{\prime}_{2}}$ by ordering the mass eigenstates.

II.2 Neutrino mass matrix

In this model, the neutrino mass arises at the 2-loop level. To facilitate the calculation, the Lagrangian should be transformed into the mass basis:

[TABLE]

Here we assume that all the Yukawa couplings are real for simplicity. The active neutrino mass matrix $m_{\nu_{ij}}$ are generated at two-loop level as shown in Figure 2, with their formulas given by

[TABLE]

where $m_{\nu_{ab}}^{(I)}$ and $m_{\nu_{ab}}^{(II)}$ respectively correspond to the left and right plots in Figure 2. The constraint on the neutrino matrix is from the neutrino oscillation data, since $({m}_{\nu})_{ab}$ have to be diagonalized by the Pontecorvo-Maki-Nakagawa-Sakata mixing matrix $V_{\rm MNS}$ (PMNS) Maki:1962mu as $({m}_{\nu})_{ij}=(V_{\rm MNS}^{*}D_{\nu}V_{\rm MNS}^{\dagger})_{ij}$ with $D_{\nu}=\rm{diag}(m_{\nu_{1}},m_{\nu_{2}},m_{\nu_{3}})$. The PMNS matrix is parametrised as:

[TABLE]

with $s_{ij}=\sin\theta_{ij}$ being three mixing angles. In the following analysis, we will also neglect the Majorana CP violation phase $\alpha_{21}$ and $\alpha_{31}$ as well as Dirac CP violation phase $\delta$. By assuming the normal mass order $m_{\nu 1}\ll m_{\nu 2}<m_{\nu 3}$, the global fit of the current experiments at $3\sigma$ is given by pdg2018 :

[TABLE]

Now we rewrite the neutrino mass matrix in terms of Yukawa couplings and the form factors:

[TABLE]

where the factor $\frac{1}{(4\pi)^{4}}$ comes from the loop integration and the exact expressions for these form factors $F_{I}$, $F_{II}$ are put in Appendix A. The form factors exhibit an interesting property, proportional to the product of mass differences $(m_{H_{2}}^{2}-m_{H_{1}}^{2})(m_{H^{\prime}_{2}}^{2}-m_{H^{\prime}_{1}}^{2})$. Thus the neutrino mass can be easily accommodated into the sub-eV order, if either one set of inert scalars is quasi-degenerate without tuning the Yukawa couplings. In particular, if we set $m_{H^{\prime}_{1}}\simeq m_{H^{\prime}_{2}}$, the LFV bound will not be influenced as $H^{\prime}_{1,2}$ do not mediate these processes.

Due to the symmetric property, the Eq. (II.46) can be conveniently recasted into a suitable form for the numerical analysis:

[TABLE]

where the $A$ is an arbitrary anti-symmetric matrix in the order $\lesssim 10^{-9}$ and of complex values if there is CP violation Okada:2015vwh . Therefore after we impose Eq.(II.47), the $y_{\eta}$ coupling is no longer a free parameter but as a function of $y_{S}$ and the neutrino mass form factors. This parameter will be determined by the neutrino oscillation data up to an uncertainty. Notice that $y_{\eta}\lesssim\sqrt{4\pi}$ should be satisfied in the perturbative limit.

II.3 LFV and Muon $g-2$

In this radiative neutrino mass model, the existence of charged scalars and vector-like fermions contribute to lepton flavor violation processes (see Figure. 3), which in turn will severely constrain the Yukawa couplings and masses of heavy scalars and fermions. The relevant Lagrangian for LFV can be expressed as:

[TABLE]

We can calculate the branching ratio for LFV decay process $\ell_{i}\to\ell_{j}\gamma$ in terms of amplitude $a_{L/R}$, which encodes the loop integration of the Feynman diagrams:

[TABLE]

where $G_{\rm F}\approx 1.166\times 10^{-5}$ GeV*-2* is the Fermi constant, $\alpha_{\rm em}(m_{Z})\approx{1/128.9}$ is the fine-structure constant pdg2018 , $C_{21}\approx 1$, $C_{31}\approx 0.1784$, and $C_{32}\approx 0.1736$. In this specific model $a_{R}$ is formulated as:

[TABLE]

where we can see that the loop contributions from two resources (Figure 2.a and Figure 2.b) are in opposite signs. And for the left-handed amplitude, $a_{L}$ is obtained by a mass substitution: $a_{L}=a_{R}(m_{\ell_{i}}\to m_{\ell_{j}})$.

The couplings involved in those LFV processes are $y_{\eta}$ and $y_{S}$, strongly correlated to the neutrino mass matrix. In particular the magnitude of $y_{\eta}$ along with masses $m_{\chi_{1}}$ and $m_{H_{1,2}}$, constrained by the LFV bound, will influence the DM relic density as well. To find out the allowed parameter space for this model, the following upper bounds are imposed TheMEG:2016wtm ; Aubert:2009ag

[TABLE]

where the upper bound from $\mu\to e\gamma$ is the most stringent one with the value in parentheses indicating a future reach of MEG experiment Renga:2018fpd .

The muon anomalous magnetic moment: The muon $g-2$ is a well-measured property and a large $3.6\sigma$ discrepancy of $\Delta a_{\mu}$ between the SM theory and experiment measurement was observed for a long time. For this model, one can estimate the muon $g-2$ through the amplitudes formulated above:

[TABLE]

The deviation from the SM prediction is $\Delta a_{\mu}=a_{\mu}^{\rm exp}-a_{\mu}^{\rm SM}=(2.74\pm 0.73)\times 10^{-9}$ pdg2018 with a positive value. However because our analysis shows the muon $g-2$ is too tiny after imposing other bounds, we just employ the muon $g-2$ as a model quality for reference.

II.4 $Z-Z^{\prime}$ mixing

The effect of the hidden $Z^{\prime}$ at TeV scale will actually decouple from the dark matter physics and we would like to qualify this argument in the section. After the three scalar fields developing VEVs, $U(1)_{H}$ and electroweak symmetries are spontaneously broken so that the mass terms of neutral gauge boson are obtained,

[TABLE]

where $g_{2}$, $g_{1}$ and $g_{H}$ are gauge couplings of $SU(2)_{L}$, $U(1)_{Y}$, and $U(1)_{H}$, respectively. The $Z_{0}$ and $\tilde{Z}$ are the gauge fields for $U(1)_{Y}$ and $U(1)_{H}$ with the $Z_{0}$ mostly composed of the SM Z boson. Here we assume the kinetic mixing between the two Abelian gauge bosons to be negligibly small for simplicity. In case of $x=1$, we parameterise the mass matrix to be:

[TABLE]

where we use the definition of $m_{Z_{0}}=\frac{\sqrt{(g_{1}^{2}+g_{2}^{2})\left(v_{H}^{2}+4v_{\Delta}^{2}\right)}}{2}$, $m_{\tilde{Z}}=3g_{H}v_{\varphi}$, $\epsilon_{1}=\frac{m_{Z_{0}}}{m_{\tilde{Z}}}$ and $\epsilon_{2}=\frac{v_{\Delta}}{v_{\varphi}}$, $\epsilon_{3}=\frac{v_{\Delta}}{\sqrt{v_{H}^{2}+4v_{\Delta}^{2}}}$. The mass matrix can be diagonalized by an orthogonal transformation to be ${\rm Diag}(m^{2}_{Z},m^{2}_{Z^{\prime}})$, and in an approximation of $v_{\Delta}\ll v_{H}\lesssim v_{\varphi}$ and $g_{H}=\mathcal{O}(1)$, this gives:

[TABLE]

If we fix $c_{W}^{2}=g_{2}^{2}/(g_{1}^{2}+g_{2}^{2})$ as the SM value, the $\rho$ parameter can be expressed to be:

[TABLE]

The experimental constraint from the PDG is $\rho_{0,\,\rm{exp}}=1.00039\pm 0.00019$ pdg2018 , which will translate into a requirement of $v_{\Delta}\lesssim 3.5$ GeV. In this paper, we assume the $Z^{\prime}$ boson mass to be above the TeV scale for $v_{\varphi}\gtrsim 350$ GeV. According to Eq. (II.61), this results in a extremely small $|\tan\theta_{Z}|<10^{-5}$ compared with the Yukawa coupling with DM and neutrino. Thus as long as we prefer the DM mass in $\mathcal{O}(100)$ GeV, it will be safe to neglect the the impact of $Z^{\prime}$ on either DM annihilation or DM-nucleon scattering,

II.5 Bound of Electroweak Precision Test

The Electroweak Precision Test (EWPT) on low energy observables can set limits for deviations from the SM. The new physics effects are mainly encoded the oblique parameters $S$, $T$ and $U$, expressed in terms of the transverse part of gauge boson’s self-energy amplitudes. For this model, since the $U$ parameter is suppressed by an additional factor of $v^{2}/M_{\mbox{new}}^{2}$, its effect is neglected. Due to the vector-like nature and degeneracy, the exotic particles of $L^{\prime}_{i}$ have no impact on oblique parameters, i.e. $\Delta S_{f}=\frac{2}{3\pi}(t_{3L}-t_{3R})^{2}=0$ and $\Delta T_{f}=0$ pdg2018 . However the inert scalars $(\eta,s)$ are possible to cause notable deviation to $S=-16\pi\Pi^{\prime}(0)_{W_{3}B}$ and $T=\frac{4\pi}{m_{Z}^{2}s^{2}_{W}c^{2}_{W}}\left[2\Pi_{W_{1}W_{1}}(0)-\Pi_{W_{3}W_{3}}(0)\right]$ Peskin:1991sw , we will discuss their constraints on the mass splitting among $(m_{\eta^{+}},m_{H_{i}})$ and the mixing angle $\sin(\alpha)$. After evaluating the relevant self-energy correlations, we find out the effects from the inert scalars are described by:

[TABLE]

where the loop functions $G(m_{1},m_{2},m_{3})$, $\chi(m_{1},m_{2})$ and $F(m_{1},m_{2})$ are defined as:

[TABLE]

with $\chi(m_{1},m_{2})$ and $F(m_{1},m_{2})$ being symmetric for $m_{1}\leftrightarrow m_{2}$ and vanishing for equal masses, i.e. $\chi(m,m)=F(m,m)=0$. During the calculation, the divergences inherent in the two-point functions are properly cancelled 222For the $T$ parameter, if we calculate it in terms of the gauge boson’s self-energy amplitudes, the divergence is fully captured in the loop function $A_{0}(m^{2})=\frac{1}{i\pi^{2}}\int\frac{d^{4}k}{(k^{2}-m^{2})}$ Haber:2010bw and should be cancelled after counting all the diagrams. The cancellation due to the mixing neutral inert scalars $((\eta_{R}+i\eta_{I})/\sqrt{2}),s)$ (precisely speaking, $(H_{1,2},A_{1,2})$ in the mass basis) demonstrates in the following pattern: $2s^{2}_{\alpha}A_{0}(m_{H_{1}}^{2})+2c^{2}_{\alpha}A_{0}(m_{H_{2}}^{2})-2s^{4}_{\alpha}A_{0}(m_{H_{1}}^{2})-2c^{4}_{\alpha}A_{0}(m_{H_{2}}^{2})-2c^{2}_{\alpha}s^{2}_{\alpha}(A_{0}(m_{H_{1}}^{2})+A_{0}(m_{H_{2}}^{2}))=0$.. In case of the SM Higgs $h$ barely mixing with $(\Delta,\varphi)$, the $\Delta T$ is exactly the wave-function renormalisation of the goldstone bosons $G^{+},G^{0}$ with $(\eta^{\pm},((\eta_{R}+i\eta_{I})/\sqrt{2}),s)$ running inside the loops(referring to Appendix B for detail)  Barbieri:2006dq . While for the $\Delta S$, the function $G(m_{1},m_{2},m_{3})$ is related to $\frac{d}{dp^{2}}\left[B_{22}(p^{2},m_{1}^{2},m_{2}^{2})-B_{22}(p^{2},m_{3}^{2},m_{3}^{2})\right]|_{p^{2}=0}=\frac{1}{2}\int_{0}^{1}dx~{}x(1-x)\ln[\frac{xm_{1}^{2}+(1-x)m_{2}^{2}}{m_{3}^{2}}]$, using the Passarino-Veltman function $B_{22}$ defined in Passarino:1978jh .

The bound for the $S$ and $T$ parameters is obtained from the precision electroweak data, such as $M_{Z}$ and $\Gamma_{Z}$, at the $1\sigma$ deviation by fixing $U=0$ pdg2018 :

[TABLE]

with an off-diagonal correlation of $92\%$. In Figure 4, combining all the discussed parts, we translate the $S$-$T$ $\chi^{2}$ fit into the bounds of inert scalar masses and the mixing angle. Since both $\Delta S$ and $\Delta T$ are symmetric functions of $(m_{H_{1}},m_{H_{2}})$, we focus on the case of $m_{\eta^{+}}=m_{H_{2}}$ for simplicity. The operation of switching $H_{1}\leftrightarrow H_{2}$ is to shift the mixing angle by $\alpha\to(\alpha-\frac{\pi}{4})$. the $1\sigma$ EWPT fit prefers the mass splitting $(m_{H_{2}}-m_{H_{1}})>0$ in a small range of $0.92<\cos\alpha<1.0$ , i.e. $H_{2}(A_{2})$ dominantly composed of $\eta_{R}(\eta_{I})$ should be heavier. However at $3\sigma$ fit, $(m_{H_{2}}-m_{H_{1}})$ is permitted in either sign $(+/-)$ for $0<\alpha<\frac{\pi}{2}$, with its magnitude decreasing with $\cos\alpha$. In the right plot, we show that assuming $m_{H_{2}}>m_{H_{1}}$, the $S$-$T$ bound requires $(m_{H_{2}}-m_{H_{1}})\subset(30,120)$ GeV at $1\sigma$ and $(m_{H_{2}}-m_{H_{1}})<150$ GeV at $3\sigma$ for $100<m_{H_{1}}<500$ GeV under the condition specified in the caption.

III Dark matter

The relic density for a DM specie $X$ is determined by its energy density, $\propto m_{X}n_{X}(T_{0})$ in the present universe, where the number density $n_{X}$ is governed by the Boltzmann equation during the decoupling phase plus the afterwards expansion effect. For a Dirac fermion DM stabilised by a $\mathbb{Z}_{3}$ symmetry, semi-annihilation modes in addition to annihilation are expected to contribute. The Boltzmann equation can be recasted into an evolution in terms of a yield by defining $Y_{X}=n_{X}/s$ with $s$ to be entropy density and $x=M_{X}/T$ where the temperature is scaled by the DM mass. The redefined equation reads:

[TABLE]

where $A,S$ stand for annihilation and semi-annihilation, $H(x=1)$ is the Hubble constant at $T=M_{X}$, $g_{*}$ is the effective total number of relativistic degrees of freedom and $M_{\rm pl}=1.22\times 10^{19}[{\rm GeV}]$ is the Planck mass. The $\frac{1}{2}$ factor in the second term of Eq. (III.1) is due to the identical initial particles 333For the semi-annihilation, considering the evolution of number density for one specie $X$, we need take into account the processes of $XX\to\bar{X}\nu_{i}$ and $\bar{X}\bar{X}\to X\bar{\nu}_{i}$, where the number of the specie $X$ is only depleted by  $1$ in the forward direction, same as in the particle-antiparticle annihilation. Thus the Boltzmann equation with only semi-annihilation mode should be: $\frac{dn_{X}}{dt}+3Hn_{X}=-\frac{1}{2}\langle\sigma v\rangle_{\rm Semi}[n_{X}^{2}-n_{X}n_{X}^{\rm{eq}}]$. This is different from the DM annihilation of Majorana fermions, where the depletion number is $2$, and compensates the phase space factor $\frac{1}{2}$ from identical particles. and $\langle\sigma v_{\rm rel}\rangle$ is the thermal average of velocity weighted cross section which represents the DM interaction rate. This equation can be analytically solved in a proper approximation by matching the results from two regions at the freeze-out point. A brief review for this approach will be presented here in order to clarify the missing $1/2$ in some literature. We will start by defining a quality $\Delta=Y_{X}-Y_{X}^{\rm eq}$, so that the original equation is transformed into:

[TABLE]

where the Maxwell-Boltzmann distribution will be used for the yield in equilibrium so that $Y_{X}^{\rm eq}(x)\propto x^{3/2}e^{-x}$. For $x\ll x_{f}$, we can obtain:

[TABLE]

and for $x\gg x_{f}$, the integration of Boltzmann equation gives:

[TABLE]

Thus the relic density at the present universe is found as:

[TABLE]

where $\Omega h^{2}$ is rescaled by the critical density $\rho_{c}=3H^{2}/8\pi G$. We times a factor $2$ for the relic density in order to count the contribution from the antiparticle $\bar{X}$ and set $g_{*}(x_{f})\approx 100$ at the point of freeze-out. Here $\langle\sigma v\rangle_{\rm A}$ is the thermal average for annihilation, while $\langle\sigma v\rangle_{\rm S}$ is for semi-annihilation. Then the freeze-out temperature $x_{f}$ is determined by the boundary condition $\Delta(x_{f})=c~{}Y_{X}^{\rm eq}(x_{f})$ with $c=\sqrt{2}-1$ to be:

[TABLE]

which is up to a $1/2$ factor for the semi-annihilation part as given by DEramo:2010keq and we set $g=2$ for a fermion DM of two degrees of freedom without counting its antiparticle Kolb .

As we can see that in order to estimate the relic density, one has to calculate the thermal average of cross section times the relative velocity $\langle\sigma v_{rel}\rangle$. Generally the thermal average is approximated by an expansion in order of $x^{-n}$ ($\langle v^{2}\rangle\sim\frac{6}{x}$ in the non-relativistic limit). However in our case, the dominant DM cross section proceeds through an $S$-channel with one very narrow resonance $\Gamma_{M}/M_{X}\ll v_{\rm rel}$ and one wider resonance $\Gamma_{M}/M_{X}\sim v_{\rm rel}$. Also for a $S$-channel interaction mediated by a scalar, the s-wave is vanishing for the velocity averaged cross section, thus the expansion in terms of $v_{\rm rel}^{2}$ is complicated to handle for two resonances interfering with each other. We prefer to use the integration approach to evaluate $\langle\sigma v_{\rm rel}\rangle$ which is given by Gondolo:1990dk ; Edsjo:1997bg

[TABLE]

where $s=(k_{1}+k_{2})^{2}$ is a Mandelstam variable and $K_{1,2}$ are the modified Bessel functions of order 1 and 2 respectively.

[TABLE]

Here $\sigma_{XX}^{i}$ is the cross section of the $2\to 2$ process (denoting $\bf{k}_{1}$ as 3-momentum of the first out-going particle) and with the amplitude squared $|{\cal\bar{M}}_{1,2}|^{2}$ corresponding to $X\bar{X}\to\nu_{i}\bar{\nu}_{j}$ and $X\bar{X}\to\ell_{i}\bar{\ell}_{j}$ in Fig. 5(a-b) and the third $|{\cal\bar{M}}_{3}|^{2}$ standing for $XX\to\bar{X}\nu_{i}$, i.e. the semi-annihilation as depicted in Fig. 5(c)-(e).

We derive the analytic expression for each amplitude squared present in Eq.(III.11). Let us consider the case that only the lightest flavor of $\chi_{i}$ is the DM candidate. By defining $X=\chi_{1}$ and assuming $y^{L}_{\chi}=y^{R}_{\chi}$, the DM-scalar interaction in this model is written as:

[TABLE]

For the annihilation processes, $|{\cal\bar{M}}_{1,2}|$ are the usual amplitude squared with the spin averaged for the initial states and summed for the final states. However a special treatment is needed for $|{\cal\bar{M}}_{3}|$ because of the identical incoming particles. As illustrated in Fig. 5(c)-(e), the semi-annihilation proceeds in $S$, $T$ and $U$ channels after counting the momentum exchanging for the identical initial particles. In particular, there is a symmetry factor $2$ for the $S$-channel amplitude 444We need consider the momentum exchanging for the identical initial particles due to the phase space integration in thermal average. For semi-annihilation $X(p_{1})X(p_{2})+X(p_{2})X(p_{1})\to\bar{X}(k_{1})v_{i}(k_{2})$, the $S$-channel amplitude is proportional to $[\bar{u^{c}}(p_{1})u(p_{2})-\bar{u^{c}}(p_{2})u(p_{1})][\bar{v}(k_{1})u(k_{2})]=2~{}[\bar{v}(p_{1})u(p_{2})]~{}[\bar{v}(k_{1})u(k_{2})]$, where we use the identities $u^{c}=C\bar{u}^{T}=v$ and $\bar{v}(p_{2})u(p_{1})=u^{T}(p_{1})C^{-1}C\bar{v}^{T}(p_{2})=-\bar{v}(p_{1})u(p_{2})$, with $C=i\gamma^{0}\gamma^{2}$ being the charge conjugate operator. This is similar to the identical scalar case $\phi\phi\to H_{a}$, the symmetry factor is normally encoded in the vertex.. Combining all channels, we can arrive the following analytic expressions:

[TABLE]

[TABLE]

In the ${\cal M}_{3}$ amplitude of semi-annihilation, we define $S_{inv}^{a}=1/(s-m_{a}^{2}+im_{a}\Gamma_{a})$, $T_{inv}^{a}=1/(2M_{X}^{2}-m_{a}^{2}-2p_{1}\cdot k_{1})$, $U_{inv}^{a}=1/(M_{X}^{2}-m_{a}^{2}-2p_{1}\cdot k_{2})$ and the index $a=1,2$ corresponds to $H_{1}(A_{1})$, $H_{2}(A_{2})$ respectively. The inner products are given in Appendix C.

For the $S$-channel amplitude, the widths of inert scalars $H_{1,2}(A_{1,2})$ enter into the Breit-Wigner propagator $S_{inv}^{a}$, whose magnitude near two on-shell poles $m_{H_{1}}=2M_{X}$ or $m_{H_{2}}=2M_{X}$ is determined by the $\Gamma_{H_{1}}$ or $\Gamma_{H_{2}}$. Under this consideration we will only be interested in the parameter region $m_{H_{1}}<m_{H_{2}}<{\rm min}(m_{E^{\prime}_{i}},m_{\chi_{2}},m_{\chi_{3}})$ to ensure a narrow resonance. Therefore the decay widths of $\Gamma_{H_{1}}$ ($=\Gamma_{A_{1}}$) and $\Gamma_{H_{2}}$ ($=\Gamma_{A_{2}}$) are formulated as:

[TABLE]

and for $H_{2}$, one more decay channel $H_{2}\to H_{1}h_{0}$, with a coupling vertex of $\frac{1}{2}\lambda_{0}v_{\varphi}(c_{\alpha}^{2}-s_{\alpha}^{2})=s_{\alpha}c_{\alpha}(c_{\alpha}^{2}-s_{\alpha}^{2})(m_{H_{1}}^{2}-m_{H_{2}}^{2})/v_{H}$ and $h_{0}$ being the SM Higgs boson, will be open if it is permitted by kinematics.

[TABLE]

where the step function is defined as $\theta(x)=1$ only for $x>0$ otherwise being zero.

III.1 Relic density analysis

In this section, we will show the numerical analysis to satisfy all the constraints discussed in Section II. We find out that after imposing the LFV bounds and neutrino oscillation data, one DM-neutrino-scalar coupling $y_{\eta_{i1}}$ populates in the range of $(10^{-3},1.0)$, so that the annihilation process in this model can not account for a correct relic density. However a large enhancement for $\langle\sigma v\rangle$ could be achieved if the semi-annihilation proceeds through a $S$-channel and in the vicinity of one narrow-width resonance. Since Eq. (III.15) indicates two resonances of complex scalars are deconstructive, one condition $100<(m_{H_{2}}-m_{H_{1}})<150$ GeV is imposed in the analysis, with the upper limit from the $3\sigma$ EWPT fit at $\cos\alpha=\frac{1}{\sqrt{2}}$. Thus for a given DM mass, only one resonance can effectively be on-shell. On the other hand, we will require $m_{H^{\prime}_{1}}\simeq m_{H^{\prime}_{2}}$, i.e. quasi-degenerate, in order to satisfy the neutrino oscillation data. This condition can be easily fulfilled if we set the mixing term $\lambda^{\prime}_{0}H^{\dagger}\eta^{\prime}s^{\prime*}\varphi^{*}$ to be tiny. In order to simplify the analysis, we adopt several assumptions as below:

[TABLE]

We set $m_{\eta^{\pm}}=m_{H_{2}}$ which is consistent with the EWPT bound as shown in Figure 4 and $\ y^{\prime}_{\eta^{\prime}},\ y^{\prime}_{s^{\prime}}$ are taken to be diagonal matrices. Under these assumptions, a numerical scan is conducted for the parameter space by imposing the relevant neutrino and LFV bounds and limiting the relic density to be $0.117<\Omega h^{2}<0.123$. We explore the two on-shell scenarios in two overlapping DM mass regions with $m_{H_{1}}=2M_{X}$ for $80<M_{X}<350$ GeV and $m_{H_{2}}=2M_{X}$ for $200<M_{X}<400$ GeV. Furthermore, in order to work well under the Breit-Wigner narrow width prescription, we remove the points with $\rm{max}(\frac{\Gamma_{H_{1}}}{m_{H_{1}}},\frac{\Gamma_{H_{2}}}{m_{H_{2}}})>0.2$. For the latter case of $m_{H_{2}}=2M_{X}$, we will impose a smaller splitting $(m_{H_{2}}-m_{H_{1}})\simeq(115,125)$ GeV, thus $m_{H_{1}}\gg m_{X}$. This condition will ensure $\Gamma_{H_{1}}\ll m_{H_{1}}$ and avoid co-annihilation from scalars. From the left plot in Fig: 6 we can see that the observed relic density dominantly comes from the semi-annihilation. At the time of freeze out $x_{f}\approx 21.0$ (calculated by Eq.(III.8)), the thermal average of cross section is within the range of $5.98\times 10^{-10}~{}{\rm GeV^{-2}}\lesssim\langle\sigma v_{Semi}\rangle\lesssim 8.83\times 10^{-10}~{}{\rm GeV^{-2}}$, where the larger value normally corresponds to a larger DM mass. In the right plot we show the allowed region in the $(M_{X},y_{\chi_{11}})$ plan with other parameters randomly scanned. The plot demonstrates that a small DM mass $M_{X}<200$ GeV is more sensitive to the lighter $H_{1}+iA_{1}$ resonance and permits a DM Yukawa coupling $y_{\chi_{11}}\gtrsim 0.1$. However for $M_{X}>200$ GeV, our fitting analysis indicates a larger DM coupling $y_{\chi_{11}}\gtrsim 0.5$, which is close to the perturbative limit $\sqrt{4\pi}$ regardless of the lighter or heavier resonance scenario.

Fig. 7 presents the mass ranges for $M_{E^{\prime}}(=M_{N^{\prime}})$ and $M_{N}$ which enter into the numerator of neutrino mass form factors as well as values of $|\Delta a_{\mu}|$ versus $M_{X}$. The typical value for the lightest vector-like fermions $L^{\prime}_{i}$ lies in $0.5-2.5$ TeV, but the degeneracy results in no effect on EWPT. Also this mass range of $M_{E^{\prime}}$ is not sensitive to the LHC bound for charged lepton pairs plus missing transverse energy Cai:2018upp . While after enforcing all the bounds, the maximum value for $|\Delta a_{\mu}|$ is of order $\lesssim 10^{-14}$, even lower for most benchmark points, is negligible compared with the $3.6~{}\sigma$ deviation of order $10^{-9}$ as measured by the experiment. Thus this model can not simultaneously account for the large discrepancy in muon $g-2$.

Direct detection: In our case, there are no direction interactions between $H_{1,2}/A_{1,2}$ and quarks at the tree level, therefore the constraints of direct detection searches should be satisfied without difficulty.

IV Conclusions and discussions

We have constructed a neutrino mass model based on hidden local $U(1)_{H}$ symmetry which gives rise to a Dirac fermion type of Dark matter. The neutrino masses are generated at the two-loop level due to the symmetry and particle content. Furthermore because the form factor of the neutrino mass is proportional to the mass squared differences of inert scalars, we require one set of inert scalars to be quasi-degenerate so that a sub-eV scale neutrino mass can be achieved without large fine-tuning for the Yukawa couplings. As variation to this model, we illustrate that the heavy $Z^{\prime}$ associated with the $U(1)$ will not impact the DM annihilation because its mixing with SM $Z$ boson is induced by a complex triplet field $\Delta$, whose VEV is severely constrained by $\rho$-parameter. Particularly, the presence of inert scalars $(\eta,s)$ gives rise to notable $S$ and $T$ deviations. Note that the impact of singlet $s$ on oblique parameters is via the mixing with doublet $\eta$. The $3\sigma$ EWPT fit prefers the mass splitting of $|m_{H_{2}}-m_{H_{1}}|\lesssim 150$ GeV provided $\cos\alpha=\frac{1}{\sqrt{2}}$ and $m_{\eta}^{+}=m_{H_{2}}$.

Our DM is is the lightest neutral particle stabilised by a discrete $\mathbb{Z}_{3}$ parity which is a residual symmetry of $U(1)_{H}$ after spontaneous symmetry breaking. Therefore, in addition to the standard DM annihilation process, DM semi-annihilation is induced in this model. After imposing the LFV bounds and neutrino oscillation data and assuming no specific flavour structure in Yukawa couplings, we find out that the $S-$channel semi-annihilation plays an important role to determine the observed relic density with a DM mass of $\mathcal{O}(100)$ GeV. Our analysis demonstrates that the lighter and heavier resonances can contribute significantly when either one is actually put on-shell and the allowed DM-scalar Yukawa coupling is in the range of ($0.1$-$\sqrt{4\pi}$) depending on the DM mass region.

Acknowledgments

The research of H.C. is supported by the Ministry of Science, ICT and Future Planning of Korea, the Pohang City Government, and the Gyeongsangbuk-do Provincial Government.

Appendix A Loop functions for neutrino mass

The neutrino mass in this radiative seesaw model is generated by the two-loop Feynman diagrams in Figure 2. It is convenient to decompose the mass matrix as $(m_{\nu})_{ij}=m_{\nu_{ij}}^{(I)}+m_{\nu_{ij}}^{(II)}+[m_{\nu_{ij}}^{(I)}]^{T}+[m_{\nu_{ij}}^{(II)}]^{T}$, with $m_{\nu_{ij}}^{(I)}$ and $m_{\nu_{ij}}^{(II)}$ calculated to be:

[TABLE]

[TABLE]

For clarity, we can redefine $m_{\nu_{ij}}^{(I/II)}=\frac{1}{(4\pi)^{4}}~{}y_{\eta_{ia}}F_{I/II}(H_{1,2},H^{\prime}_{1,2})_{ab}~{}y_{bj}^{T}$ by extracting out a loop factor and Yukawa couplings in the outer loop of Feynman diagrams. After imposing the Feynman parametrisation, the two loop functions $F_{I/II}(H_{1,2},H^{\prime}_{1,2})$ are given by:

[TABLE]

where we use the definitions: $[da]_{3}\equiv\int_{0}^{1}db\int_{0}^{1-b}dc$ with $a=1-b-c$, and $[d\alpha]_{5}\equiv\int_{0}^{1}d\alpha\int_{0}^{1-\alpha}d\beta\int_{0}^{1-\alpha-\beta}d\gamma\int_{0}^{1-\alpha-\beta-\gamma}d\rho$ with $\sigma=1-\alpha-\beta-\gamma-\rho$. Note that these form factors are finite and will be numerically evaluated.

Appendix B T parameter from mixing inert scalars

Since the longitude modes of $W,Z$ gauge bosons are $\partial_{\mu}G^{\pm,0}$, the $T$ parameter is easily calculated from the wave-function renormalisation of Goldstone bosons. We are going to show that two approaches are matching with each other. The relevant terms from the scalar potential are:

[TABLE]

Due to the $\mathbb{Z}^{3}$ parity, there is no mass splitting among the imaginary and real parts of inert neutral scalars. The masses can be read off from Eq.(B.1):

[TABLE]

with the diagonal parts to be

[TABLE]

The following identities will hold for the mass eigenstates and rotating angle:

[TABLE]

since $\delta\rho=\delta Z_{G^{+}}-\delta Z_{G^{0}}$, the two-point self-energy diagrams in Fig. 8 give us:

[TABLE]

with the function $f(m_{1},m_{2})=-i\frac{d\Pi(p^{2})}{dp^{2}}|_{p^{2}=0}$, and $\hat{\alpha}=\frac{e^{2}}{4\pi}$. The $\Pi(p^{2})$ is defined as:

[TABLE]

with $\Delta=-p^{2}(1-x)x+xm_{1}^{2}+(1-x)m_{2}^{2}$. Thus we can obtain the analytic formula:

[TABLE]

Using Eqs.(B.2), (B.6), (B.7), the coefficients in Eq.(B.8) are related to be:

[TABLE]

Then after substituting those identities back to Eq.(B.8), we obtain the $\Delta T$ expression in Eq.(II.64).

Appendix C Inner products for the amplitudes

[TABLE]

where $s\equiv(p_{1}+p_{2})^{2}$ is a Mandelstam valuable, $m_{1,2}(p_{1,2})$ are initial state masses(momenta), while $n_{1,2}(k_{1,2})$ are final state masses(momenta).