#  Three loop neutrino model with isolated $k^{\pm\pm}$

Kenji Nishiwaki<EMAIL_ADDRESS>School of Physics, Korea Institute for
Advanced Study,
Seoul 130-722, Republic of Korea  Hiroshi Okada<EMAIL_ADDRESS>School of Physics, Korea Institute for Advanced Study,
Seoul 130-722, Republic of Korea  Physics Division, National Center for
Theoretical Sciences, National Tsing-Hua University, Hsinchu 30013, Taiwan
Yuta Orikasa<EMAIL_ADDRESS>School of Physics, Korea Institute for
Advanced Study,
Seoul 130-722, Republic of Korea  Department of Physics and Astronomy, Seoul
National University, Seoul 151-742, Republic of Korea

###### Abstract

We propose a three loop radiative neutrino mass scenario with an isolated
doubly charged singlet scalar $k^{\pm\pm}$ without couplings to the charged
leptons, while two other singly charged scalars $h_{1}^{\pm}$ and
$h_{2}^{\pm}$ attach to them. In this setup, the lepton flavor violation
originating from $k^{\pm\pm}$ exchanges is suppressed and the model is less
constrained, where some couplings can take sizable values. As reported in our
previous work Hatanaka:2014tba , the loop suppression factor at the three loop
level would be too strong and realized neutrino masses in a three loop
scenario could be smaller than the observed minuscule values. The sizable
couplings can help us to enhance neutrino masses without drastically large
scalar trilinear couplings appearing in a neutrino mass matrix, which tends to
drive the vacuum stability to become jeopardized at the one loop level. Now
the doubly charged scalar $k^{\pm\pm}$ has less constraint via lepton flavor
violation and the vacuum can be quite stable, and thus a few hundred GeV mass
in $k^{\pm\pm}$ is possible, which is within the LHC reach and this model can
be tested in the near future. Note that the other $h_{1}^{\pm}$ and
$h_{2}^{\pm}$ should be heavy at least around a few TeV. We suitably arrange
the charges of an additional global $U(1)$ symmetry, where the decay constant
of the associated Nambu-Goldstone boson can be around a TeV scale
consistently. Also, this model is indirectly limited through a global analysis
on results of the LHC Higgs search and issues on a dark matter candidate, the
lightest Majorana neutrino. After $h_{1}^{\pm}$ and $h_{2}^{\pm}$ are
decoupled, this particle couples to the standard model particles only through
two charge parity even scalars in theory and thus information on this scalar
sector is important. Consistent solutions are found, but a part of them is now
on the edge.

Neutrinos, Dark Matter

††preprint: KIAS-P15035

##  I Introduction

Recently, the second round of the physics run at the Large Hadron Collider
(LHC) started and a magnificent operation for exploring the scalar sector
describing the electroweak (EW) scale was launched. The greatest achievement
at the first run of the LHC is the observation of the new Higgs-like scalar
boson, which was the last missing piece of the standard model (SM), around
$125\,\text{GeV}$ by the ATLAS and the CMS experiments Aad:2012tfa ;
Chatrchyan:2012ufa .

The Yukawa couplings to heavy fermions of the SM, namely the top and bottom
quarks and the tau lepton have been surveyed with good precision, whereas the
lighter states are still mysterious from both the experimental and theoretical
point of view. The extremely tiny observed masses of active components of the
neutrinos would be key to investigating the scalar sector theoretically
because we should accept at least $10^{12}$-order hierarchy in the Higgs
Yukawa couplings for describing the neutrino nature within the SM.

One of the most stimulating ideas for addressing this issue is radiative
generations of the neutrino profiles. Loop suppression factors should appear
in the neutrino masses in this type of scenario, which help to alleviate the
hierarchy in couplings. Another motivation for this direction is that the
continuous and/or discrete symmetries ensuring the loop origin could also
guarantee the existence of a (or multiple) dark matter (DM) candidate(s).
Following the landmarks Zee ; Li ; zee-babu ; Krauss:2002px ; Ma:2006km ,
recently a variety of works on radiative seesaw model have been done
Hambye:2006zn –Kashiwase:2015pra , where we can also find studies emphasizing
non-Abelian discrete symmetry Ahn:2012cg ; Ma:2012ez ; Kajiyama:2013lja ;
Hernandez:2013dta ; Ma:2014eka ; Aoki:2014cja ; Ma:2014yka ; Ma:2015pma ,
radiative generation of quark/charged lepton masses Hernandez:2013hea ;
Ma:2013mga ; radlepton1 ; radlepton2 ; Okada:2014nsa , operator analysis
Bonnet:2012kz ; Sierra:2014rxa , radiative models accompanying conformal EW
symmetry breaking Davoudiasl:2014pya ; Lindner:2014oea ; Okada:2014nea , and
others Babu:2002uu ; AristizabalSierra:2006ri ; AristizabalSierra:2006gb ;
Nebot:2007bc ; Bouchand:2012dx ; Kajiyama:2013sza ; McDonald:2013hsa ;
Ma:2014cfa ; Schmidt:2014zoa ; Herrero-Garcia:2014hfa ; Ahriche:2014xra ;
Long1 ; Long2 .

From the naturalness point of view, higher-loop generation is better. The
first three loop model for a natural explanation of the neutrino profiles was
proposed in Ref. Krauss:2002px  and the following works continue Aoki:2008av ;
Gustafsson ; Kajiyama:2013lja ; Ahriche:2014xra ; Ahriche:2014cda ;
Ahriche:2014oda ; Chen:2014ska ; Okada:2014oda ; Hatanaka:2014tba ;
Jin:2015cla ; Culjak:2015qja ; Geng:2015coa ; Ahriche:2015wha , where in such
situations, couplings related with the neutrino masses can be close to unity
compared with those in models with one or two loop level generation.

On the other hand, three loop generations could face problems owing to the
largeness of couplings. As discussed in Ref. Hatanaka:2014tba , the three loop
suppression factor $1/(4\pi)^{6}$ sometimes looks very strong and we probably
enhance part of couplings for a suitable realization of the neutrino masses.
Considering large Yukawa-type couplings with lepton flavor violation (LFV)
tend to result in unacceptable enhancements in LFV processes. Hence, choosing
sizable scalar trilinear couplings appearing in the neutrino mass matrix,
which do not generate LFV directory, seems to be a reasonable prescription.
But this diagnosis could still be at least partly a misjudgment. When we go
for one loop level, as shown in Ref. Hatanaka:2014tba , these substantial
trilinear couplings give drastic negative contributions to quartic couplings
of charged scalar bosons and the vacuum can be threatened with being
destabilized. The last option would be to accept decouplings of additional
charged bosons, namely around $10\,\text{TeV}$, where the neutrino profile
itself can be suitably generated. But finding clear signals at colliders
becomes very difficult even at the updated $13$ or $14\,\text{TeV}$ LHC.

To circumvent the situation, we propose a refined three loop model by use of a
global $U(1)$ symmetry without additional discrete symmetry. A key point is
introducing additional Majorana neutrinos, which can violate lepton flavors in
the fermion line inside diagrams describing the neutrino masses. It is
important that with the above setup, a doubly charged scalar $k^{\pm\pm}$ no
longer needs to have direct couplings to the charged leptons. We mention that
in our previous model without Majorana neutrinos in Ref. Hatanaka:2014tba ,
like in the Zee-Babu model zee-babu , a doubly charged scalar should attach
the charged leptons to generate violation in lepton flavors, where tree-level
lepton flavor violating processes are generated and we cannot put large values
in the corresponding couplings consistently.

In the present $k^{\pm\pm}$-isolated scenario, the doubly charged scalar
$k^{\pm\pm}$ is quarantined from the charged lepton sector at the leading
order by a suitable choice of global $U(1)$ charges. Now, $k^{\pm\pm}$ cannot
contribute to phenomena with LFV at the leading order and constraints on
couplings are weakened. Consequently, the scalar trilinear couplings can take
smaller values and the vacuum stability would not be so serious even when we
consider a few hundred GeV $k^{\pm\pm}$, the detection of which can be an
evident signal for probing this scenario at the LHC experiments.

This paper is organized as follows. In Sec. II, first we introduce our setups
and subsequently, we discuss miscellaneous issues in this model, namely, forms
of scalar masses and mixings, properties of the Nambu-Goldstone (NG) boson
associated with the breakdown of the global $U(1)$ symmetry, sizable
correction via charged scalars to vacuum stability, form of the active
neutrino mass matrix at the three loop level, and details on processes
accompanying LFV in order. In the following Sec. III, after having a
discussion on analogies with the Zee-Babu model zee-babu , we execute
parameter scans both in the normal and inverted hierarchies. In Sec. IV, we
make global fits of signal strengths of the Higgs production in various
channels announced by the ATLAS and the CMS experiments, which restrict
possible values of the mass of the doubly charged scalar and the mixing angle
between the SM Higgs boson and an additional charge parity ($CP$) even
scalars, whose vacuum expectation value (VEV) breaks the global $U(1)$
symmetry. In Sec. V, we discuss properties of the dark matter candidate of
this scenario, which is the lightest right-hand neutrino, through a relic
density calculation, an estimation of the invisible decay width of the
observed $125\,\text{GeV}$ scalar and an evaluation of spin-independent direct
detection cross section. Section VI is devoted to summarizing results and
making conclusions. In Appendix A, we give analytic forms of the loop
functions describing lepton-flavor-violating processes. In Appendix B, a part
of partial decay widths of the two $CP$ even scalars with nontrivial forms is
described. In Appendix C, we summarize the averaged matrix elements squared
for relic density calculation.

##  II Basic issues on the scenario

###  II.1 Model setup

We discuss a three loop induced radiative neutrino model with a $U(1)$ global
symmetry. The particle contents and their charges are shown in Table 1. We
introduce three Majorana fermions $N_{R_{1,2,3}}$ and new bosons; one gauge-
singlet neutral boson $\Sigma_{0}$, two singly charged singlet scalars
($h^{\pm}_{1},h^{\pm}_{2}$), and one gauge-singlet doubly charged boson
$k^{\pm\pm}$ to the SM. We assume that only the SM-like Higgs $\Phi$ and the
additional neutral scalar $\Sigma_{0}$ have VEVs, which are symbolized as
$\langle\Phi\rangle\equiv v/\sqrt{2}$ and $\langle\Sigma_{0}\rangle\equiv
v^{\prime}/\sqrt{2}$, respectively. Here, we set $v$ as $\simeq
246\,\text{GeV}$. It is natural that $v^{\prime}$ is greater than $v$ to some
extent to evade a large mixing, which is constrained by the LHC data; see Sec.
IV for details. $x\,(\neq 0)$ is an arbitrary number of the charge of the
hidden $U(1)$ symmetry,111Here, we assume that this $U(1)$ symmetry is global.
However, one can straightforwardly move to the local one by introducing
additional fermions and under the assignments, neutrino mass matrix is
generated at the three loop level. Note that a remnant $Z_{2}$ symmetry
remains after the hidden $U(1)$ symmetry breaking and the particles
$N_{R_{1,2,3}}$ and $h^{\pm}_{2}$ have negative parities. Then, when a
Majorana neutrino is the lightest among them, the stability is accidentally
ensured.

| Lepton Fields | Scalar Fields
---|---|---
Characters |  $L_{L_{i}}$ |  $e_{R_{i}}$ |  $N_{R_{i}}$ |  $\Phi$ |  $\Sigma_{0}$ |  $h^{+}_{1}$ |  $h^{+}_{2}$ |  $k^{++}$
$SU(3)_{C}$ | $\bm{1}$ | $\bm{1}$ | $\bm{1}$ | $\bm{1}$ | $\bm{1}$ | $\bm{1}$ | $\bm{1}$ | $\bm{1}$
$SU(2)_{L}$ | $\bm{2}$ | $\bm{1}$ | $\bm{1}$ | $\bm{2}$ | $\bm{1}$ | $\bm{1}$ | $\bm{1}$ | $\bm{1}$
$U(1)_{Y}$ | $-1/2$ | $-1$ | $0$ | $1/2$ | $0$ | $1$ | $1$ | $2$
$U(1)$ | $0$ | $0$ | $-x$ | $0$ | $2x$ | $0$ | $x$ | $2x$
Table 1: Contents of lepton and scalar fields and their charge assignment
under $SU(3)_{C}\times SU(2)_{L}\times U(1)_{Y}\times{U(1)}$, where $U(1)$ is
an additional global symmetry and $x\neq 0$. The subscripts found in the
lepton fields $i\,(=1,2,3)$ indicate generations of the fields. The bold
letters emphasize that these numbers correspond to representations of the Lie
groups of the NonAbelian gauge interactions.

The relevant Lagrangian for Yukawa sector $\mathcal{L}_{Y}$ and scalar
potential $\mathcal{V}$ allowed under the global symmetry is given as

$\displaystyle{-}\mathcal{L}_{Y}$
$\displaystyle=(y_{\ell})_{ij}\bar{L}_{L_{i}}\Phi
e_{R_{j}}+\frac{1}{2}(y_{L})_{ij}\bar{L}^{c}_{L_{i}}L_{L_{j}}h^{+}_{1}+(y_{R})_{ij}\bar{N}_{R_{i}}e^{c}_{R_{j}}h_{2}^{-}+\frac{1}{2}(y_{N})_{ij}\Sigma_{0}\bar{N}_{R_{i}}^{c}N_{R_{j}}+\rm{h.c.}{,}$
(II.1) $\displaystyle\mathcal{V}$
$\displaystyle=m_{\Phi}^{2}|\Phi|^{2}+m_{\Sigma}^{2}|\Sigma_{0}|^{2}+m_{h_{1}}^{2}|h^{+}_{1}|^{2}+m_{h_{2}}^{2}|h_{2}^{+}|^{2}+m_{k}^{2}|k^{++}|^{2}$
$\displaystyle+\Bigl{[}\lambda_{11}\Sigma_{0}^{*}h^{-}_{1}h^{-}_{1}k^{++}+\mu_{22}h^{+}_{2}h^{+}_{2}k^{--}+{\rm
h.c.}\Bigr{]}+\lambda_{\Phi}|\Phi|^{4}+\lambda_{\Phi\Sigma}|\Phi|^{2}|\Sigma_{0}|^{2}+\lambda_{\Phi
h_{1}}|\Phi|^{2}|h^{+}_{1}|^{2}$ $\displaystyle+\lambda_{\Phi
h_{2}}|\Phi|^{2}|h^{+}_{2}|^{2}+\lambda_{\Phi
k}|\Phi|^{2}|k^{++}|^{2}+\lambda_{\Sigma}|\Sigma_{0}|^{4}+\lambda_{\Sigma
h_{1}}|\Sigma_{0}|^{2}|h^{+}_{1}|^{2}+\lambda_{\Sigma
h_{2}}|\Sigma_{0}|^{2}|h^{+}_{2}|^{2}$ $\displaystyle+\lambda_{\Sigma
k}|\Sigma_{0}|^{2}|k^{++}|^{2}+\lambda_{h_{1}}|h_{1}^{+}|^{4}{+}\lambda_{h_{1}h_{2}}|h_{1}^{+}|^{2}|h^{+}_{2}|^{2}+\lambda_{h_{1}k}|h_{1}^{+}|^{2}|k^{++}|^{2}$
$\displaystyle+\lambda_{h_{2}}|h_{2}^{+}|^{4}+\lambda_{h_{2}k}|h_{2}|^{2}|k^{++}|^{2}+\lambda_{k}|k^{++}|^{4},$
(II.2)

where the indices $i,j$ indicate matter generations and the superscript $``c"$
means charge conjugation.222For $SU(2)_{L}$ doublets, charge conjugation is
defined with the $SU(2)_{L}$ rotation described by a Pauli matrix as
$i\sigma_{2}$. $y_{L}$, $y_{R}$ and $y_{N}$ are antisymmetric, general,
symmetric three-by-three matrices, respectively. The first term of
$\mathcal{L}_{Y}$ generates the charged-lepton masses following the SM manner.
Majorana mass terms are derived from the fourth one after $\Sigma_{0}$ obtains
a VEV. Note that this VEV also generates an effective trilinear coupling
$\mu_{11}\left({\text{in front of }}h_{1}^{-}h_{1}^{-}k^{++}\right)$ from the
quartic coupling $\lambda_{11}\left({\text{in front of
}}\Sigma_{0}^{\ast}h_{1}^{-}h_{1}^{-}k^{++}\right)$, where the coefficient is
given as

$\displaystyle\mu_{11}=\lambda_{11}\langle\Sigma_{0}^{\ast}\rangle=\lambda_{11}{\frac{v^{\prime}}{\sqrt{2}},}$
(II.3)

where we use the parametrization of $\langle\Sigma_{0}\rangle$ declared in Eq.
(II.7) in the next subsection. We assume the following two things: (i)
$\lambda_{11}$ and $\mu_{22}$ are positive real; (ii) $y_{N}$ is diagonal and
obeys the hierarchy $(y_{N})_{11}<(y_{N})_{22}<(y_{N})_{33}$ among positive-
real parameters, which means that a generated Majorana mass matrix for $N_{R}$
is also diagonal one and the mass ordering is $M_{N_{1}}<M_{N_{2}}<M_{N_{3}}$.
The concrete forms of the masses are

$\displaystyle M_{N_{1}}=\frac{v^{\prime}}{\sqrt{2}}(y_{N})_{11},\quad
M_{N_{2}}=\frac{v^{\prime}}{\sqrt{2}}(y_{N})_{22},\quad
M_{N_{3}}=\frac{v^{\prime}}{\sqrt{2}}(y_{N})_{33}.$ (II.4)

###  II.2 Mass eigenvalues and eigenstates of scalars

The neutral scalar fields are parametrized as

$\displaystyle\Phi=\left[\begin{array}[]{c}w^{+}\\\
\frac{v+\phi+iz}{\sqrt{2}}\end{array}\right],\quad\
\Sigma_{0}=\frac{v^{\prime}+\sigma}{\sqrt{2}}e^{iG/v^{\prime}}.$ (II.7)

where $v\simeq 246\,\text{GeV}$ is the VEV of the Higgs doublet field, and
$w^{\pm}$ and $z$ are (would-be) NG bosons that are absorbed as the
longitudinal components of the $W$ and $Z$ bosons, respectively. Requiring the
tadpole conditions, $\partial\mathcal{V}/\partial\phi|_{v}=0$ and
$\partial\mathcal{V}/\partial\sigma|_{v^{\prime}}=0$, the resultant mass
matrix squared of the $CP$ even components $(\phi,\sigma)$ is given by

$m^{2}(\phi,\sigma)=\left[\begin{array}[]{cc}2\lambda_{\Phi}v^{2}&\lambda_{\Phi\Sigma}vv^{\prime}\\\
\lambda_{\Phi\Sigma}vv^{\prime}&2\lambda_{\Sigma}v^{\prime 2}\\\
\end{array}\right]=\left[\begin{array}[]{cc}\cos\alpha&\sin\alpha\\\
-\sin\alpha&\cos\alpha\end{array}\right]\left[\begin{array}[]{cc}m^{2}_{h}&0\\\
0&m^{2}_{H}\end{array}\right]\left[\begin{array}[]{cc}\cos\alpha&-\sin\alpha\\\
\sin\alpha&\cos\alpha\end{array}\right],$ (II.8)

where $h$ is the SM-like Higgs ($m_{h}=125\,\text{GeV}$) and $H$ is an
additional $CP$ even Higgs mass eigenstate. The mixing angle $\alpha$ is
determined as

$\sin 2\alpha=\frac{2\lambda_{\Phi\Sigma}vv^{\prime}}{{m^{2}_{H}-m_{h}^{2}}}.$
(II.9)

The neutral bosons $\phi$ and $\sigma$ are rewritten in terms of the mass
eigenstates $h$ and $H$ as

$\displaystyle\phi=h\cos\alpha+H\sin\alpha,\quad\sigma=-h\sin\alpha+H\cos\alpha.$
(II.10)

A NG boson $G$ emerges due to the spontaneous symmetry breaking of the global
$U(1)$ symmetry. The mass eigenvalues for the singly charged bosons
$h_{1}^{\pm}$, $h_{2}^{\pm}$ and the doubly charged boson $k^{\pm\pm}$ are
given as

$\displaystyle m^{2}_{h^{\pm}_{1}}=m_{h_{1}}^{2}+\frac{1}{2}(\lambda_{\Phi
h_{1}}v^{2}+\lambda_{\Sigma h_{1}}v^{\prime 2}),\quad
m^{2}_{h^{\pm}_{2}}=m_{h_{2}}^{2}+\frac{1}{2}(\lambda_{\Phi
h_{2}}v^{2}+\lambda_{\Sigma h_{2}}v^{\prime 2}),$ $\displaystyle
m^{2}_{k^{\pm\pm}}=m_{k}^{2}+\frac{1}{2}(\lambda_{\Phi k}v^{2}+\lambda_{\Sigma
k}v^{\prime 2}),$ (II.11)

where the three charged particles are not mixed due to the symmetries of the
system and thus they themselves are mass eigenstates.

###  II.3 Issues on the Goldstone boson

Accompanying the spontaneous breakdown of a $U(1)$ global symmetry, a NG boson
emerges as an almost massless state in theory, which could play significant
roles in particle physics and cosmology Weinberg:2013kea . Like the usual
Majoron case majoron , our NG boson $G$ communicates with the lepton sector.

An important characteristic of $G$ is that, as described in Table 1, the
lepton doublets and the charged lepton singlets do not hold nonzero charges of
the global $U(1)$. This means that no anomaly-induced interaction to two
photons is generated in our setup, which puts a significant constraint on the
decay constant of NG bosons Cadamuro:2011fd . Thus, in the present scenario,
we can choose “lower” values around a TeV scale without doing any harm.

Another route for constraining models via NG bosons is through the active-
sterile component mixing as through the lepton-flavor-violating transition
like $\mu^{-}\to e^{-}\,\text{(NG)}$ seen in Majoron seesaw scenarios, e.g.,
discussed in Ref. Pilaftsis:1993af . Different from such a situation in our
case, the active and the sterile components cannot mix with each other since
this mixture is prohibited by the residual $Z_{2}$ symmetry after the global
$U(1)$ breaking shown in Table 1. Then the absence of this type of constraint
is assured via the accidentally remaining symmetry. The neutrinoless double
beta decay via $W$ exchanges does not restrict our scenario since the sequence
with $W$ boson also requires the active-sterile mixing. Note that the three
additional charged singlet scalars have no direct coupling to the quarks and
are therefore ineffective.

In contrast, the NG boson $G$ couples to the corresponding $CP$ even component
$\sigma$, which should mix with the Higgs component of the doublet $\Phi$.
This means that $G$ can contribute to physics associated with the $CP$ even
scalars. As we see in Sec. V, the pair annihilation process of the dark matter
candidate $N_{R_{1}}$ is just an example.333Another interesting topics is
collider searches for a NG boson through invisible channels (subsequent decays
from $CP$ even scalars) Cheung:2013oya .

Finally, we briefly comment on possible bounds from cosmological issues. For
example, an effect on cosmic microwave background via cosmic string generated
by the spontaneous breakdown of the global $U(1)$ symmetry possibly put a
constraint on our scenario. The bound discussed in Ref. Battye:2010xz  can be
interpreted as $v^{\prime}<10^{15}\,\text{GeV}$, and thus this issue is
harmless. On the other hand, as we discuss later in Sec. III.2, at least part
of scalar self couplings tends to be $\mathcal{O}(1)$ (at around the EW scale)
owing to the requirements via coexistence of the observed active neutrino
profiles and the null observation in lepton-flavor-violating currents. This
trend would lead to blowups of the self couplings a little bit above the lower
scale that we focus on in this paper. Then it might not be so fruitful to
discuss issues originating from physics at a higher scale.

###  II.4 Vacuum stability against charged scalar trilinear couplings

Vacuum stability has to be especially assured in the Higgs potential against
contributions from electrically charged bosons
($h_{1}^{\pm},h_{2}^{\pm},k^{\pm\pm}$). However, our model has some loop
contributions to leading-order values of these quartic couplings via the
scalar trilinear couplings $\mu_{11}$ and $\mu_{22}$. When they are large, we
should examine the vacuum stability against the effect. Here, we examine this
issue at the one loop level. Let us describe these quartic couplings as
follows,

$\displaystyle 0\ {<}\
\lambda_{h_{1}}=\lambda^{(0)}_{h_{1}}+\lambda^{(1)}_{h_{1}},$ (II.12)
$\displaystyle 0\ {<}\
\lambda_{h_{2}}=\lambda^{(0)}_{h_{2}}+\lambda^{(1)}_{h_{2}},$ (II.13)
$\displaystyle 0\ {<}\ \lambda_{k}=\lambda^{(0)}_{k}+\lambda^{(1)}_{k},$
(II.14)

where the upper indices denote the numbers of the order in loops, and the one
loop contributions are given as

$\displaystyle\lambda^{(1)}_{h_{1}}=-8|\mu_{11}|^{4}F_{0}(m_{h^{\pm}_{1}},m_{k^{\pm\pm}}),$
(II.15)
$\displaystyle\lambda^{(1)}_{h_{2}}=-8|\mu_{22}|^{4}F_{0}(m_{h^{\pm}_{2}},m_{k^{\pm\pm}}),$
(II.16)
$\displaystyle\lambda^{(1)}_{k}=-4|\mu_{11}|^{4}F_{0}(m_{h^{\pm}_{1}},m_{h^{\pm}_{1}})-4|\mu_{22}|^{4}F_{0}(m_{h^{\pm}_{2}},m_{h^{\pm}_{2}}),$
(II.17)

with

$\displaystyle F_{0}(m_{a},m_{b})$
$\displaystyle=\frac{1}{(4\pi)^{2}}\int_{0}^{1}dxdy\delta(x+y-1)\frac{xy}{(xm_{a}^{2}+ym_{b}^{2})^{2}}$
$\displaystyle=\begin{cases}\displaystyle\frac{1}{(4\pi)^{2}}\frac{m_{a}^{2}\left(\log\left(\frac{m_{a}^{2}}{m_{b}^{2}}\right)-2\right)+m_{b}^{2}\left(\log\left(\frac{m_{a}^{2}}{m_{b}^{2}}\right)+2\right)}{\left(m_{a}^{2}-m_{b}^{2}\right){}^{3}}&(\text{for
}m_{a}\not=m_{b}),\\\
\displaystyle\frac{1}{(4\pi)^{2}}\frac{1}{6m_{a}^{4}}&(\text{for
}m_{a}=m_{b}),\end{cases}$ (II.18)

where the form of $\mu_{11}$ is shown in Eq. (II.3) and each of $m_{1}$ and
$m_{2}$ in $F_{0}$ represents a mass of propagating fields in the loops. We
include these constraints in the numerical analysis later. To avoid the global
minimum accompanying charge breaking, the following condition should at least
be satisfied:

$\displaystyle|\mu_{22}|<\sqrt{\Lambda}\left[m^{2}_{\Phi}+m^{2}_{h_{1}}+m^{2}_{h_{2}}+m^{2}_{k}+m^{2}_{\Sigma}\right]^{1/2},\quad\Lambda\equiv\sum_{i={\rm
all\ quartic\ couplings\ including\ \lambda_{11}}}\lambda_{i},$ (II.19)

where we assume the simplified configuration,
$r\equiv|\Phi|=|h^{+}_{1}|=|h^{+}_{2}|=|k^{++}|=|\Sigma_{0}|$ and the above
inequality comes from the requirement that $r$ does not have a finite nonzero
value. The summation is taken over the coefficients of the $17$ quartic terms
in Eq. (II.2) including $\Sigma_{0}^{\ast}h_{1}^{-}h_{1}^{-}k^{++}$ and its
Hermitian conjugate. When all of these quartic couplings are assumed to take
the same value $\lambda$, the above condition is rewritten as

$\displaystyle|\mu_{22}|\lesssim
4\sqrt{\lambda}\left[m^{2}_{h_{1}}+m^{2}_{h_{2}}+m^{2}_{k}+m^{2}_{\Sigma}\right]^{1/2},$
(II.20)

where the contributions via $m^{2}_{\Phi}$ and $\lambda_{\Phi}$ are
insignificant and thus neglected.

###  II.5 Neutrino mass matrix

Figure 1: A schematic description for the radiative generation of neutrino
masses.

A Majorana neutrino mass matrix $m_{\nu}$ is generated at the three loop level
via the diagram shown in Fig. 1, which corresponds to the coefficients of the
effective operators,
$-\frac{1}{2}(m_{\nu})_{ab}\times\overline{(\nu_{L_{a}})^{c}}\nu_{L_{b}}$. The
form of $(m_{\nu})_{ab}$ is evaluated by a straightforward calculation as

$\displaystyle(m_{\nu})_{ab}$
$\displaystyle=\frac{\mu_{11}\mu_{22}}{(4\pi)^{6}}\sum_{i,j,k=1}^{3}{\frac{1}{M_{k}^{4}}}\left[(y_{L})_{ai}m_{\ell_{i}}(y^{T}_{R})_{ik}(M_{N_{k}})(y_{R})_{kj}m_{\ell_{j}}(y^{T}_{L})_{jb}\right]$
$\displaystyle\quad\quad\times
F_{1}\left(\frac{m_{h^{+}_{1}}^{2}}{{M_{k}^{2}}},\frac{m_{h^{+}_{2}}^{2}}{{M_{k}^{2}}},\frac{m_{\ell_{i}}^{2}}{{M_{k}^{2}}},\frac{m_{\ell_{j}}^{2}}{{M_{k}^{2}}},\frac{{M_{N_{k}}}^{2}}{{M_{k}^{2}}},\frac{m_{k^{\pm\pm}}^{2}}{{M_{k}^{2}}}\right),$
(II.21)

where the mass scale ${M_{k}}={\rm
max}[m_{h^{\pm}_{1}},m_{h^{\pm}_{2}},m_{\ell_{i}},m_{\ell_{j}},M_{N_{k}},m_{k^{\pm\pm}}]$
is used for factorizing the loop function $F_{1}$ as a dimensionless variable.
Here, we take the relationship in masses in our setup,
$m_{\ell_{i}},m_{\ell_{j}}<M_{N_{k}}$, into consideration and then $M_{k}$ has
only the index $k$. $F_{1}$ is symbolically calculated as follows:

$\displaystyle
F_{1}\left(X_{1},X_{2},X_{3},X_{4},X_{5},X_{6}\right)=\int\mathbf{dX}\frac{1}{\Delta_{1}}\frac{1}{(\Delta_{2})^{2}}\frac{\rho}{(\Delta_{3})^{2}},$
(II.22) $\displaystyle\int\mathbf{dX}$
$\displaystyle=\int_{0}^{1}\\!\\!\\!dxdydz\,\delta(x+y+z-1)\int_{0}^{1}\\!\\!\\!d\alpha
d\beta d\gamma
d\delta\,\delta(\alpha+\beta+\gamma+\delta-1)\int_{0}^{1}\\!\\!\\!d\rho
d\sigma d\omega\,\delta(\rho+\sigma+\omega-1),$ (II.23)

with

$\displaystyle\Delta_{1}$ $\displaystyle=y(y-1)+z(z-1)+2yz,$ (II.24)
$\displaystyle\Delta_{2}$ $\displaystyle=(\alpha Y+\delta)^{2}-\delta-\alpha
Y^{2}-\alpha X,$ (II.25) $\displaystyle\Delta_{3}$ $\displaystyle=\rho
A\left(X_{1},X_{2},X_{3},X_{5},X_{6}\right)-\sigma X_{4}-\omega X_{1},$
(II.26) $\displaystyle A\left(X_{1},X_{2},X_{3},X_{5},X_{6}\right)$
$\displaystyle=-\frac{\alpha((x+y)X_{2}+zX_{5})}{((\alpha
Y+\delta)^{2}-\delta-\alpha Y^{2}-\alpha X)(y(y-1)+z(z-1)+2yz)}$
$\displaystyle\quad+\frac{\beta X_{1}+\gamma X_{3}+\delta X_{6}}{((\alpha
Y+\delta)^{2}-\delta-\alpha Y^{2}-\alpha X)},$ (II.27) $\displaystyle X$
$\displaystyle=-\left(\frac{y}{y+z}\right)^{2}+\frac{y(y-1)}{y(y-1)+z(z-1)+2yz},\quad
Y=\frac{y}{y+z}.$ (II.28)

Here, note that the shape of $F_{1}$ is completely the same as that in Ref.
Hatanaka:2014tba  except for the content of $X_{5}$ in Eq. (II.21).

Neutrino mass eigenstates and their mixings are evaluated by reflecting on
similarities to the Zee-Babu model Herrero-Garcia:2014hfa . The structure of
the fermion line is similar to the Zee-Babu model zee-babu , that is, a rank
two model of the neutrino mass matrix with a massless eigenstate due to the
antisymmetricity of $y_{L}$. Let us describe the neutrino mass matrix as

$\displaystyle(m_{\nu})_{ab}=(U_{{\text{PMNS}}}\,m^{\text{diag}}_{\nu}\,U^{T}_{{\text{PMNS}}})_{ab}\equiv{\sum_{i,j=1}^{3}}(y_{L})_{ai}\omega_{ij}(y^{T}_{L})_{jb},$
(II.29)
$\displaystyle\omega_{ij}={\sum_{k=1}^{3}}\,m_{\ell_{i}}(y^{T}_{R})_{ik}(\zeta_{k}M_{N_{k}})(y_{R})_{kj}m_{\ell_{j}}{,}$
(II.30)
$\displaystyle\zeta_{k}=\frac{\mu_{11}\mu_{22}}{(4\pi)^{6}{M^{4}_{k}}}F_{1}\left(\frac{m_{h^{+}_{1}}^{2}}{{M^{2}_{k}}},\frac{m_{h^{+}_{2}}^{2}}{{M^{2}_{k}}},0,0,\frac{{M_{N_{k}}}^{2}}{{M^{2}_{k}}},\frac{m_{k^{\pm\pm}}^{2}}{{M^{2}_{k}}}\right),$
(II.31)

where $i,j,k$ run over $1$ to $3$,
$m^{\text{diag}}_{\nu}\equiv(m_{1},m_{2},m_{3})$ are the neutrino mass
eigenvalues, and $U_{{\text{PMNS}}}$ (Pontecorvo-Maki-Nakagawa-Sakata matrix
Maki:1962mu ; Pontecorvo:1967fh ) is the mixing matrix to diagonalize the
neutrino mass matrix, which is parametrized as Herrero-Garcia:2014hfa

$\displaystyle
U_{{\text{PMNS}}}=\left[\begin{array}[]{ccc}{c_{13}}c_{12}&c_{13}s_{12}&s_{13}e^{-i\delta}\\\
-c_{23}s_{12}-s_{23}s_{13}c_{12}e^{i\delta}&c_{23}c_{12}-s_{23}s_{13}s_{12}e^{i\delta}&s_{23}c_{13}\\\
s_{23}s_{12}-c_{23}s_{13}c_{12}e^{i\delta}&-s_{23}c_{12}-c_{23}s_{13}s_{12}e^{i\delta}&c_{23}c_{13}\\\
\end{array}\right]\left[\begin{array}[]{ccc}1&0&0\\\ 0&e^{i\phi/2}&0\\\
0&0&1\\\ \end{array}\right],$ (II.38)

with $c_{ij}\equiv\cos\theta_{ij}$ and $s_{ij}\equiv\sin\theta_{ij}$. $\delta$
and $\phi$ represent the Dirac $CP$ phase and the Majorana one, respectively.
Here, we treat the two ratios in $F_{1}$, ${m_{\ell_{i}}}/{M_{k}}$ and
${m_{\ell_{j}}}/{M_{k}}$ as zero values approximately, which means that the
loop factor has no dependence on $i$ or $j$ as $\zeta_{k}$.

Requirements for the observed neutral profiles depend on the ordering of the
neutrino masses, which are normal [$m_{1}\,(=0)<m_{2}<m_{3}$] or inverted
[$m_{3}\,(=0)<m_{1}<m_{2}$].444More details can be given in Ref. Herrero-
Garcia:2014hfa  for both of the cases. When we consider the normal ordering,
the following relations should hold for realizing the observed neutrino
profiles Herrero-Garcia:2014hfa ,

$\displaystyle(y_{L})_{13}$
$\displaystyle=(s_{12}c_{23}/(c_{12}c_{13})+s_{13}s_{23}e^{-i\delta}/c_{13})(y_{L})_{23},$
$\displaystyle(y_{L})_{12}$
$\displaystyle=(s_{12}s_{23}/(c_{12}c_{13})-s_{13}c_{23}e^{-i\delta}/c_{13})(y_{L})_{23},$
$\displaystyle{\left(\frac{(y_{L})_{23}}{2}\right)^{2}}\omega_{33}$
$\displaystyle\approx
m_{3}c^{2}_{13}s^{2}_{23}+m_{2}e^{i\phi}(c_{12}c_{23}-e^{i\delta}s_{12}s_{13}s_{23})^{2},$
$\displaystyle{\left(\frac{(y_{L})_{23}}{2}\right)^{2}}\omega_{23}$
$\displaystyle\approx-
m_{3}c^{2}_{13}c_{23}s_{23}+m_{2}e^{i\phi}(c_{12}s_{23}+e^{i\delta}c_{23}s_{12}s_{13})(c_{12}c_{23}-e^{i\delta}s_{12}s_{13}s_{23}),$
$\displaystyle{\left(\frac{(y_{L})_{23}}{2}\right)^{2}}\omega_{22}$
$\displaystyle\approx
m_{3}c^{2}_{13}c^{2}_{23}+m_{2}e^{i\phi}(c_{12}s_{23}{+}e^{i\delta}c_{23}s_{12}s_{13})^{2},$
(II.39)

where we use $m_{e}\,{\ll}\,m_{\mu},m_{\tau}$. In the case of the inverted
neutrino mass hierarchy, the conditions are deformed Herrero-Garcia:2014hfa ,

$\displaystyle(y_{L})_{13}$
$\displaystyle=-(c_{13}s_{23}e^{-i\delta}/s_{13})(y_{L})_{23},$
$\displaystyle(y_{L})_{12}$
$\displaystyle=+(c_{13}c_{23}e^{-i\delta}/s_{13})(y_{L})_{23},$
$\displaystyle{\left(\frac{(y_{L})_{23}}{2}\right)^{2}}\omega_{33}$
$\displaystyle\approx
m_{1}(c_{23}s_{12}+e^{i\delta}c_{12}s_{13}s_{23})^{2}+m_{2}e^{i\phi}(c_{12}c_{23}-e^{i\delta}s_{12}s_{13}s_{23})^{2},$
$\displaystyle{\left(\frac{(y_{L})_{23}}{2}\right)^{2}}\omega_{23}$
$\displaystyle\approx
m_{1}(s_{12}s_{23}-e^{i\delta}c_{12}c_{23}s_{13})(c_{23}s_{12}+e^{i\delta}c_{12}s_{13}s_{23})$
$\displaystyle\qquad\qquad\quad+m_{2}e^{i\phi}(c_{12}s_{23}+e^{i\delta}c_{23}s_{12}s_{13})(c_{12}c_{23}-e^{i\delta}s_{12}s_{13}s_{23}),$
$\displaystyle{\left(\frac{(y_{L})_{23}}{2}\right)^{2}}\omega_{22}$
$\displaystyle\approx
m_{1}(s_{12}s_{23}-e^{{+i}\delta}c_{12}c_{23}s_{13})^{2}+m_{2}e^{i\phi}(c_{12}s_{23}{+}e^{i\delta}c_{23}s_{12}s_{13})^{2}.$
(II.40)

In the numerical analysis in the following Sec. III.2, we adopt the latest
values of $\theta_{12}$, $\theta_{23}$, $\theta_{13}$, $\Delta m_{21}^{2}$ and
$\Delta m_{3\ell}^{2}$ ($\ell=1$ for the normal hierarchy and $\ell=2$ for the
inverted one) found in Ref. Gonzalez-Garcia:2014bfa . The two $CP$-violating
phases $\delta$ and $\phi$ are treated as free parameters.555The best fit
values are as follows: $\sin^{2}{\theta_{12}}=0.304$,
$\sin^{2}{\theta_{23}}=0.452$, $\sin^{2}{\theta_{13}}=0.0218$, $\Delta
m_{21}^{2}=7.50\times 10^{-5}\,\text{eV}^{2}$, $\Delta m_{31}^{2}=+2.457\times
10^{-3}\,\text{eV}^{2}$ (in the normal hierarchy);
$\sin^{2}{\theta_{12}}=0.304$, $\sin^{2}{\theta_{23}}=0.579$,
$\sin^{2}{\theta_{13}}=0.0219$, $\Delta m_{21}^{2}=7.50\times
10^{-5}\,\text{eV}^{2}$, $\Delta m_{32}^{2}=-2.449\times
10^{-3}\,\text{eV}^{2}$ (in the inverted hierarchy) Gonzalez-Garcia:2014bfa .

###  II.6 Lepton flavor violations and the universality of charged currents

In the present model, owing to interactions containing charged scalars, new
contributions to several lepton-flavor-violating processes are found at the
tree or the one loop level. Also, universalities of charged currents are
threatened by the vertices. Some of them are very similar to the case of the
Zee-Babu model discussed, e.g., in Ref. Nebot:2007bc , while situations are
changed in the others. Here, we briefly summarize important points:

  * •

$\ell_{i}\to\ell_{j}\nu\bar{\nu}$: like in the Zee-Babu model, the processes
receive additional contributions that are never found in the SM via the
exchange of the singly charged scalar coupling to the charged leptons
$h_{1}^{\pm}$, which have an influence on the value of the Fermi constant via
charged-lepton decay. If the masses of the right-hand neutrinos are small,
there is a possibility that the channels $\ell_{i}\to\ell_{j}N_{1}\bar{N_{1}}$
via the exchange of $h_{2}^{\pm}$ are open and they fake contributions to
$\ell_{i}\to\ell_{j}\nu\bar{\nu}$. In this paper, we do not consider such a
light $N_{1}$. Then, we simply ignore this effect in the following analysis
part.

  * •

$\ell_{i}^{-}\to\ell_{f}^{-}\gamma$: like in the Zee-Babu model, these decay
sequences (when $m_{\ell_{i}^{-}}>m_{\ell_{f}^{-}}$) are radiatively generated
at the one loop level via the diagrams containing charged scalars which the
photon couples to. In these processes, both of the interactions,
$\bar{L}^{c}_{L_{i}}L_{L_{j}}h^{+}_{1}$ and
$\bar{N}_{R_{i}}e^{c}_{R_{j}}h_{2}^{-}$ (also with their Hermitian
conjugations included) make nonzero contributions to the decay widths. The
doubly charged scalar $k^{\pm\pm}$ has no contribution at the one loop level
due to the absence of tree-level interactions with the charged leptons.

  * •

$\ell_{i}^{-}\to\ell_{j}^{+}\ell_{k}^{-}\ell_{l}^{-}$: being different from
the Zee-Babu model, absence of direct interactions between the charged leptons
and the doubly charged scalar $k^{\pm\pm}$ generates the situation that these
processes are induced at the one loop level. Two types of box diagrams
contributing to the processes are there: containing (i)
$\bar{L}^{c}_{L_{i}}L_{L_{j}}h^{+}_{1}$ interactions; and containing (ii)
$\bar{N}_{R_{i}}e^{c}_{R_{j}}h_{2}^{-}$ ones (where their Hermitian
conjugations are also included), while two types of photon and $Z$ penguin
diagrams are also generated. Here, the contribution from the penguins is
directory related to the dipole diagrams of
$\ell_{i}^{-}\to\ell_{f}^{-}\gamma$ in the language of branching ratios, e.g.,
when we consider $\mu\to 3e$ Arganda:2005ji ; Toma:2013zsa ; Chowdhury:2015sla
as

$\displaystyle\text{Br}(\mu\to
3e)\simeq\frac{\alpha_{\text{EM}}}{3\pi}\left(\log\left(\frac{m_{\mu}^{2}}{m_{e}^{2}}\right)-\frac{11}{4}\right)\text{Br}(\mu\to
e\gamma),$ (II.41)

where $\alpha_{\text{EM}}$ is the electromagnetic fine structure constant and
the penguin-type contributions to $\mu\to 3e$ are suppressed compared with
those to $\mu\to e\gamma$. Since the experimental bounds on
$\ell_{i}^{-}\to\ell_{f}^{-}\gamma$ and
$\ell_{i}^{-}\to\ell_{j}^{+}\ell_{k}^{-}\ell_{l}^{-}$ are roughly comparable,
we simply ignore the contributions from the penguin-type diagrams in the
calculation of the constraints in this category.

Now we show the explicit forms of constraints on the present model. We recast
the results of the analysis on the Zee-Babu model in Ref. Herrero-
Garcia:2014hfa  for our case, where the three $\ell\to\ell\gamma$ processes
($\mu^{-}\to e^{-}\gamma$, $\tau^{-}\to e^{-}\gamma$,
$\tau^{-}\to\mu^{-}\gamma$), the four types of gauge universalities
(lepton/hadron, $\mu/e$, $\tau/\mu$, $\tau/e$), and the seven $\ell\to 3\ell$
processes ($\mu^{-}\to e^{+}e^{-}e^{-}$, $\tau^{-}\to e^{+}e^{-}e^{-}$,
$\tau^{-}\to e^{+}e^{-}\mu^{-}$, $\tau^{-}\to e^{+}\mu^{-}\mu^{-}$,
$\tau^{-}\to\mu^{+}e^{-}e^{-}$, $\tau^{-}\to\mu^{+}e^{-}\mu^{-}$,
$\tau^{-}\to\mu^{+}\mu^{-}\mu^{-}$) are included.

Process | $(i,f)$ | Experimental bounds ($90\%$ C.L.) | $C_{if}$
---|---|---|---
$\mu^{-}\to e^{-}\gamma$ | $(2,1)$ | $\text{Br}(\mu\to e\gamma)<5.7\times 10^{-13}$ | $1.6\times 10^{-6}$
$\tau^{-}\to e^{-}\gamma$ | $(3,1)$ | $\text{Br}(\tau\to e\gamma)<3.3\times 10^{-8}$ | $0.52$
$\tau^{-}\to\mu^{-}\gamma$ | $(3,2)$ | $\text{Br}(\tau\to\mu\gamma)<4.4\times 10^{-8}$ | $0.7$
Table 2: Summary of the coefficient $C_{if}$ in $\ell\to\ell\gamma$ processes and experimental data used in the analysis in Herrero-Garcia:2014hfa . Type of universality | Experimental bounds ($90\%$ C.L.)
---|---
Lepton/hadron universality | $\sum_{q=d,s,b}|V_{uq}^{\text{exp}}|=0.9999\pm 0.0006$
$\mu/e$ universality | $G_{\mu}^{\text{exp}}/G_{e}^{\text{exp}}=1.0010\pm 0.0009$
$\tau/\mu$ universality | $G_{\tau}^{\text{exp}}/G_{\mu}^{\text{exp}}=0.9998\pm 0.0013$
$\tau/e$ universality | $G_{\tau}^{\text{exp}}/G_{e}^{\text{exp}}=1.0034\pm 0.0015$
Table 3: Summary of the experimental data on universality of charged currents used in the analysis in Herrero-Garcia:2014hfa . Process | $(i,j,k,l)$ | Experimental bounds ($90\%$ C.L.) | $A_{ijkl}$
---|---|---|---
$\mu^{-}\to e^{+}e^{-}e^{-}$ | $(2,1,1,1)$ | $\text{Br}<1.0\times 10^{-12}$ | $2.3\times 10^{-5}$
$\tau^{-}\to e^{+}e^{-}e^{-}$ | $(3,1,1,1)$ | $\text{Br}<2.7\times 10^{-8}$ | $0.009$
$\tau^{-}\to e^{+}e^{-}\mu^{-}$ | $(3,1,1,2)$ | $\text{Br}<1.8\times 10^{-8}$ | $0.005$
$\tau^{-}\to e^{+}\mu^{-}\mu^{-}$ | $(3,1,2,2)$ | $\text{Br}<1.7\times 10^{-8}$ | $0.007$
$\tau^{-}\to\mu^{+}e^{-}e^{-}$ | $(3,2,1,1)$ | $\text{Br}<1.5\times 10^{-8}$ | $0.007$
$\tau^{-}\to\mu^{+}e^{-}\mu^{-}$ | $(3,2,1,2)$ | $\text{Br}<2.7\times 10^{-8}$ | $0.007$
$\tau^{-}\to\mu^{+}\mu^{-}\mu^{-}$ | $(3,2,2,2)$ | $\text{Br}<2.1\times 10^{-8}$ | $0.008$
Table 4: Summary of the coefficient $A_{ijkl}$ in $\ell\to 3\ell$ processes
and experimental data used in the analysis in Herrero-Garcia:2014hfa .

  * •

$\ell_{i}^{-}\to\ell_{f}^{-}\gamma$ processes: in this case, the result of
recasting is

$\displaystyle\frac{\left|\sum_{a=1}^{3}\left[(y_{L}^{\dagger})_{af}(y_{L})_{ia}\right]^{2}\left(I_{1,a}I_{2,a}+I_{1,a}^{2}\right)+\sum_{a=1}^{3}\left[(y_{R})_{af}(y_{R}^{\dagger})_{ia}\right]^{2}\left(I^{\prime}_{1,a}I^{\prime}_{2,a}+{I^{\prime}_{1,a}}^{2}\right)\right|}{16m_{h_{1}^{\pm}}^{4}\left|\sum_{a=1}^{3}\left(I_{1,a}I_{2,a}+I_{1,a}^{2}\right)\right|}<\frac{C_{if}}{\text{TeV}^{4}},$
(II.42)

with the loop functions

$\displaystyle I_{1,a}$
$\displaystyle={\frac{1}{(4\pi)^{2}}}{\int_{0}^{1}\\!\\!dx\int_{0}^{1-x}\\!\\!\\!\\!\\!\\!\\!dy}\
\frac{x(2x-1)}{(x+y)m_{h_{1}^{\pm}}^{2}+(1-x-y)m_{\nu_{a}}^{2}}\simeq-{\frac{1}{(4\pi)^{2}}}\frac{1}{36m_{h_{1}^{\pm}}^{2}},$
(II.43) $\displaystyle I_{2,a}$
$\displaystyle={\frac{1}{(4\pi)^{2}}}{\int_{0}^{1}\\!\\!dx\int_{0}^{1-x}\\!\\!\\!\\!\\!\\!\\!dy}\
\frac{x(2y-1)}{(x+y)m_{h_{1}^{\pm}}^{2}+(1-x-y)m_{\nu_{a}}^{2}}\simeq-{\frac{1}{(4\pi)^{2}}}\frac{5}{36m_{h_{1}^{\pm}}^{2}},$
(II.44) $\displaystyle I^{\prime}_{1,a}$
$\displaystyle={\frac{1}{(4\pi)^{2}}}{\int_{0}^{1}\\!\\!dx\int_{0}^{1-x}\\!\\!\\!\\!\\!\\!\\!dy}\
\frac{x(2x-1)}{(x+y)m_{h_{2}^{\pm}}^{2}+(1-x-y)M_{N_{a}}^{2}},$ (II.45)
$\displaystyle I^{\prime}_{2,a}$
$\displaystyle={\frac{1}{(4\pi)^{2}}}{\int_{0}^{1}\\!\\!dx\int_{0}^{1-x}\\!\\!\\!\\!\\!\\!\\!dy}\
\frac{x(2y-1)}{(x+y)m_{h_{2}^{\pm}}^{2}+(1-x-y)M_{N_{a}}^{2}},$ (II.46)

where we use $m_{\nu_{a}}\simeq 0$. Here, we treat the final-state lepton
$\ell^{-}_{f}$ as a massless particle. Concrete forms of the integrals are
summarized in Appendix A. The dimensionless coefficient $C_{if}$ representing
the digits in Herrero-Garcia:2014hfa  (before recasting) are summarized in
Table 2. The factor $16$ comes from the difference in the coupling convention
of the interaction $\bar{L}^{c}_{L_{i}}L_{L_{j}}h^{+}_{1}$
($\bar{L}^{c}_{L_{i}}L_{L_{j}}h^{+}$ in the Zee-Babu model). These decay
processes are one loop induced ones in both of the models, and thus the loop
factor $1/(4\pi)^{2}$ in the integrals is canceled out in the final form in
Eq. (II.42).

  * •

Gauge coupling universalities: in this category, recasting is just
straightforward by the replacement $h^{\pm}\to h^{\pm}_{1}$,

$\displaystyle\left|{\frac{({y}_{L})_{12}}{2}}\right|^{2}$
$\displaystyle<0.007\ \left(\frac{{m_{h_{1}^{\pm}}}}{{\rm
TeV}}\right)^{2}\quad(\text{lepton/hadron universality}),$
$\displaystyle\left|\left|{\frac{({y}_{L})_{23}}{2}}\right|^{2}-\left|{\frac{({y}_{L})_{13}}{2}}\right|^{2}\right|$
$\displaystyle<0.024\ \left(\frac{{m_{h_{1}^{\pm}}}}{{\rm
TeV}}\right)^{2}\quad(\text{$\mu/e$ universality}),$
$\displaystyle\left|\left|{\frac{({y}_{L})_{13}}{2}}\right|^{2}-\left|{\frac{({y}_{L})_{12}}{2}}\right|^{2}\right|$
$\displaystyle<0.035\ \left(\frac{{m_{h_{1}^{\pm}}}}{{\rm
TeV}}\right)^{2}\quad(\text{$\tau/\mu$ universality}),$
$\displaystyle\left|\left|{\frac{({y}_{L})_{23}}{2}}\right|^{2}-\left|{\frac{({y}_{L})_{12}}{2}}\right|^{2}\right|$
$\displaystyle<0.04\ \left(\frac{{m_{h_{1}^{\pm}}}}{{\rm
TeV}}\right)^{2}\quad\ \,(\text{$\tau/e$ universality}).$ (II.47)

The corresponding experimental bounds are summarized in Table 3.

  * •

$\ell\to 3\ell$ processes: all of the cases are summarized symbolically as

$\displaystyle\frac{1}{4}{\Big{|}}\left|A_{\nu}+B_{\nu}\right|^{2}+\left|A_{N}+B_{N}\right|^{2}-2\,\text{Re}\left[A_{N}C_{N}^{\ast}\right]-2\,\text{Re}\left[B_{N}C_{N}^{\ast}\right]+\frac{1}{2}\left|C_{N}\right|^{2}{\Big{|}^{1/2}}<\frac{A_{ijkl}}{\text{TeV}^{2}},$
(II.48)

with the effective couplings

$\displaystyle A_{\nu}$
$\displaystyle=(y_{L}y_{L}^{\dagger})_{ik}(y_{L}y_{L}^{\dagger})_{jl}\,J_{1,0},$
(II.49) $\displaystyle B_{\nu}$
$\displaystyle=(y_{L}y_{L}^{\dagger})_{il}(y_{L}y_{L}^{\dagger})_{jk}\,J_{1,0},$
(II.50) $\displaystyle A_{N}$
$\displaystyle=\sum_{a,b=1}^{3}(y_{R})_{al}(y_{R})_{bk}(y_{R}^{\dagger})_{ja}(y_{R}^{\dagger})_{ib}\,J_{1,ab},$
(II.51) $\displaystyle B_{N}$
$\displaystyle=\sum_{a,b=1}^{3}(y_{R})_{al}(y_{R})_{bk}(y_{R}^{\dagger})_{jb}(y_{R}^{\dagger})_{ia}\,J_{1,ab},$
(II.52) $\displaystyle C_{N}$
$\displaystyle=2\sum_{a,b=1}^{3}(y_{R})_{al}(y_{R})_{ak}(y_{R}^{\dagger})_{jb}(y_{R}^{\dagger})_{ib}M_{N_{a}}M_{N_{b}}J_{2,ab}.$
(II.53)

The loop functions $J_{1,ab}$ and $J_{2,ab}$ are given as

$\displaystyle J_{1,ab}$
$\displaystyle=\frac{1}{(4\pi)^{2}}{\int_{0}^{1}\\!\\!dx\int_{0}^{1-x}\\!\\!\\!\\!\\!\\!\\!dy}\
\frac{1-x-y}{xM_{N_{a}}^{2}+yM_{N_{b}}^{2}+(1-x-y)m_{h_{2}^{\pm}}^{2}},$
(II.54) $\displaystyle J_{2,ab}$
$\displaystyle=\frac{1}{(4\pi)^{2}}{\int_{0}^{1}\\!\\!dx\int_{0}^{1-x}\\!\\!\\!\\!\\!\\!\\!dy}\
\frac{1-x-y}{\left[xM_{N_{a}}^{2}+yM_{N_{b}}^{2}+(1-x-y)m_{h_{2}^{\pm}}^{2}\right]^{2}},$
(II.55)

where the form of $J_{1,0}$ is obtained by the replacements $M_{N_{a}}\to
m_{\nu_{a}}(\simeq 0)$, $M_{N_{b}}\to m_{\nu_{b}}(\simeq 0)$ and
$m_{h_{2}^{\pm}}\to m_{h_{1}^{\pm}}$ as

$\displaystyle
J_{1,0}\simeq\frac{1}{2}\frac{1}{(4\pi)^{2}}\frac{1}{m_{h_{1}^{\pm}}^{2}}.$
(II.56)

Also in this calculation, we treat the final-state leptons $\ell^{+}_{j}$,
$\ell^{-}_{k}$ and $\ell^{-}_{l}$ as massless particles. Other concrete forms
of the integrals are summarized in Appendix A. Here, the dimensionless
valuables $A_{ijkl}$ indicate the numerical values (in the analysis in
Herrero-Garcia:2014hfa  before recasting) and their concrete values are
summarized in Table 4. The factor $1/4$ originates from recasting. The loop
suppression factor $1/(4\pi)^{2}$ should appear in Eq. (II.48) since these
processes are generated at the tree level in the original Zee-Babu model.

Finally, we briefly have a discussion on the LFV via $k^{\pm\pm}$ exchange in
our model. When we focus on the LFV process accompanying two neutrinos,
$\ell^{-}\to k^{--}\ell^{+}\to
2(h_{1}^{-})^{\ast}+\ell^{+}\to(2\ell^{-}+2\nu)+\ell^{+}$, where the first
step is one loop induced and the two $h_{1}^{-}$ as intermediate states are
usually off shell since the mass of $h_{1}^{-}$ is bounded from below, where
$m_{h_{1}^{\pm}}$ should be more than a few TeV; see Sec. III for details. The
order of the loops is the same as that of $\ell\to 3\ell$. But, the off
shellness suppresses the decay width of this LFV process. Therefore we can
simply ignore this possibility.

##  III Valid parameter region with neutrino sector

###  III.1 Form of $y_{R}$

Different from the Zee-Babu model, $y_{R}$ is not a symmetric matrix in
general but $y_{R}^{T}\,(\zeta M_{N})\,y_{R}\,(\equiv Y_{R})$ is symmetric
one. Hence, we can utilize fitting relations in the Zee-Babu model with
modifications. Here, we adopt the following forms in $y_{R}$ [and in $y_{N}$
as shown in Eq. (II.4)] for simplicity,

$\displaystyle y_{R}=\left[\begin{array}[]{ccc}\ast&\ast&\ast\\\ \ast&a&b\\\
\ast&b&c\\\ \end{array}\right],\quad M_{N}={\rm
diag}\,(M_{N_{1}},M_{N_{2}},M_{N_{3}}),$ (III.4)

where $a,b,c$ are arbitrary complex numbers. The correspondence to the factors
$\omega_{ij}$ in Eq. (II.30) is as follows:

$\displaystyle\omega_{ij}=m_{\ell_{i}}(Y_{R})_{ij}m_{\ell_{j}}.$ (III.5)

Here, because of the mass hierarchy of the charged leptons $m_{\ell_{1}}\ll
m_{\ell_{2}}<m_{\ell_{3}}$, the terms in $\omega_{ij}$ including
$m_{\ell_{1}}$ are negligible. Combining this issue with the hierarchy assumed
in the right-hand neutrinos $M_{N_{1}}<M_{N_{2}}<M_{N_{3}}$, we conclude that
the elements of $y_{R}$ expressed by $\ast$ in Eq. (III.4) do not affect the
values of $\omega_{ij}$ significantly. Therefore in the following analysis, we
only consider the couplings $a,b,c$ of $y_{R}$ in fitting the neutrino
profiles.

Now, we rewrite the relations in Eq. (III.5) as follows:

$\displaystyle\frac{\omega_{22}}{m_{\ell_{2}}^{2}}$
$\displaystyle=(Y_{R})_{22}\simeq\left[a^{2}(\zeta_{2}M_{N_{2}})+b^{2}(\zeta_{3}M_{N_{3}})\right],$
$\displaystyle\frac{\omega_{23}}{m_{\ell_{2}}m_{\ell_{3}}}$
$\displaystyle=(Y_{R})_{23}\simeq\left[ab\,(\zeta_{2}M_{N_{2}})+bc\,(\zeta_{3}M_{N_{3}})\right],$
$\displaystyle\frac{\omega_{33}}{m_{\ell_{3}}^{2}}$
$\displaystyle=(Y_{R})_{33}\simeq\left[b^{2}(\zeta_{2}M_{N_{2}})+c^{2}(\zeta_{3}M_{N_{3}})\right].$
(III.6)

As shown in Eqs. (II.39)–(II.40), when we fix the values of $(y_{L})_{23}$,
$\delta$ and $\phi$, the values of $\omega_{22}$, $\omega_{23}$, $\omega_{33}$
[and also $(y_{L})_{12}$, $(y_{L})_{13}$] are automatically determined through
the relations. In each scanning in the following section, we pick up a
solution on $a,b,c$ of the above simultaneous equations.

###  III.2 Parameter scanning

Adopting the structure of $y_{R}$ in the previous subsection, we execute
parameter scans to search for consistent regions in the parameter space. In
this model, $k^{\pm\pm}$ does not contribute to the processes with LFV
significantly. Thus, we consider the two possibilities $m_{k^{\pm\pm}}=250$
and $500\,\text{GeV}$, while the other two singly charged scalars
$h_{1^{\pm}},h_{2^{\pm}}$ have a few TeV masses. Here, we select the four
choices: $m_{h_{1}^{\pm}}=m_{h_{2}^{\pm}}=3,\,3.5,\,4,$ and $4.8\,\text{TeV}$.
For brevity, we fix the three right-hand neutrinos as follows:
$M_{N_{1}}=m_{h}/2$, $M_{N_{2}}=5\,\text{TeV}$ and $M_{N_{3}}=10\,\text{TeV}$,
where we mention that $N_{1}$ is a “Higgs-portal” dark matter candidate. The
mass is assumed to be around the $125\,\text{GeV}$ Higgs resonant region.
Detailed discussion on this topic is found in Sec. V.

We mention that compliance with the relations in the case of the normal
hierarchy in Eq. (II.39) or the inverted hierarchy in Eq. (II.40) leads to the
situation that only the parameters $\mu_{11}$, $\mu_{22}$, $\delta$, $\phi$,
$(y_{L})_{23}$ are free to be chosen (in addition to the masses of particles
that we fixed in the above paragraph). Note that the five matrix elements of
$y_{R}$, namely $(y_{R})_{11}$, $(y_{R})_{12}$, $(y_{R})_{13}$,
$(y_{R})_{21}$, and $(y_{R})_{31}$, are ineffective in the determination of
the active neutrino profiles, while their nonzeroness possibly modifies the
strengths of the lepton-flavor-violating processes significantly. Then we
consider nonzero values in the five couplings for the sake of completeness. In
each scanning, we randomly select values of the ten parameters within the
corresponding ranges

$\displaystyle\mu_{12}=\mu_{22}\
{(\equiv\mu)}\in[1\,\text{TeV},20\,\text{TeV}],\quad\delta\in[0,2\pi],\quad\phi\in[0,2\pi],\quad(y_{L})_{23}\in[-1,1],$
$\displaystyle(y_{R})_{ij}\in[-0.01,-0.1]\,\cup\,[0.01,0.1]\quad\Big{(}\text{for
}(i,j)=(1,1),(1,2),(1,3),(2,1),(3,1)\Big{)}.$ (III.7)

Our definition of an allowed point is a set of parameters where all the
following requirements are satisfied:

  * •

Observed values in masses and mixings of the three active neutrinos are
generated. As we mentioned above, realization of this requirement is equal to
the compliance with the relations in Eq. (II.39) or (II.40).

  * •

Predictions for LFVs do not exceed the bounds shown in Tables 2
($\ell\to\ell\gamma$), 3 (gauge universalities), and 4 ($\ell\to 3\ell$).
Concretely, evading bounds corresponds to meeting the inequalities in Eqs.
(II.42) and (II.47)–(II.48), respectively.

  * •

Fulfilling the requirements on vacuum stability in Eqs. (II.12)–(II.14) and
(II.20). In the last condition, we adopt the following modified form with
physical masses,

$\displaystyle|\mu_{22}|\lesssim
4\sqrt{\lambda}\left[m^{2}_{h_{1}^{\pm}}+m^{2}_{h_{2}^{\pm}}+m^{2}_{k^{\pm\pm}}\right]^{1/2},$
(III.8)

where we simply change the parameters $m_{h_{1}},\,m_{h_{2}},\,m_{k}$ to their
physical masses and ignore the mass parameter $m_{\Sigma}$. The value of the
right-hand side would be near the original one and it would be useful for
roughly estimating this type of bound. In the following analysis, we adopt the
initial conditions for the quartic couplings,

$\displaystyle\lambda^{(0)}_{h_{1}}=\lambda^{(0)}_{h_{2}}=\lambda^{(0)}_{k}=4\pi,\quad\lambda=4\pi.$
(III.9)

  * •

Ensuring perturbativity, all the couplings should be equal to or less than
$4\pi$.

$m_{k^{\pm\pm}}$ | $(m_{h_{1}^{\pm}},m_{h_{2}^{\pm}})$ | $F_{1}^{(\zeta_{2})}$ | $F_{1}^{(\zeta_{3})}$ | Number of allowed points
---|---|---|---|---
| $(4.8\,\text{TeV},4.8\,\text{TeV})$ | $-1.16818$ | $-11.5428$ | $2804$
$500$ | $(4.0\,\text{TeV},4.0\,\text{TeV})$ | $-2.14198$ | $-20.0287$ | $422$
GeV | $(3.5\,\text{TeV},3.5\,\text{TeV})$ | $-3.30195$ | $-29.6525$ | $89$
| $(3.0\,\text{TeV},3.0\,\text{TeV})$ | $-5.37518$ | $-46.0932$ | $16$
| $(4.8\,\text{TeV},4.8\,\text{TeV})$ | $-1.21466$ | $-11.9695$ | $1644$
$250$ | $(4.0\,\text{TeV},4.0\,\text{TeV})$ | $-2.24746$ | $-20.9444$ | $190$
GeV | $(3.5\,\text{TeV},3.5\,\text{TeV})$ | $-3.49141$ | $-31.2321$ | $31$
| $(3.0\,\text{TeV},3.0\,\text{TeV})$ | $-5.74153$ | $-49.0141$ | $2$
Table 5: The numbers of consistent points in the normal hierarchy obtained by
parameter scanning. $F_{1}^{(\zeta_{2})}$ and $F_{1}^{(\zeta_{3})}$ are the
loop functions in $\zeta_{2}$ and $\zeta_{3}$ shown in Eq. (II.31),
respectively. Kindly refer to the main body of this subsection for details of
this scanning.

Figure 2: Histograms showing distributions of $\delta$, $\phi$,
$(y_{L})_{23}$, $|(y_{R})_{22}|$ and $\mu$ of the result of a scanning with
the charged scalar masses
$m_{k^{\pm\pm}}=500\,\text{GeV},\,m_{h_{1}^{\pm}}=m_{h_{2}^{\pm}}=4.8\,\text{TeV}$.
The total number of the consistent data points is $2804$.

Figure 3:  Correlations between $\phi$ and $\delta$, $(y_{L})_{23}$ and
$|(y_{R})_{22}|$, and $\mu$ and $|(y_{R})_{22}|$ are shown. The data set is
the same as in Fig. 2.

Our result in the normal hierarchy is summarized in Table 5. In each of the
eight combinations of the charged scalar masses, we randomly scan $10^{5}$
points in the ten parameters, where the range we consider is shown in Eq.
(III.7). Apart from the previous work Hatanaka:2014tba , the result says that
the mass of the doubly charged scalar can be light around a few hundred GeV,
whereas the other two singly charged ones should be heavy (at least) around a
few TeV. In this scenario, it is very hard to produce $h_{1}^{\pm}$ and
$h_{2}^{\pm}$ in colliders. On the other hand, detecting $k^{\pm\pm}$ could be
a reasonable option for probing this model in present and future collider
experiments.666Lower bounds at $95\%$ CL on $m_{k^{\pm\pm}}$ via
$8\,\text{TeV}$ data were provided by the ATLAS group in Ref. ATLAS:2014kca
as $374\,\text{GeV}$, $402\,\text{GeV}$, $438\,\text{GeV}$ when assuming a
$100\%$ branching ratio to $e^{\pm}e^{\pm}$, $e^{\pm}\mu^{\pm}$,
$\mu^{\pm}\mu^{\pm}$ pairs, respectively. In our model, the decay sequence
$k^{\pm\pm}\to\ell^{\pm}\ell^{\pm}$ is one-loop induced and many parameters
contribute. In this paper, we skip to have detailed discussions on prospects
in colliders.

Next, we look into properties of the parameters in the allowed region. Here,
we only show the case of
$m_{k^{\pm\pm}}=500\,\text{GeV},\,m_{h_{1}^{\pm}}=m_{h_{2}^{\pm}}=4.8\,\text{TeV}$
as an example, where we checked that other cases are not so different. In Fig.
2, we display the distributions of $\delta$, $\phi$, $(y_{L})_{23}$,
$|(y_{R})_{22}|$ and $\mu$ based on $2804$ data points of the allowed region
as histograms.777The value of $(y_{R})_{22}$ is complex in general obtained as
a part of a solution of the simultaneous equations in Eq. (III.6). The
correlations between $\phi$ and $\delta$, $(y_{L})_{23}$ and $|(y_{R})_{22}|$,
and $\mu$ and $|(y_{R})_{22}|$ are also shown in Fig. 3. In the present three
loop scenario, the loop factor $1/(4\pi)^{6}$ in the neutrino masses tends to
suppress the realized masses very much. Thus, at least a part of the
parameters related with the masses would be sufficiently large. In the
following, we investigate details.

  * •

We can see that this model could give some preference to $\delta$ around
$\pi$. On the other hand, $\delta$ not around $\pi$ is also realized. Few
trends are seen in the distribution of $\phi$.

  * •

A typical digit of the absolute value of $(y_{L})_{23}$ is in $0.5\sim 0.6$,
which would be greater compared with other radiative neutrino models in one or
two loop level. Since $|\omega_{22}|\simeq|\omega_{23}|\simeq|\omega_{33}|$,
the other two nonzero components of $y_{L}$, $(y_{L})_{12}$ and
$(y_{L})_{13}$, tend to have almost the same order of magnitude.

  * •

Like the Zee-Babu model, via the relation
$|\omega_{22}|\simeq|\omega_{23}|\simeq|\omega_{33}|$, the three components of
the effective symmetric Yukawa couplings $Y_{R}$ meet the definite hierarchy,

$\displaystyle{|(Y_{R})_{33}|:|(Y_{R})_{23}|:|(Y_{R})_{22}|}\sim\frac{m_{\mu}^{2}}{m_{\tau}^{2}}:\frac{m_{\mu}}{m_{\tau}}:1,$
(III.10)

which means that the original $|(y_{R})_{22}|$ tends to hold a significantly
large value. The peak of the distribution is around $9$, which is still rather
small compared with the perturbative upper bound $4\pi$. Because of this
characteristic, the masses of the singly charged scalars $m_{h_{1}^{\pm}}$ and
$m_{h_{2}^{\pm}}$ should be greater than around $3\,\text{TeV}$ to circumvent
the bounds.

  * •

The common coefficients of the trilinear terms $\mu$ should be large at around
$14\sim 15\,\text{TeV}$ to compensate for the suppression factor in the
neutrino masses. From the perturbativity in the couplings $\lambda_{11}$ in
the form of the trilinear coupling $\mu_{11}$ in Eq. (II.3) and $y_{N}$ found
in the masses of the right-hand neutrinos in Eq. (II.4), $v^{\prime}$ cannot
be so small. A reasonable value of $v^{\prime}$ is
$\mathcal{O}(1)\,\text{TeV}$.

Finally, we briefly comment on the possibility with the inverted mass
hierarchy. In our model, like the Zee-Babu model, the components of the active
neutrino mass matrix $(m_{\nu})_{ab}$ contain the charged lepton masses, and
also the Majorana masses with the assumed mass hierarchy
$M_{N_{1}}<M_{N_{2}}<M_{N_{3}}$. Then the normal hierarchy would be
preferable. Within our search among $10^{5}$ points in the inverted case, even
in the choice of
$(m_{k^{\pm\pm}},m_{h_{1}^{\pm}},m_{h_{2}^{\pm}})=(500\,\text{GeV},4.8\,\text{TeV},4.8\,\text{TeV})$,
no solution is found.

##  IV Constraint from LHC Higgs search

Process | ATLAS | CMS | Reference
---|---|---|---
$h\to\gamma\gamma$ | $1.17^{+0.27}_{-0.27}$ | $1.12\pm 0.24$ | Aad:2014eha ; Khachatryan:2014jba
$h\to ZZ^{\ast}\to 4\ell$ | $1.44^{+0.40}_{-0.33}$ | $1.00\pm 0.29$ | Aad:2014eva ; Khachatryan:2014jba
$h\to WW^{\ast}\to\ell\nu\ell\nu$ | $1.09^{+0.23}_{-0.21}$ | $0.83\pm 0.21$ | ATLAS:2014aga ; Khachatryan:2014jba
$h\to b\bar{b}$ | $0.5^{+0.4}_{-0.4}$ | $0.84\pm 0.44$ | Aad:2014xzb ; Khachatryan:2014jba
$h\to\tau\bar{\tau}$ | $1.4^{+0.4}_{-0.4}$ | $0.91\pm 0.28$ | ATLAS_Higgs_tautau ; Khachatryan:2014jba
Table 6: Summary of the latest LHC Higgs search data as (combined) signal
strengths in the five Higgs decay patterns ($h\to\gamma\gamma,h\to
ZZ^{\ast}\to 4\ell,h\to WW^{\ast}\to\ell\nu\ell\nu,h\to
b\bar{b},h\to\tau\bar{\tau}$). The analyses are based on the data sets
accumulated in the LHC first run, whose details of the ATLAS and the CMS are
$4.5$ –
$4.7\,\text{fb}^{-1}\,(7\,\text{TeV})+20.3\,\text{fb}^{-1}\,(8\,\text{TeV})$,
$5.1\,\text{fb}^{-1}\,(7\,\text{TeV})+19.7\,\text{fb}^{-1}\,(8\,\text{TeV})$,
respectively.

In this part, we evaluate constraints on the parameter space of scalars in
this model by use of the latest results of LHC Higgs searches by the ATLAS and
CMS experiment groups. First, we describe the method we use for global
analysis. Like in the papers Giardino:2012dp ; Carmi:2012in ; Ellis:2012hz ,
we adopt the following form of signal strength of the single Higgs production
channel with the subsequent Higgs decay to the particles $f$,

$\displaystyle\mu_{f}=\frac{\sigma_{\text{total}}}{\sigma_{\text{total}}^{\text{SM}}}\times\frac{\text{Br}_{h\to
f}}{\text{Br}_{h\to f}^{\text{SM}}},$ (IV.1)

where $\sigma_{\text{total}}$ ($\sigma_{\text{total}}^{\text{SM}}$) represents
the total production cross section of a Higgs boson in this model (SM),
respectively. The Higgs branching ratios to the particles $f$ of the SM and
this model are discriminated by having or not having the superscript SM.

Here, note that in our scenario, the observed Higgs boson is a mixture of the
$SU(2)_{L}$ doublet $\Phi$ and the singlet $\Sigma_{0}$ as shown in Eq.
(II.10) and no additional colored particle is introduced, which means the
absence of new contributions to the Higgs production via the gluon fusion
process. Since the SM gauge bosons and quarks do not couple to $\Sigma_{0}$,
the ratio of the total cross sections is easily calculated as

$\displaystyle\frac{\sigma_{\text{total}}}{\sigma_{\text{total}}^{\text{SM}}}=\cos^{2}{\alpha}.$
(IV.2)

Whereas, evaluating the ratio of the branding ratios is rather complicated.
First, we look at the following decomposition,

$\displaystyle\frac{\text{Br}_{h\to f}}{\text{Br}_{h\to
f}^{\text{SM}}}=\left(\frac{\Gamma_{h}}{\Gamma^{\text{SM}}_{h}}\right)^{-1}\times\frac{\Gamma_{h\to
f}}{\Gamma^{\text{SM}}_{h\to f}},$ (IV.3)

where $\Gamma_{h}^{(\text{SM})}$ represents the corresponding Higgs total
width. In the Higgs decay, the presence of the charged scalars in our model
($h_{1}^{\pm},h_{2}^{\pm},k^{\pm\pm}$) also modifies the partial decay widths
of the photon-associated decay processes, $h\to\gamma\gamma$ and $h\to
Z\gamma$, in addition to the doublet-singlet mixing effect. As we discussed in
Sec. III.2, the two singly charged scalars tend to be very massive at around a
few TeV at least. In such a situation, we can completely neglect the effect
via the two particles and then only the contribution from $k^{\pm\pm}$ is
included as a loop effect in the following part.

The ratio of the total widths is rewritten as follows:

$\displaystyle\frac{\Gamma_{h}}{\Gamma^{\text{SM}}_{h}}=\text{Br}^{\text{SM}}_{h\to\text{SM\,others}}\times\frac{\Gamma_{h\to\text{SM\,others}}}{\Gamma_{h\to\text{SM\,others}}^{\text{SM}}}+\text{Br}^{\text{SM}}_{h\to\gamma\gamma}\times\frac{\Gamma_{h\to\gamma\gamma}}{\Gamma^{\text{SM}}_{h\to\gamma\gamma}}+\text{Br}^{\text{SM}}_{h\to
Z\gamma}\times\frac{\Gamma_{h\to Z\gamma}}{\Gamma^{\text{SM}}_{h\to
Z\gamma}}+\frac{\Gamma_{h\to\text{inv}}}{\Gamma^{\text{SM}}_{h}},$ (IV.4)

where $\Gamma_{h\to\text{inv}}$ expresses the Higgs partial decay width to
invisible pairs, which is written as

$\displaystyle\Gamma_{h\to\text{inv}}=\Gamma_{h\to GG}+\Gamma_{h\to
N_{R_{1}}\bar{N}_{R_{1}}}.$ (IV.5)

Concrete forms of partial decay width are summarized in Appendix B.

$\text{Br}^{\text{SM}}_{h\to\text{SM\,others}}$ can be described as
$1-\text{Br}^{\text{SM}}_{h\to\gamma\gamma}-\text{Br}^{\text{SM}}_{h\to
Z\gamma}$ very precisely. In the following analysis, we use the values
$\text{Br}^{\text{SM}}_{h\to\gamma\gamma}=2.28\times
10^{-3},\,\text{Br}^{\text{SM}}_{h\to Z\gamma}=1.54\times 10^{-3}$ and
$\Gamma^{\text{SM}}_{h}=4.07\,\text{MeV}$ in $m_{h}=125\,\text{GeV}$ reported
by the LHC Higgs Cross Section Working Group Heinemeyer:2013tqa . The ratios
of the partial decay widths are described with the help of the formulas for
$h\to\gamma\gamma$ and $h\to Z\gamma$ in Refs. Ellis:1975ap ; Shifman:1979eb ;
Djouadi:2005gi ; Carena:2012xa ; Chen:2013vi  as

$\displaystyle\frac{\Gamma_{h\to f}}{\Gamma^{\text{SM}}_{h\to f}}$
$\displaystyle=\cos^{2}{\alpha}\quad(f\in(\text{others in
SM}))\quad\to\quad\frac{\Gamma_{h\to\text{SM\,others}}}{\Gamma_{h\to\text{SM\,others}}^{\text{SM}}}=\cos^{2}{\alpha},$
(IV.6)
$\displaystyle\frac{\Gamma_{h\to\gamma\gamma}}{\Gamma^{\text{SM}}_{h\to\gamma\gamma}}$
$\displaystyle=\left|\cos{\alpha}+\frac{1}{2}\frac{v^{2}}{m_{k^{\pm\pm}}^{2}}Q_{k}^{2}C_{hkk}\frac{A_{0}^{\gamma\gamma}(\tau_{k})}{A_{1}^{\gamma\gamma}(\tau_{W})+N_{C}Q_{t}^{2}A_{1/2}^{\gamma\gamma}(\tau_{t})}\right|^{2},$
(IV.7) $\displaystyle\frac{\Gamma_{h\to Z\gamma}}{\Gamma^{\text{SM}}_{h\to
Z\gamma}}$
$\displaystyle=\left|\cos{\alpha}-\frac{v^{2}}{m_{k^{\pm\pm}}^{2}}\left(2Q_{k}g_{Zkk}\right)C_{hkk}\frac{A_{0}^{Z\gamma}(\tau_{k},\lambda_{k})}{vA_{\text{SM}}^{Z\gamma}}\right|^{2},$
(IV.8)

with effective couplings

$\displaystyle C_{hkk}=\cos{\alpha}\lambda_{\Phi
k}-\sin{\alpha}\left(\frac{v^{\prime}}{v}\right)\lambda_{\Sigma k},\quad
g_{Zkk}=-Q_{k}\left(\frac{s_{W}}{c_{W}}\right)$ (IV.9)

and the form factor

$\displaystyle
A_{\text{SM}}^{Z\gamma}=\frac{2}{v}\left[\cot{\theta_{W}}A_{1}^{Z\gamma}(\tau_{W},\lambda_{W})+N_{C}\frac{(2Q_{t})(T_{3}^{(t)}-2Q_{t}s_{W}^{2})}{s_{W}c_{W}}A_{{1/2}}^{Z\gamma}(\tau_{t},\lambda_{t})\right].$
(IV.10)

$N_{C}\,(=3)$, $Q_{t}\,(=2/3)$, $Q_{k}\,(=2)$, $T_{3}^{(t)}\,(=1/2)$, $c_{W}$
and $s_{W}$ are the QCD color factor for quarks, the electric charges of the
top quark, the doubly charged scalar in unit of the positron’s one, the weak
isospin of the top quark, and the cosine and the sine of the Weinberg angle
$\theta_{W}$, respectively. The loop factors take the following forms,

$\displaystyle A_{1}^{\gamma\gamma}(x)$
$\displaystyle=-x^{2}\left[2x^{-2}+3x^{-1}+3(2x^{-1}-1)f(x^{-1})\right],$
(IV.11) $\displaystyle A_{1/2}^{\gamma\gamma}(x)$
$\displaystyle=2x^{2}\left[x^{-1}+(x^{-1}-1)f(x^{-1})\right],$ (IV.12)
$\displaystyle A_{0}^{\gamma\gamma}(x)$
$\displaystyle=-x^{2}\left[x^{-1}-f(x^{-1})\right],$ (IV.13) $\displaystyle
A_{1}^{Z\gamma}(x,y)$
$\displaystyle=4(3-\tan^{2}{\theta_{W}})I_{2}(x,y)+\left[(1+2x^{-1})\tan^{2}{\theta_{W}}-(5+2x^{-1})\right]I_{1}(x,y),$
(IV.14) $\displaystyle A_{1/2}^{Z\gamma}(x,y)$
$\displaystyle=I_{1}(x,y)-I_{2}(x,y),$ (IV.15) $\displaystyle
A_{0}^{Z\gamma}(x,y)$ $\displaystyle=I_{1}(x,y),$ (IV.16)

with the functions

$\displaystyle I_{1}(x,y)$
$\displaystyle=\frac{xy}{2(x-y)}+\frac{x^{2}y^{2}}{2(x-y)^{2}}\left[f(x^{-1})-f(y^{-1})\right]+\frac{x^{2}y}{(x-y)^{2}}\left[g(x^{-1})-g(y^{-1})\right],$
(IV.17) $\displaystyle I_{2}(x,y)$
$\displaystyle=-\frac{xy}{2(x-y)}\left[f(x^{-1})-f(y^{-1})\right].$ (IV.18)

In the above formulas, forms of the input variables to the loop factors
$\tau_{i}$ and $\lambda_{i}$ are defined as fractions by the Higgs boson mass
($m_{h}$) or the Z boson mass ($m_{Z}$)

$\displaystyle\tau_{i}=\frac{4m_{i}^{2}}{m_{h}^{2}},\quad\lambda_{i}=\frac{4m_{i}^{2}}{m_{Z}^{2}}\quad(i=t,W,k).$
(IV.19)

The two ratios usually take values above one ($m_{h}\leq 2m_{i},\,m_{Z}\leq
2m_{i}$). The two functions $f(z)$ and $g(z)$ ($z\equiv x^{-1}\ \text{or}\
y^{-1}$) are formulated as

$\displaystyle f(z)$
$\displaystyle=\begin{cases}\arcsin^{2}{\sqrt{z}}&\text{for }z\leq 1,\\\
-\frac{1}{4}\left[\log\left(\frac{1+\sqrt{1-z^{-1}}}{1-\sqrt{1-z^{-1}}}\right)-i\pi\right]^{2}&\text{for
}z>1,\end{cases}$ (IV.20) $\displaystyle g(z)$
$\displaystyle=\begin{cases}\sqrt{z^{-1}-1}\arcsin{\sqrt{z}}&\text{for }z\leq
1,\\\
\frac{\sqrt{1-z^{-1}}}{2}\left[\log\left(\frac{1+\sqrt{1-z^{-1}}}{1-\sqrt{1-z^{-1}}}\right)-i\pi\right]&\text{for
}z>1,\end{cases}$ (IV.21)

where the situation $m_{h}\leq 2m_{i},\,m_{Z}\leq 2m_{i}$ corresponds to
$z\leq 1$.

Figure 4: The $1\sigma$ (blue) and $2\sigma$ (red) allowed regions of the
global analysis based on the LHC data summarized in Table 6. From top left to
bottom, we choose the parameters $(\lambda_{\Phi k},\lambda_{\Sigma
k}v^{\prime})=(1.0,v),\,(\pi,v),\,(0.1,v)$, respectively.

Figure 5: The $1\sigma$ (blue) and $2\sigma$ (red) allowed regions of the
global analysis based on the LHC data summarized in Table 6. From left to
right, we choose the parameters $(\lambda_{\Phi k},\lambda_{\Sigma
k}v^{\prime})=(1.0,5v),\,(1.0,10v)$, respectively.

To estimating the consistent parameter region with the latest results of the
Higgs search at the LHC, we define the $\chi^{2}$ valuable by use of the
signal strength in Eq. (IV.1) as follows:

$\displaystyle\chi^{2}=\sum_{f=\gamma\gamma,ZZ^{\ast},WW^{\ast},\atop
b\bar{b},\tau\bar{\tau}\,(\text{ATLAS})}\left(\frac{\mu_{f}-\hat{\mu}_{f}}{\hat{\sigma}_{f}}\right)^{2}+\sum_{f=\gamma\gamma,ZZ^{\ast},WW^{\ast},\atop
b\bar{b},\tau\bar{\tau}\,(\text{CMS})}\left(\frac{\mu_{f}-\hat{\mu}_{f}}{\hat{\sigma}_{f}}\right)^{2}.$
(IV.22)

Here, we take the results of the five Higgs decay channels reported by the
ATLAS and the CMS experiments into consideration, which are
$h\to\gamma\gamma,h\to ZZ^{\ast}\to 4\ell,h\to WW^{\ast}\to\ell\nu\ell\nu,h\to
b\bar{b},h\to\tau\bar{\tau}$ summarized in Table 6. The two hatted symbols
$\hat{\mu}_{f}$ and $\hat{\sigma}_{f}$ represent the corresponding central
value and error, respectively. We assumed that each of ten experimental inputs
follows the Gaussian distribution and there are no correlations among them.
Also, when an error is asymmetric, we adopt its simple average as an input
value of the corresponding $\hat{\sigma}_{f}$ for analysis. These
simplifications are justified for our purpose of roughly estimating survived
regions of the parameter space of this model.

In the following analysis, for simplicity, we focus on a reasonable situation
in which the mass of the DM is around $m_{h}/2$ where the Higgs invisible
channel to a pair of the DMs is near the threshold and negligible. Detailed
discussions on the DM candidate are given in Sec. V. Now, apparently, we see
the five parameters $m_{k^{\pm\pm}}$, $\lambda_{\Phi k}$, $\lambda_{\Sigma
k}$, $v^{\prime}$ and $\alpha$ govern the signal strengths of the Higgs boson.
But in fact, as shown in Eq. (IV.9), $\lambda_{\Sigma k}$ and $v^{\prime}$
contribute to the Higgs physics only as the combination of $\lambda_{\Sigma
k}v^{\prime}$. Thus, four independent degrees of freedom are relevant in
total.

We fix the two parameters $\lambda_{\Phi k}$ and $\lambda_{\Sigma
k}v^{\prime}$ in each of the following global fits for simplicity. Here, we
consider the five possibilities, $(\lambda_{\Phi k},\lambda_{\Sigma
k}v^{\prime})=(1.0,v)$, $(\pi,v)$, $(0.1,v)$, $(1.0,5v)$ and $(1.0,10v)$ and
search for their global minima on the remaining two variables
$(m_{k^{\pm\pm}},\sin{\alpha})$ at first. The positions of each corresponding
minimum are $(m_{k^{\pm\pm}},\sin{\alpha})=({74.7}\,\text{GeV},-0.0759)$,
$(109\,\text{GeV},-0.0774)$, $(65.0\,\text{GeV},0.102)$,
$(80.3\,\text{GeV},-0.0768)$ and $(87.5\,\text{GeV},-0.0804)$ with
$\chi^{2}_{\text{min}}\simeq 5.08$ (commonly in all the cases), respectively.
The $1\sigma$ and $2\sigma$ boundaries are defined as the positions with
$\chi^{2}=\chi^{2}_{\text{min}}+\Delta\chi^{2}_{1\sigma,2\sigma}$, where the
values of $\Delta\chi^{2}_{1\sigma,2\sigma}$ are calculated by the cumulative
distribution function of the $\chi^{2}$-distribution with two degrees of
freedom as $\Delta\chi^{2}_{1\sigma}=2.296$, $\Delta\chi^{2}_{2\sigma}=6.180$.
The results are shown in Figs. 4–5. Small disconnected regions are found
possibly due to accidental cancellations when $(\lambda_{\Phi
k},\lambda_{\Sigma k}v^{\prime})=(\pi,v)$ in Fig. 4 and $(\lambda_{\Phi
k},\lambda_{\Sigma k}v^{\prime})=(1.0,10v)$ in Fig. 5. When we consider the
doubly charged scalars with $m_{k^{\pm\pm}}=500\,\text{GeV}$, the range of
$\alpha$

$\displaystyle|\sin{\alpha}|\ {\lesssim}\ {0.3}$ (IV.23)

is roughly allowed within $2\sigma$ confidence levels, where actual upper
bounds depend on contents of the two prefixed parameters.

##  V Issues on $N_{R_{1}}$ dark matter

In this section, we have discussions on dark matter related issues of this
model. After the breakdown of the new global $U(1)$ symmetry, an accidental
$Z_{2}$ symmetry still remains, which ensures the existence of a dark matter
candidate. Under the assignment of the $U(1)$ charge in Table 1,
$N_{R_{i}}\,(i=1,2,3)$ and $h^{\pm}_{2}$ hold negative parities, where the
lightest one among them gets to be absolutely stable. After considering the
aptness discussed in Sec. III.2 that $h^{\pm}_{2}$ (and also $h^{\pm}_{1}$)
should be sufficiently heavy at around a few TeV to evade the constraints,
e.g., via the processes with LFV, the lightest right-hand neutrino, namely
$N_{R_{1}}$, in our setup is usually stabilized by the symmetry and plays a
significant role as dark matter.

An important point of our dark sector is in the mechanism for generating
Majorana mass terms of $N_{R}$ that originates from the spontaneous global
$U(1)$ breaking accompanying a nonzero VEV of the field $\Sigma_{0}$. Through
the mixing between $\Sigma_{0}$ and the Higgs doublet $\Phi$, $N_{R_{1}}$ can
communicate with the SM particles, where the two neutral scalars $h$ and $H$
work as mediators. Here, we expect that a suitable amount of dark matter
relics are left in the Universe through resonant effects via the $h$ or $H$
pole, where the mass of $N_{R_{1}}$ is fixed around $m_{h}/2$ or $m_{H}/2$,
respectively.

The relic abundance and the spin-independent cross section for direct
detection of the candidate $N_{R_{1}}$ are calculable by following a standard
method. In the following part, we use the shorthand notation $\chi$ for
showing the DM $N_{R_{1}}$ (in the mass eigenstate as a Majorana fermion).
Within the freeze-out approximation, the present-day relic density is
evaluated as Griest:1990kh

$\displaystyle\Omega_{\chi}h^{2}\simeq\frac{1.07\times
10^{9}\,\text{GeV}^{-1}}{\sqrt{g_{\ast}(x_{f})}M_{\text{Pl}}J(x_{f})},$ (V.1)

where $\Omega_{\chi}$, $h\,(\simeq 0.7)$, $M_{\text{Pl}}$ and $g_{\ast}$
express the present-day energy density of $\chi$, the present-day scaled
Hubble parameter, the Planck mass and the number of the relativistic degrees
of freedom, respectively. Note that the latest value of $\Omega_{\chi}h^{2}$
reported by the Planck experiment is $0.1196\pm 0.0031\,(68\%\,\text{{C.L.}})$
Ade:2013zuv . The valuable $x_{f}$ is defined by the dark matter mass
$m_{\chi}$ and the temperature at the freeze-out $T_{f}$ as $m_{\chi}/T_{f}$.
The efficiency of the annihilation after the freeze-out is described through
the integral $J$:

$\displaystyle J(x_{f})=\int_{x_{f}}^{\infty}dx\frac{\langle\sigma
v_{\text{rel}}\rangle}{x^{2}},$ (V.2)

where $\langle\sigma v_{\text{rel}}\rangle$ stands for the thermally averaged
annihilation cross section multiplied with the relativistic relative velocity
$v_{\text{rel}}$. In this work, we adopt the sophisticated way for taking
thermal average relativistically discussed in Refs. Gondolo:1990dk ;
Edsjo:1997bg , where $\langle\sigma v_{\text{rel}}\rangle$ is estimated as

$\displaystyle\langle\sigma
v_{\text{rel}}\rangle=\frac{\displaystyle\int_{4m_{\chi}^{2}}^{\infty}ds\frac{\sqrt{s-4m_{\chi}^{2}}}{16}W_{\chi\chi}K_{1}\left(\frac{\sqrt{s}}{T}\right)}{\displaystyle
m_{\chi}^{4}T\left[K_{2}\left(\frac{m_{\chi}}{T}\right)\right]^{2}}=\frac{\displaystyle\int_{4m_{\chi}^{2}}^{\infty}ds{\sqrt{s-4m_{\chi}^{2}}}W_{\chi\chi}K_{1}\left(\frac{\sqrt{s}}{m_{\chi}}x\right)}{\displaystyle
16m_{\chi}^{5}x^{-1}\left[K_{2}\left(x\right)\right]^{2}},$ (V.3)

where $x$ is defined as $m_{\chi}/T$ (like $x_{f}$) at the temperature $T$,
$K_{1,2}$ are the modified Bessel functions of the second kind of order $1$
and $2$, respectively, and $W_{\chi\chi}$ is a Lorentz invariant variable
describing the cross section multiplied with the Lorentz invariant flux factor
$4E_{\chi}^{2}v_{\text{rel}}$. $W_{\chi\chi}$ is formulated in the center of
mass system with the integration over the solid angle as

$\displaystyle
W_{\chi\chi}=\sum_{f}\frac{1}{32\pi^{2}}\sqrt{1-\frac{4m_{f}^{2}}{s}}\int
d\Omega\left|\overline{\mathcal{M}}(\chi{\bar{\chi}}\to f\bar{f})\right|^{2}$
(V.4)

where we sum over all possible two-body final states consisting of the same
particle $f$. The detail of the amplitudes is found in Appendix C. In the
following numerical calculation, we simply use the fixed reasonable values
$x_{f}=20$ and $g_{\ast}(x_{f})=100$ throughout the calculations for brevity.

Our evaluation of the $\chi$-nucleon spin-independent cross section is based
on the discussions in Refs. Hisano:2010ct ; Abe:2015rja . We consider the
following effective Lagrangian for calculating the cross section at the
leading order,

$\displaystyle\mathcal{L}_{\text{eff}}=\sum_{{q}=u,d,s}f_{q}m_{q}\bar{\chi}\chi\bar{q}q-\frac{\alpha_{s}}{4\pi}f_{G}\bar{\chi}\chi
G^{a}_{\mu\nu}G^{a\mu\nu},$ (V.5)

where $q$, $m_{q}$, $\alpha_{s}$, and $G^{a}_{\mu\nu}$ represent the
corresponding quark fields, the quark masses, the QCD coupling strength, and
the field strength of the gluon, respectively. The coefficients $f_{q}$
($f_{G}$) determine the effective interactions between the quarks (gluon) and
the DM $\chi$. The corresponding values in our model are

$\displaystyle
f_{q}=f_{G}=\frac{1}{2}\left(\frac{m_{\chi}}{vv^{\prime}}\right)\left(-\frac{1}{m_{h}^{2}}+\frac{1}{m_{H}^{2}}\right)\cos{\alpha}\sin{\alpha}.$
(V.6)

The spin-independent cross section with the target nucleus $T$ is formulated
by use of the spin-independent coupling of $\chi$ with nucleon
$f_{N}\,(N=p,n)$ as

$\displaystyle\sigma_{\text{SI}}^{T}=\frac{4}{\pi}\mu_{T\chi}^{2}\left|n_{p}f_{p}+n_{n}f_{n}\right|^{2},$
(V.7)

where $\mu_{T\chi}$ is the reduced mass among the nucleus and the dark matter
is defined as $\mu_{T\chi}\equiv m_{T}m_{\chi}/(m_{T}+m_{\chi})$. $m_{T}$
shows the mass of the target nucleus $T$, which contains $n_{p}$ and $n_{n}$
numbers of protons ($p$) and neutrons ($n$), respectively. The effective
$\chi$-nucleon ($N$) coupling $f_{N}$ can be written down by the coefficients
of the effective operators in $\mathcal{L}_{\text{eff}}$ ($f_{q},f_{G}$) and
matrix elements ($f_{Nq},f_{NG}$) as

$\displaystyle f_{N}/m_{N}=\sum_{q=u,d,s}f_{q}f_{Nq}+\frac{2}{9}f_{G}f_{NG},$
(V.8)

with the concrete forms of the matrix elements,

$\displaystyle\bra{N}m_{q}\bar{q}q\ket{N}=f_{Nq}m_{N},\quad
1-\sum_{q=u,d,s}f_{Nq}=f_{NG}{,}$ (V.9)

where the value of $f_{NG}$ is calculable through the latter relation by the
use of $f_{Nq}$. Like in Ref. Abe:2015rja , we adopt the following default
values in the program micrOMEGAs Belanger:2013oya :
$f_{nu}=0.0110,\,f_{nd}=0.0273,\,f_{ns}=0.0447$.

The latest bound on the spin-independent scattering process was reported by
the LUX experiment as an upper limit on the spin-independent (elastic) DM-
nucleon cross section, which is approximately $10^{-45}\,\text{cm}^{2}$ (when
$m_{\chi}\sim 10^{2}\,\text{GeV}$) with the $90\%$ confidence level
Akerib:2013tjd . Here, we choose the neutron as the nucleon for putting a
bound on our parameter space via the LUX result. Now, the (spin-independent)
$\chi$-neutron cross section is calculated with ease as

$\displaystyle\sigma_{\text{SI}}^{N=n}$
$\displaystyle=\frac{4}{\pi}\mu_{n\chi}^{2}\left|m_{n}\left(\sum_{q=u,d,s}f_{nq}+\frac{2}{9}f_{nG}\right)f_{q}\right|^{2}$
$\displaystyle=\frac{1}{\pi}\mu_{n\chi}^{2}m_{n}^{2}\,C^{2}\left(\frac{m_{\chi}\cos{\alpha}\sin{\alpha}}{vv^{\prime}}\right)^{2}\left(-\frac{1}{m_{h}^{2}}+\frac{1}{m_{H}^{2}}\right)^{2},$
(V.10)

with

$\displaystyle C=\frac{2}{9}+\frac{7}{9}\sum_{q=u,d,s}f_{nq}\simeq 0.287,$
(V.11)

where similar calculations were done, e.g., in Refs. Baek:2011aa ; Baek:2012se
.

Figure 6: Realized present-day relic densities of the dark matter candidate
$\chi$ are shown as functions of $m_{\chi}$. The other effective parameters
are fixed as
$\sin{\alpha}=0.4,\,v^{\prime}=800\,\text{GeV},\,m_{H}=250\,\text{GeV},\,m_{k^{\pm\pm}}=500\,\text{GeV}$
(left panel) and
$\sin{\alpha}=0.5,\,v^{\prime}=5\,\text{TeV},\,m_{H}=250\,\text{GeV},\,m_{k^{\pm\pm}}=500\,\text{GeV}$
(right panel) for demonstration purposes, respectively. The area inside the
two horizontal dashed lines (where the splitting is almost invisible) suggests
the $2\sigma$ consistent region with the Planck data [$0.1196\pm
0.0031\,(68\%\,\text{C.L.})$ Ade:2013zuv ]. In the current setup, the relic
density sharply drops around the two resonant regions around $m_{h}/2\simeq
62.5\,\text{GeV}$ and $m_{H}/2\simeq 125\,\text{GeV}$, where in the latter
case, shown in the right panel, dropping is not sufficient for obtaining a
proper amount of dark matter even around $m_{H}/2$ at present. We neglect the
three- and four-body final states via virtual $W$ and $Z$ boson decays, which
gives sizable modifications near the thresholds for producing gauge boson
pairs Cline:2012hg ; Cline:2013gha , since our interest is only around
$m_{h}/2$ and $m_{H}/2$ where this correction is expected to be insignificant.

Figure 7: The matrix plots indicate suitable choices of the two parameters
$v^{\prime}$ and $\sin{\alpha}$ when we consider the situation of
$m_{H}=250\,\text{GeV},\,m_{k^{\pm\pm}}=250\,\text{GeV}$. We obtain the
observed relic density in the red (blue) region in the left (right) panel when
$m_{\chi}$ is around $m_{h}/2\,\simeq 62.5\,\text{GeV}$ ($m_{H}/2\,\simeq
125\,\text{GeV}$), respectively. The green region is excluded via excess of
the invisible decay channel of the observed Higgs boson. No excluded region is
found by the direct detection in the shown parameter range.

Figure 8: The matrix plots for showing the region where a proper amount of the
relic density is realized when
$m_{H}=250\,\text{GeV},\,m_{k^{\pm\pm}}=500\,\text{GeV}$. Conventions are the
same with ones in Fig. 7. No excluded region is found by the direct detection
in the shown parameter range.

Figure 9: Similar plots as Figs. 7 and 8 when
$m_{H}=500\,\text{GeV},\,m_{k^{\pm\pm}}=250\,\text{GeV}$ (left) and
$m_{H}=500\,\text{GeV},\,m_{k^{\pm\pm}}=500\,\text{GeV}$ (right),
respectively. Here, no solution is found around $m_{\chi}\simeq m_{H}/2\simeq
250\,\text{GeV}$, while one around $m_{\chi}\simeq m_{h}/2\simeq
62.5\,\text{GeV}$ is still there. The green region is excluded by the
invisible decay of the $125\,\text{GeV}$ Higgs. No excluded region is found by
the direct detection in the shown parameter range.

In the following calculation, we treat the variables $\alpha$, $m_{H}$, and
$v^{\prime}$ independently, which determines the coefficients
$\lambda_{\Phi\Sigma}$, $\lambda_{\Phi}$, and $\lambda_{\Sigma}$ through the
relations in Eqs. (II.8)–(II.9) as

$\displaystyle\lambda_{\Phi\Sigma}$
$\displaystyle=\frac{\sin{\alpha}\cos{\alpha}\,(m_{H}^{2}-m_{h}^{2})}{vv^{\prime}},$
$\displaystyle\lambda_{\Phi}$
$\displaystyle=\frac{\cos^{2}{\alpha}\,m_{h}^{2}+\sin^{2}{\alpha}\,m_{H}^{2}}{2v^{2}},$
$\displaystyle\lambda_{\Sigma}$
$\displaystyle=\frac{\sin^{2}{\alpha}\,m_{h}^{2}+\cos^{2}{\alpha}\,m_{H}^{2}}{2{v^{\prime}}^{2}}.$
(V.12)

We set two quartic couplings $\lambda_{\Phi k}$ and $\lambda_{\Sigma k}$ as
$1.0$ and $0.1$. In our scenario, $v^{\prime}$ tends to be
$\mathcal{O}(1)\,\text{TeV}$ leading to suppression of the thermally averaged
cross section describing pair annihilation processes; see Appendix C for
detail forms. It would require a mechanism for enhancing the cross section.
Hence, in this manuscript, we focus on the two resonant regions of $m_{\chi}$
around $m_{h}/2$ or $m_{H}/2$. Our requirement for the amount of the relics is
that it should be within the $\pm 2\sigma$ range of the latest value reported
by the Planck experiment [$0.1196\pm 0.0031\,(68\%\,\text{C.L.})$ Ade:2013zuv
]. In the following matrix plots, the values of $v^{\prime}$ and
$\sin{\alpha}$ at the center of a square represent the inputs in the whole
region shown by the square.

An important point is that all the effective diagrams of this pair
annihilation of $\chi$ are Higgs mediated.888Note that the processes
$N_{R_{1}}\bar{N}_{R_{1}}\to\ell^{+}_{i}\ell^{-}_{j}\,(i,j=1,2,3)$ with
exchanging $h^{\pm}_{2}$ in $t$-channel, which are important in similar
situations Krauss:2002px ; Cheung:2004xm , are ineffective and neglected in
the present scenario because the lower bound on $m_{h^{\pm}_{2}}$ is around a
few TeV and decoupled as shown in subsection III.2. Like the couplings between
the observed fermions and the $125\,\text{GeV}$ Higgs boson in the SM, the DM-
portal scalar couplings are proportional to the factor $m_{\chi}/v^{\prime}$.
In Fig. 6, we present realized present-day relic densities of the dark matter
candidate $\chi$ as functions of $m_{\chi}$, where the other effective
parameters are fixed as
$\sin{\alpha}=0.4,\,v^{\prime}=800\,\text{GeV},\,m_{H}=250\,\text{GeV},\,m_{k^{\pm\pm}}=500\,\text{GeV}$
(left panel) and
$\sin{\alpha}=0.5,\,v^{\prime}=5\,\text{TeV},\,m_{H}=250\,\text{GeV},\,m_{k^{\pm\pm}}=500\,\text{GeV}$
(right panel), respectively. Here, we can see that a suppression due to a huge
$v^{\prime}$ in the latter case and the point around the heavier resonance of
$H$ is not suitable for a correct amount of DM today. Around the lighter one
of $h$, less $s$-channel propagator suppression could be expected and thus the
solution around this pole still remains.

Another critical variable in the relic calculation is the mixing angle
$\alpha$ among the two $CP$ even scalars $h$ and $H$. If $\sin{\alpha}$ is
close to zero, contributions from the mediator $h$ are diminished
significantly due to almost zeroness of the couplings between $\chi$ and the
SM particles in the final state. Thereby, a suitable amount of the mixing is
requested for the solution near $m_{h}/2$. In the following analysis, we focus
on non-negative values in $\sin\alpha$ since the sign of $\sin\alpha$ is
insignificant in the resonant regions.999We mention details of this point. The
sign of $\sin\alpha$ affects results only through the amplitude of
$\chi\bar{\chi}\to hh$ and the partial decay width of $h/H\to\gamma\gamma$ and
$h/H\to Z\gamma$ in the $s$-channel propagators. In the DM masses which we
consider, the amplitude gives only a subleading contribution. The effect
coming from the widths is also negligible since the branching ratios of
$h/H\to\gamma\gamma$ and $h/H\to Z\gamma$ are pretty small.  Also, we ignore
the range of $|\sin\alpha|>0.5$ since no allowed region was found in the
results of the global analyses shown in Figs. 4–5.

As actual examples, we investigate the four possibilities in
$(m_{H},m_{k^{\pm\pm}})$ of $(250\,\text{GeV},250\,\text{GeV})$ in Fig. 7;
$(250\,\text{GeV},500\,\text{GeV})$ in Fig. 8; and
$(500\,\text{GeV},250\,\text{GeV})$ and $(500\,\text{GeV},500\,\text{GeV})$ in
Fig. 9, respectively. The regions covered with the color green are excluded,
where the branching ratio of the observed Higgs boson to invisible pairs
exceeds the $95\%$ confidence level upper bound $0.29$ reported by the ATLAS
experiment ATLAS:2015yda .101010In this estimation, the Standard Model
production cross section was assumed in $m_{h}=125\,\text{GeV}$. Here, we do
not care about this point seriously and simply utilize this result for putting
a bound on the present scenario. The value of $v^{\prime}$ and $\alpha$ almost
completely determines the profiles since the primary final state with
invisibleness is a pair of the NG bosons $G$. We mention that no exclusion
from the direct detection process is found in the shown parameter ranges since
the elastic scattering cross section is suppressed by the large VEV
$v^{\prime}$. We observe small differences originating from the value of
$m_{k^{\pm\pm}}$ since in the DM mass range of $m_{\chi}\lesssim
m_{k^{\pm\pm}}$, $k^{\pm\pm}$ can contribute to the relic density only
indirectly through the total widths $\Gamma_{h}$ and $\Gamma_{H}$.111111We
simply ignore the process $\chi\chi\to k^{++}k^{--}$ in the case near the
kinetic boundary, $m_{H}=500\,\text{GeV}$ and
$m_{k^{\pm\pm}}=250\,\text{GeV}$, where the heavier resonant region
corresponds to $m_{\chi}\simeq 250\,\text{GeV}$. Here, we do not consider the
lighter $k^{\pm\pm}$ with $100\sim 250\,\text{GeV}$ mass, the possibility of
which was examined in the global fits in section IV.

A suggestive aspect via the results is that only around the $m_{h}/2$ is
possible when $m_{H}$ gets large as $500\,\text{GeV}$ because the suppression
via the heavy resonant particle could eliminate the solution near $m_{H}/2$.
On the other hand, as discussed in Sec. IV, the $7$ and $8\,\text{TeV}$ LHC
results restrict the scalar mixing angle $\alpha$, where we quote a typical
bound we obtained, $|\sin{\alpha}|\lesssim 0.3$. For making the DM solution
around $m_{h}/2$ realizable, we should select the absolute value of the mixing
angle $|\sin{\alpha}|$ around or more than $0.3$ as shown in Figs. 7–9. Let us
remind the reader that the sign of $\sin\alpha$ does not affect the result
(calculated in $\sin\alpha\geq 0$) significantly. In addition, here, the value
of $v^{\prime}$ is also constrained as around or less than $9\,\text{TeV}$ in
the case of $m_{\chi}\simeq m_{h}/2$. In conclusion, the present scenario is
still viable, but a part of the parameter space generating a suitable amount
of relics comes to be on the edge. Note that our calculation would be
applicable for a rather general setup of a fermion dark matter communicating
the SM sector through two scalars with renormalizable interactions and
$\mathcal{O}(1)$ TeV $v^{\prime}$ after the decoupling of $k^{\pm\pm}$.

##  VI Summary and Conclusion

In this paper, we proposed a three loop radiative neutrino mass scenario with
an isolated doubly charged singlet scalar $k^{\pm\pm}$ without couplings to
the charged leptons at the leading order for solving an undesirable situation
in our previous study Hatanaka:2014tba . The previous scenario has the same
particle contents of charged singlet scalars with the present model, two
singly charged ones ($h_{1}^{\pm}$, $h_{2}^{\pm}$) and one doubly charged one
($k^{\pm\pm}$), but assignments of additional global symmetries are different,
which makes a difference in interactions associated with the charged bosons.

In the previous case without Majorana neutrinos, $k^{\pm\pm}$ should attach
the charged leptons at the leading order to generate a neutrino mass matrix at
the three loop level. But a side effect due to this type of couplings with LFV
could be serious from a phenomenological point of view. Because of resultant
tree-level contributions to LFV processes, related couplings are severely
constrained, where realized neutrino masses get to be smaller than the
observed values. Increasing scalar trilinear couplings with charged scalars
appearing in the numerator of elements of the loop-induced neutrino mass
matrix seems to be an option. However, the largeness in the trilinear
couplings tends to destabilize the vacuum rapidly. Making all of the charged
scalars very heavy is a solution since the scalar trilinear couplings do not
contribute to LFV at the leading order and the ratio of the trilinear
couplings divided by the masses of the charged scalars can remain “sizable”
for the observed neutrino masses. Here, however, a typical mass scale of the
charged bosons should be around $10\,\text{TeV}$, where we cannot expect a
signal suggesting the existence of the particles at the LHC.

Our modified setup in the assignment of a global $U(1)$ symmetry brings us a
better situation. Introducing Majorana neutrinos allows $k^{\pm\pm}$ to be
away from the charged leptons. Now the LFV is relatively relaxed, and at least
a part of the couplings contributing to the neutrino mass can take larger
values. Thus, scalar trilinear couplings need not to be drastically large,
and, as a result, vacuum destabilization is alleviated. We found that a few
hundred GeV $k^{\pm\pm}$ is realizable, which can be tested at the ongoing LHC
experiment, even though the two singly charged scalars are still very heavy at
around a few TeV. It is noteworthy that the associated NG boson $G$ does not
generate dangerous interactions that severely restrict the decay constant
$v^{\prime}$. The constant $v^{\prime}$ can takes a value around the TeV scale
with no harm.

Another constraint on the system comes from the mixing of two $CP$ even
neutral scalars, which are mixed states between the component of the
$SU(2)_{L}$ doublet $\Phi$ in the SM and the corresponding part of a complex
$SU(2)_{L}$ singlet scalar $\Sigma_{0}$ for the spontaneous symmetry breaking
of the global $U(1)$. After the breakdown, Majorana fermions obtain masses via
the Higgs mechanism and a residual discrete parity ensures the stability of an
accidental DM candidate of the lightest Majorana neutrino $N_{R_{1}}$. This
mixing angle plays a very important role in the DM-related phenomena in this
model since the candidate can couple to the SM sector only through the two
mass eigenstates of the $CP$ even scalars $h\,(=125\,\text{GeV})$ and $H$.
Note that the two singly charged scalars should be decoupled in our scenario
and diagrams associated with them are negligible.

The latest LHC result of the Higgs search by the ATLAS and the CMS experiments
puts a stringent bound on the mixing angle $\alpha$,
$|\sin{\alpha}|\,\lesssim\,0.3$ as a roughly estimated bound through a global
analysis. A typical scale of the VEV of $\Sigma_{0}$ would be around a few TeV
scale, where the DM pair annihilation cross section is significantly
suppressed since the Yukawa couplings between $N_{R_{1}}$ and $\Sigma_{0}$ are
inversely proportional to the VEV. When the DM mass is around the two resonant
regions, $m_{h}/2$ and $m_{H}/2$, the cross section is enhanced very much. But
we found that when $M_{H}$ is heavy at around $500\,\text{GeV}$, the
possibility near $m_{H}/2$ is closed due to the propagator suppression of the
massive mediator. On the other hand, to activate the solution around
$m_{h}/2$, the mixing angle $\alpha$ should be large to a certain degree to
ensure a sizable interaction between the DM and SM particles. Typically, a
required value in this direction is $|\sin{\alpha}|\,\gtrsim\,0.3$, which
would not have large overlaps with the preferred area obtained via the global
analysis on the Higgs search results.

In spite of the insignificance of the bounds via the invisible
$125\,\text{GeV}$ Higgs decay and the DM direct detection because of the
coupling suppression via the huge VEV of $\Sigma_{0}$, experimental data in
Higgs physics put sizable constraints on our scenario. In the near future,
updated Higgs results could declare validity of our scenario much more
precisely. Also, searching for a suitable way to discriminate the present
scenario from other models with charged scalars will be an important task at
the LHC and other future colliders.

## Acknowledgments

K. N. is grateful for Tomohiro Abe, Motoi Endo, Masahiro Ibe, Shinya Kanemura,
Tetsuo Shindou, Hiroaki Sugiyama, Koji Tsumura and Toshifumi Yamashita for
kind suggestions and useful discussions. H. O. is sincerely grateful for all
the KIAS members, Korean cordial persons, foods, culture, weather, and all the
other things. This work is supported in part by NRF Research No. 2009-0083526
(Y. O.) of the Republic of Korea.

##  Appendix A Analytic forms of loop functions for LFV

In this appendix, we show analytic forms of the loop functions for processes
with LFV. The functions for $\ell\to\ell\gamma$ in Eq. (II.42) are obtained as

$\displaystyle I^{\prime}_{1,a}$
$\displaystyle={\frac{1}{(4\pi)^{2}}}{\int_{0}^{1}\\!\\!dx\int_{0}^{1-x}\\!\\!\\!\\!\\!\\!\\!dy}\
\frac{x(2x-1)}{(x+y)m_{h_{2}^{\pm}}^{2}+(1-x-y)M_{N_{a}}^{2}}$
$\displaystyle=\begin{cases}{\frac{1}{(4\pi)^{2}}}\frac{9m_{h_{2}^{\pm}}^{4}M_{N_{a}}^{2}+9m_{h_{2}^{\pm}}^{2}M_{N_{a}}^{4}\left(2\log\left(\frac{M_{N_{a}}^{2}}{m_{h_{2}^{\pm}}^{2}}\right)+1\right)+M_{N_{a}}^{6}\left(6\log\left(\frac{M_{N_{a}}^{2}}{m_{h_{2}^{\pm}}^{2}}\right)-17\right)-m_{h_{2}^{\pm}}^{6}}{36\left(m_{h_{2}^{\pm}}^{2}-M_{N_{a}}^{2}\right)^{4}}&(\text{for
}m_{h_{2}^{\pm}}\not=M_{N_{a}}),\\\ 0&(\text{for
}m_{h_{2}^{\pm}}=M_{N_{a}}),\end{cases}$ (A.1) $\displaystyle
I^{\prime}_{2,a}$
$\displaystyle={\frac{1}{(4\pi)^{2}}}{\int_{0}^{1}\\!\\!dx\int_{0}^{1-x}\\!\\!\\!\\!\\!\\!\\!dy}\
\frac{x(2y-1)}{(x+y)m_{h_{2}^{\pm}}^{2}+(1-x-y)M_{N_{a}}^{2}}$
$\displaystyle=\begin{cases}{\frac{1}{(4\pi)^{2}}}\frac{27m_{h_{2}^{\pm}}^{4}M_{N_{a}}^{2}-9m_{h_{2}^{\pm}}^{2}M_{N_{a}}^{4}\left(2\log\left(\frac{m_{h_{2}^{\pm}}^{2}}{M_{N_{a}}^{2}}\right)+3\right)+M_{N_{a}}^{6}\left(6\log\left(\frac{m_{h_{2}^{\pm}}^{2}}{M_{N_{a}}^{2}}\right)+5\right)-5m_{h_{2}^{\pm}}^{6}}{36\left(m_{h_{2}^{\pm}}^{2}-M_{N_{a}}^{2}\right)^{4}}&(\text{for
}m_{h_{2}^{\pm}}\not=M_{N_{a}}),\\\
-{\frac{1}{(4\pi)^{2}}}\frac{1}{12m_{h_{2}^{\pm}}^{2}}&(\text{for
}m_{h_{2}^{\pm}}=M_{N_{a}}).\end{cases}$ (A.2)

Note that when we take the limit $M_{N_{a}}\to 0$, $I^{\prime}_{1,a}$ and
$I^{\prime}_{2,a}$ are reduced to $I_{1,a}$ and $I_{2,a}$ (up to the
difference of $m_{h_{1}^{\pm}}$ and $m_{h_{2}^{\pm}}$), respectively,

$\displaystyle
I^{\prime}_{1,a}\to-{\frac{1}{(4\pi)^{2}}}\frac{1}{36m_{h_{2}^{\pm}}^{2}},\quad
I^{\prime}_{2,a}\to-{\frac{1}{(4\pi)^{2}}}\frac{5}{36m_{h_{2}^{\pm}}^{2}}.$
(A.3)

The functions $J_{1,ab}$ and $J_{2,ab}$ for describing the $\ell\to 3\ell$
processes in Eq. (II.48) take the following forms

$\displaystyle J_{1,ab}$
$\displaystyle=\frac{1}{(4\pi)^{2}}{\int_{0}^{1}\\!\\!dx\int_{0}^{1-x}\\!\\!\\!\\!\\!\\!\\!dy}\
\frac{1-x-y}{\left[xM_{N_{a}}^{2}+yM_{N_{b}}^{2}+(1-x-y)m_{h_{2}^{\pm}}^{2}\right]}$
$\displaystyle={\scriptstyle{\frac{1}{(4\pi)^{2}}\bigg{\\{}\frac{m_{h_{2}^{\pm}}^{6}M_{N_{a}}^{2}+m_{h_{2}}^{4}M_{N_{a}}^{4}\left(\log\left(\frac{M_{N_{a}}^{2}}{m_{h_{2}^{\pm}}^{2}}\right)-1\right)-m_{h_{2}^{\pm}}^{6}M_{N_{b}}^{2}+m_{h_{2}^{\pm}}^{2}M_{N_{a}}^{4}M_{N_{b}}^{2}\left(2\log\left(\frac{m_{h_{2}^{\pm}}^{2}}{M_{N_{a}}^{2}}\right)+1\right)}{2\left(M_{N_{a}}^{2}-M_{N_{b}}^{2}\right)\left(m_{h_{2}^{\pm}}^{2}-M_{N_{a}}^{2}\right)^{2}\left(m_{h_{2}^{\pm}}^{2}-M_{N_{b}}^{2}\right)^{2}}}}$
$\displaystyle\qquad{\scriptstyle{+\frac{m_{h_{2}^{\pm}}^{4}M_{N_{b}}^{4}\left(\log\left(\frac{m_{h_{2}^{\pm}}^{2}}{M_{N_{b}}^{2}}\right)+1\right)+m_{h_{2}^{\pm}}^{2}M_{N_{a}}^{2}M_{N_{b}}^{4}\left(2\log\left(\frac{M_{N_{b}}^{2}}{m_{h_{2}^{\pm}}}\right)-1\right)+M_{N_{a}}^{4}M_{N_{b}}^{4}\log\left(\frac{M_{N_{a}}^{2}}{M_{N_{b}}^{2}}\right)}{2\left(M_{N_{a}}^{2}-M_{N_{b}}^{2}\right)\left(m_{h_{2}^{\pm}}^{2}-M_{N_{a}}^{2}\right)^{2}\left(m_{h_{2}^{\pm}}^{2}-M_{N_{b}}^{2}\right)^{2}}\bigg{\\}}}}$
$\displaystyle\hskip 261.76535pt(\text{for
}m_{h_{2}^{\pm}}\not=M_{N_{a}}\not=M_{N_{b}}),$ (A.4)
$\displaystyle{J_{2,ab}}={\frac{1}{(4\pi)^{2}}}{\int_{0}^{1}\\!\\!dx\int_{0}^{1-x}\\!\\!\\!\\!\\!\\!\\!dy}\
\frac{1-x-y}{\left[xM_{N_{a}}^{2}+yM_{N_{b}}^{2}+(1-x-y)m_{h_{2}^{\pm}}^{2}\right]^{2}}$
$\displaystyle={\scriptstyle\frac{m_{h_{2}^{\pm}}^{4}\left(M_{N_{a}}^{2}\left(\log\left(\frac{m_{h_{2}^{\pm}}^{2}}{M_{N_{a}}^{2}}\right)-1\right)+M_{N_{b}}^{2}\left(\log\left(\frac{M_{N_{b}}^{2}}{m_{h_{2}^{\pm}}^{2}}\right)+1\right)\right)+m_{h_{2}^{\pm}}^{2}\left(2M_{N_{a}}^{2}M_{N_{b}}^{2}\log\left(\frac{M_{N_{a}}^{2}}{M_{N_{b}}^{2}}\right)+M_{N_{a}}^{4}-M_{N_{b}}^{4}\right)+M_{N_{a}}^{2}M_{N_{b}}^{4}\left(\log\left(\frac{m_{h_{2}^{\pm}}^{2}}{M_{N_{a}}^{2}}\right)+1\right)+M_{N_{a}}^{4}M_{N_{b}}^{2}\left(\log\left(\frac{M_{N_{b}}^{2}}{m_{h_{2}^{\pm}}^{2}}\right)-1\right)}{{(4\pi)^{2}}\left(m_{h_{2}^{\pm}}^{2}-M_{N_{a}}^{2}\right)^{2}\left(m_{h_{2}^{\pm}}^{2}-M_{N_{b}}^{2}\right)^{2}\left(M_{N_{a}}^{2}-M_{N_{b}}^{2}\right)}}$
$\displaystyle\hskip 284.52756pt(\text{for
}m_{h_{2}^{\pm}}\not=M_{N_{a}}\not=M_{N_{b}}),$ (A.5)

where the following specific cases are also obtained:

$\displaystyle J_{1,ab}$
$\displaystyle\to\begin{cases}\frac{1}{(4\pi)^{2}}\frac{2m_{h_{2}^{\pm}}^{2}M_{N_{a}}^{2}\log\left(\frac{M_{N_{a}}^{2}}{m_{h_{2}^{\pm}}^{2}}\right)-M_{N_{a}}^{4}+m_{h_{2}^{\pm}}^{4}}{2\left(m_{h_{2}^{\pm}}^{2}-M_{N_{a}}^{2}\right)^{3}}&(\text{for
}M_{N_{a}}=M_{N_{b}}\not=m_{h_{2}^{\pm}}),\\\
\frac{1}{(4\pi)^{2}}\frac{1}{2m_{h_{2}^{\pm}}^{2}}&(\text{for
}M_{N_{a}}=M_{N_{b}}=m_{h_{2}^{\pm}}),\end{cases}$ (A.6)
$\displaystyle{J_{2,ab}}$
$\displaystyle\to\begin{cases}{\frac{1}{(4\pi)^{2}}}\frac{m_{h_{2}^{\pm}}^{2}\left(\log\left(\frac{m_{h_{2}^{\pm}}^{2}}{M_{N_{a}}^{2}}\right)-2\right)+M_{N_{a}}^{2}\left(\log\left(\frac{m_{h_{2}^{\pm}}^{2}}{M_{N_{a}}^{2}}\right)+2\right)}{\left(m_{h_{2}^{\pm}}^{2}-M_{N_{a}}^{2}\right)^{3}}&(\text{for
}M_{N_{a}}=M_{N_{b}}\not=m_{h_{2}^{\pm}}),\\\
{\frac{1}{(4\pi)^{2}}}\frac{1}{6m_{h_{2}^{\pm}}^{4}}&(\text{for
}M_{N_{a}}=M_{N_{b}}=m_{h_{2}^{\pm}}).\end{cases}$ (A.7)

##  Appendix B Partial widths of CP even scalars

In this appendix, we represent a part of partial decay widths with nontrivial
forms of the two $CP$ even scalars $h$ and $H$ with our notation for loop
functions.

$\displaystyle\Gamma_{h\to GG}$ $\displaystyle=\frac{m_{h}^{3}}{{32}\pi
v^{\prime 2}}\sin^{2}{\alpha},$ (B.1) $\displaystyle\Gamma_{h\to
N_{R_{1}}\bar{N}_{R_{1}}}$
$\displaystyle=\frac{1}{{16}\pi}\left(\frac{M_{N_{1}}}{v^{\prime}}\right)^{2}m_{h}\left(1-\frac{4M_{N_{1}}^{2}}{m_{h}^{2}}\right)^{3/2}\sin^{2}{\alpha},$
(B.2) $\displaystyle\Gamma_{H\to GG}$ $\displaystyle=\frac{m_{H}^{3}}{{32}\pi
v^{\prime 2}}\cos^{2}{\alpha},$ (B.3) $\displaystyle\Gamma_{H\to
N_{R_{1}}\bar{N}_{R_{1}}}$
$\displaystyle=\frac{1}{{16}\pi}\left(\frac{M_{N_{1}}}{v^{\prime}}\right)^{2}m_{H}\left(1-\frac{4M_{N_{1}}^{2}}{m_{H}^{2}}\right)^{3/2}\cos^{2}{\alpha},$
(B.4) $\displaystyle\Gamma_{H\to 2h}$ $\displaystyle=\frac{1}{{32}\pi
m_{H}}|C_{Hhh}|^{2}\sqrt{1-\frac{4m_{h}^{2}}{m_{H}^{2}}},$ (B.5)
$\displaystyle\Gamma_{h\to gg}$
$\displaystyle=\frac{G_{F}\alpha_{s}^{2}m_{h}^{3}}{36\sqrt{2}\pi^{3}}\left|\frac{3}{4}\left(A_{1/2}^{\gamma\gamma}(\tau_{t})\right)\cos{\alpha}\right|^{2},$
(B.6) $\displaystyle\Gamma_{h\to\gamma\gamma}$
$\displaystyle=\frac{G_{F}\alpha_{\text{EM}}^{2}m_{h}^{3}}{128\sqrt{2}\pi^{3}}\left|\left(A_{1}^{\gamma\gamma}(\tau_{W})+N_{C}Q_{t}^{2}A_{1/2}^{\gamma\gamma}(\tau_{t})\right)\cos{\alpha}+\frac{1}{2}\frac{v^{2}}{m_{k^{\pm\pm}}^{2}}Q_{k}^{2}C_{hkk}A_{0}^{\gamma\gamma}(\tau_{k})\right|^{2},$
(B.7) $\displaystyle\Gamma_{h\to Z\gamma}$
$\displaystyle=\frac{{\alpha_{\text{EM}}^{2}}m_{h}^{3}}{512\pi^{3}}\left(1-\frac{m_{Z}^{2}}{m_{h}^{2}}\right)^{3}\left|{A_{\text{SM}}^{Z\gamma}}\cos{\alpha}-\frac{v}{m_{k^{\pm\pm}}^{2}}\left(2Q_{k}g_{Zkk}\right)C_{hkk}A_{0}^{Z\gamma}(\tau_{k},\lambda_{k})\right|^{2},$
(B.8) $\displaystyle\Gamma_{H\to gg}$
$\displaystyle=\frac{G_{F}\alpha_{s}^{2}m_{H}^{3}}{36\sqrt{2}\pi^{3}}\left|\frac{3}{4}\left(A_{1/2}^{\gamma\gamma}(\tau^{\prime}_{t})\right)\sin{\alpha}\right|^{2},$
(B.9) $\displaystyle\Gamma_{H\to\gamma\gamma}$
$\displaystyle=\frac{G_{F}\alpha_{\text{EM}}^{2}m_{H}^{3}}{128\sqrt{2}\pi^{3}}\left|\left(A_{1}^{\gamma\gamma}(\tau^{\prime}_{W})+N_{C}Q_{t}^{2}A_{1/2}^{\gamma\gamma}(\tau^{\prime}_{t})\right)\sin{\alpha}+\frac{1}{2}\frac{v^{2}}{m_{k^{\pm\pm}}^{2}}Q_{k}^{2}C_{Hkk}A_{0}^{\gamma\gamma}(\tau^{\prime}_{k})\right|^{2},$
(B.10) $\displaystyle\Gamma_{H\to Z\gamma}$
$\displaystyle=\frac{{\alpha_{\text{EM}}^{2}}m_{H}^{3}}{512\pi^{3}}\left(1-\frac{m_{Z}^{2}}{m_{H}^{2}}\right)^{3}\left|{A_{\text{SM}}^{Z\gamma}}(\tau_{W,t}\to\tau^{\prime}_{W,t})\sin{\alpha}-\frac{v}{m_{k^{\pm\pm}}^{2}}\left(2Q_{k}g_{Zkk}\right)C_{Hkk}A_{0}^{Z\gamma}(\tau^{\prime}_{k},\lambda_{k})\right|^{2},$
(B.11)

with $C_{hkk}$ in Eq. (IV.9) and the parameters

$\displaystyle\tau^{\prime}_{i}=\frac{4m_{i}^{2}}{m_{H}^{2}}\ (i=t,W,k),\quad
C_{Hkk}=\sin{\alpha}\lambda_{\Phi k}+\cos{\alpha}\lambda_{\Sigma
k}\left(\frac{v^{\prime}}{v}\right).$ (B.12)

$G_{F}$ and $\alpha_{s}$ are the Fermi constant and the QCD coupling strength,
respectively. Note that in the calculation $\Gamma_{H\to Z\gamma}$, the
replacement $\tau_{i}\to\tau^{\prime}_{i}$ in the factor
$A_{\text{SM}}^{Z\gamma}$ is required. One can refer to Sec. IV for the detail
of loop functions.

##  Appendix C Matrix elements for $N_{R_{1}}$ dark matter annihilations

In this appendix, we summarize the averaged matrix elements squared for relic
density calculation of the dark matter candidate $N_{R_{1}}$ through pair
annihilations. Note that no symmetric factor is included originating from
identical final states in the following forms, which should be considered in
the integration over the solid angle for estimating the thermal average of the
cross section $\langle\sigma v_{\text{rel}}\rangle$.

$\displaystyle|\overline{\mathcal{M}}(N_{R_{1}}\bar{N}_{R_{1}}\to
f\bar{f})|^{2}$
$\displaystyle=4{N_{C_{f}}}\left(\frac{M_{N_{1}}m_{f}s_{\alpha}c_{\alpha}}{vv^{\prime}}\right)^{2}\left|\frac{1}{s-m_{h}^{2}+im_{h}\Gamma_{h}}-\frac{1}{s-m_{H}^{2}+im_{H}\Gamma_{H}}\right|^{2}\times$
$\displaystyle\quad\left[(p_{1}\cdot
p_{2})-M_{N_{1}}^{2}\right]\left[(k_{1}\cdot k_{2})-m_{f}^{2}\right],$ (C.1)
$\displaystyle|\overline{\mathcal{M}}(N_{R_{1}}\bar{N}_{R_{1}}\to ZZ)|^{2}$
$\displaystyle=\left(\frac{M_{N_{1}}g_{2}M_{W}s_{\alpha}c_{\alpha}}{v^{\prime}c_{W}^{2}}\right)^{2}\left|\frac{1}{s-m_{h}^{2}+im_{h}\Gamma_{h}}-\frac{1}{s-m_{H}^{2}+im_{H}\Gamma_{H}}\right|^{2}\times$
$\displaystyle\quad\left[(p_{1}\cdot
p_{2})-M_{N_{1}}^{2}\right]\left[2+\frac{(k_{1}\cdot
k_{2})^{2}}{M_{Z}^{4}}\right],$ (C.2)
$\displaystyle|\overline{\mathcal{M}}(N_{R_{1}}\bar{N}_{R_{1}}\to
W^{+}W^{-})|^{2}$
$\displaystyle=\left(\frac{M_{N_{1}}g_{2}M_{W}s_{\alpha}c_{\alpha}}{v^{\prime}}\right)^{2}\left|\frac{1}{s-m_{h}^{2}+im_{h}\Gamma_{h}}-\frac{1}{s-m_{H}^{2}+im_{H}\Gamma_{H}}\right|^{2}\times$
$\displaystyle\quad\left[(p_{1}\cdot
p_{2})-M_{N_{1}}^{2}\right]\left[2+\frac{(k_{1}\cdot
k_{2})^{2}}{M_{W}^{4}}\right],$ (C.3)
$\displaystyle|\overline{\mathcal{M}}(N_{R_{1}}\bar{N}_{R_{1}}\to GG)|^{2}$
$\displaystyle=4\left(\frac{M_{N_{1}}}{vv^{\prime}}\right)^{2}(k_{1}\cdot
k_{2})^{2}\left|\frac{s_{\alpha}^{2}}{s-m_{h}^{2}+im_{h}\Gamma_{h}}+\frac{c_{\alpha}^{2}}{s-m_{H}^{2}+im_{H}\Gamma_{H}}\right|^{2}\times$
$\displaystyle\quad\left[(p_{1}\cdot p_{2})-M_{N_{1}}^{2}\right],$ (C.4)
$\displaystyle|\overline{\mathcal{M}}(N_{R_{1}}\bar{N}_{R_{1}}\to hh)|^{2}$
$\displaystyle=\left(\frac{M_{N_{1}}}{v^{\prime}}\right)^{2}\left|\frac{s_{\alpha}C_{hhh}}{s-m_{h}^{2}+im_{h}\Gamma_{h}}-\frac{c_{\alpha}C_{Hhh}}{s-m_{H}^{2}+im_{H}\Gamma_{H}}\right|^{2}\times$
$\displaystyle\quad\left[(p_{1}\cdot p_{2})-M_{N_{1}}^{2}\right],$ (C.5)

where $N_{C_{f}}$ shows the color factor $3$ ($1$) when $f$ is a quark
(lepton), $p_{1},p_{2}$ and $k_{1},k_{2}$ are initial- and final-state
momenta, respectively. $s$ is defined as
$s=(p_{1}+p_{2})^{2}=(k_{1}+k_{2})^{2}$. Here, we use the shorthand notations
of $c_{\alpha}\equiv\cos{\alpha}$ and $s_{\alpha}\equiv\sin{\alpha}$. The
forms of effective couplings describing scalar self interactions are

$\displaystyle C_{hhh}$
$\displaystyle=3\lambda_{\Phi\Sigma}\left(v(s_{\alpha}^{2}c_{\alpha})-v^{\prime}(c_{\alpha}^{2}s_{\alpha})\right){+6\lambda_{\Phi}vc_{\alpha}^{3}-6\lambda_{\Sigma}v^{\prime}s_{\alpha}^{3}},$
(C.6) $\displaystyle C_{Hhh}$
$\displaystyle=\lambda_{\Phi\Sigma}{\left(v\left(s^{3}_{\alpha}-2c^{2}_{\alpha}s_{\alpha}\right)+v^{\prime}\left(c_{\alpha}^{3}-2c_{\alpha}s^{2}_{\alpha}\right)\right)+6\lambda_{\Phi}vc_{\alpha}^{2}s_{\alpha}+6\lambda_{\Sigma}v^{\prime}s_{\alpha}^{2}c_{\alpha}}.$
(C.7)

We mention that the processes $N_{R_{1}}\bar{N}_{R_{1}}$ into $HH$, $hH$,
$k^{++}k^{--}$, $h_{1}^{+}h_{1}^{-}$ and $h_{2}^{+}h_{2}^{-}$ are ineffective
in our discussion when the mass of the DM is around $m_{h}/2$ or $m_{H}/2$
considering the bounds discussed in Sec. III–IV.

Note that the processes
$N_{R_{1}}\bar{N}_{R_{1}}\to\ell^{+}_{i}\ell^{-}_{j}\,(i,j=1,2,3)$ with the
exchange of $h^{\pm}_{2}$ in the $t$ channel, which are important in similar
situations Krauss:2002px ; Cheung:2004xm , are ineffective and neglected in
the present scenario because the lower bound on $m_{h^{\pm}_{2}}$ is around a
few TeV and decoupled as shown in Sec. III.2. We neglect the three- and four-
body final states via virtual $W$ and $Z$ boson decays, which gives sizable
modifications near the thresholds for producing gauge boson pairs Cline:2012hg
; Cline:2013gha , since our interest is only around $m_{h}/2$ and $m_{H}/2$
where this correction is expected to be ineffective.

## References

  * (1) H. Hatanaka, K. Nishiwaki, H. Okada and Y. Orikasa, Nucl. Phys. B 894, 268 (2015) [arXiv:1412.8664 [hep-ph]].
  * (2) G. Aad et al. [ATLAS Collaboration], Phys. Lett. B 716, 1 (2012) [arXiv:1207.7214 [hep-ex]].
  * (3) S. Chatrchyan et al. [CMS Collaboration], Phys. Lett. B 716, 30 (2012) [arXiv:1207.7235 [hep-ex]].
  * (4) A. Zee, Phys. Lett. B 93, 389 (1980) [Erratum-ibid. B 95, 461 (1980)].
  * (5) T. P. Cheng and L. F. Li, Phys. Rev. D 22, 2860 (1980).
  * (6) A. Zee, Nucl. Phys. B 264, 99 (1986); K. S. Babu, Phys. Lett. B 203, 132 (1988).
  * (7) L. M. Krauss, S. Nasri and M. Trodden, Phys. Rev. D 67, 085002 (2003) [arXiv:hep-ph/0210389].
  * (8) E. Ma, Phys. Rev. D 73, 077301 (2006) [hep-ph/0601225].
  * (9) T. Hambye, K. Kannike, E. Ma and M. Raidal, Phys. Rev. D 75, 095003 (2007) [hep-ph/0609228].
  * (10) P. -H. Gu and U. Sarkar, Phys. Rev. D 77, 105031 (2008) [arXiv:0712.2933 [hep-ph]].
  * (11) N. Sahu and U. Sarkar, Phys. Rev. D 78, 115013 (2008) [arXiv:0804.2072 [hep-ph]].
  * (12) P. -H. Gu and U. Sarkar, Phys. Rev. D 78, 073012 (2008) [arXiv:0807.0270 [hep-ph]].
  * (13) M. Aoki, S. Kanemura and O. Seto, Phys. Rev. Lett. 102, 051805 (2009) [arXiv:0807.0361].
  * (14) J. March-Russell, C. McCabe and M. McCullough, JHEP 1003, 108 (2010) [arXiv:0911.4489 [hep-ph]].
  * (15) M. Aoki, S. Kanemura, T. Shindou and K. Yagyu, JHEP 1007, 084 (2010) [Erratum-ibid. 1011, 049 (2010)] [arXiv:1005.5159 [hep-ph]].
  * (16) S. Kanemura, O. Seto and T. Shimomura, Phys. Rev. D 84, 016004 (2011) [arXiv:1101.5713 [hep-ph]].
  * (17) M. Lindner, D. Schmidt and T. Schwetz, Phys. Lett. B 705, 324 (2011) [arXiv:1105.4626 [hep-ph]].
  * (18) S. Kanemura, T. Nabeshima and H. Sugiyama, Phys. Lett. B 703, 66 (2011) [arXiv:1106.2480 [hep-ph]].
  * (19) M. Aoki, J. Kubo, T. Okawa and H. Takano, Phys. Lett. B 707, 107 (2012) [arXiv:1110.5403 [hep-ph]].
  * (20) S. Kanemura, T. Nabeshima and H. Sugiyama, Phys. Rev. D 85, 033004 (2012) [arXiv:1111.0599 [hep-ph]].
  * (21) D. Schmidt, T. Schwetz and T. Toma, Phys. Rev. D 85, 073009 (2012) [arXiv:1201.0906 [hep-ph]].
  * (22) S. Kanemura and H. Sugiyama, Phys. Rev. D 86, 073006 (2012) [arXiv:1202.5231 [hep-ph]].
  * (23) Y. Farzan and E. Ma, Phys. Rev. D 86, 033007 (2012) [arXiv:1204.4890 [hep-ph]].
  * (24) K. Kumericki, I. Picek and B. Radovcic, JHEP 1207, 039 (2012) [arXiv:1204.6597 [hep-ph]].
  * (25) K. Kumericki, I. Picek and B. Radovcic, Phys. Rev. D 86, 013006 (2012) [arXiv:1204.6599 [hep-ph]].
  * (26) E. Ma, Phys. Lett. B 717, 235 (2012) [arXiv:1206.1812 [hep-ph]].
  * (27) G. Gil, P. Chankowski and M. Krawczyk, Phys. Lett. B 717, 396 (2012) [arXiv:1207.0084 [hep-ph]].
  * (28) H. Okada and T. Toma, Phys. Rev. D 86, 033011 (2012) arXiv:1207.0864 [hep-ph].
  * (29) D. Hehn and A. Ibarra, Phys. Lett. B 718, 988 (2013) [arXiv:1208.3162 [hep-ph]].
  * (30) S. Baek, P. Ko, H. Okada and E. Senaha, JHEP 1409, 153 (2014) [arXiv:1209.1685 [hep-ph]].
  * (31) P. S. B. Dev and A. Pilaftsis, Phys. Rev. D 86, 113001 (2012) [arXiv:1209.4051 [hep-ph]].
  * (32) Y. Kajiyama, H. Okada and T. Toma, Eur. Phys. J. C 73, 2381 (2013) [arXiv:1210.2305 [hep-ph]].
  * (33) M. Kohda, H. Sugiyama and K. Tsumura, Phys. Lett. B 718, 1436 (2013) [arXiv:1210.5622 [hep-ph]].
  * (34) H. Okada, arXiv:1212.0492 [hep-ph].
  * (35) M. Gustafsson, J. M. No and M. A. Rivera, Phys. Rev. Lett. 110, 211802 (2013) arXiv:1212.4806 [hep-ph].
  * (36) M. Aoki, J. Kubo and H. Takano, Phys. Rev. D 87, no. 11, 116001 (2013) [arXiv:1302.3936 [hep-ph]].
  * (37) Y. Kajiyama, H. Okada and K. Yagyu, Nucl. Phys. B 874, 198 (2013) [arXiv:1303.3463 [hep-ph]].
  * (38) Y. Kajiyama, H. Okada and T. Toma, Phys. Rev. D 88, 015029 (2013) [arXiv:1303.7356].
  * (39) S. Kanemura, T. Matsui and H. Sugiyama, Phys. Lett. B 727, 151 (2013) [arXiv:1305.4521 [hep-ph]].
  * (40) B. Dasgupta, E. Ma and K. Tsumura, Phys. Rev. D 89, 041702 (2014) [arXiv:1308.4138 [hep-ph]].
  * (41) S. Baek, H. Okada and T. Toma, JCAP 1406, 027 (2014) [arXiv:1312.3761 [hep-ph]].
  * (42) H. Okada, arXiv:1404.0280 [hep-ph].
  * (43) A. Ahriche, C. S. Chen, K. L. McDonald and S. Nasri, Phys. Rev. D 90, no. 1, 015024 (2014) [arXiv:1404.2696 [hep-ph]].
  * (44) A. Ahriche, K. L. McDonald and S. Nasri, JHEP 1410, 167 (2014) [arXiv:1404.5917 [hep-ph]].
  * (45) C. -S. Chen, K. L. McDonald and S. Nasri, Phys. Lett. B 734, 388 (2014) [arXiv:1404.6033 [hep-ph]].
  * (46) S. Kanemura, T. Matsui and H. Sugiyama, Phys. Rev. D 90, 013001 (2014) [arXiv:1405.1935 [hep-ph]].
  * (47) H. Okada and Y. Orikasa, Phys. Rev. D 90, no. 7, 075023 (2014) [arXiv:1407.2543 [hep-ph]].
  * (48) S. Fraser, E. Ma and O. Popov, Phys. Lett. B 737, 280 (2014) [arXiv:1408.4785 [hep-ph]].
  * (49) H. Okada, T. Toma and K. Yagyu, Phys. Rev. D 90, no. 9, 095005 (2014) [arXiv:1408.0961 [hep-ph]].
  * (50) S. Baek, H. Okada and K. Yagyu, JHEP 1504, 049 (2015) [arXiv:1501.01530 [hep-ph]].
  * (51) L. G. Jin, R. Tang and F. Zhang, Phys. Lett. B 741, 163 (2015) [arXiv:1501.02020 [hep-ph]].
  * (52) P. Culjak, K. Kumericki and I. Picek, Phys. Lett. B 744, 237 (2015) [arXiv:1502.07887 [hep-ph]].
  * (53) H. Okada, arXiv:1503.04557 [hep-ph].
  * (54) C. Q. Geng and L. H. Tsai, arXiv:1503.06987 [hep-ph].
  * (55) H. Okada, N. Okada and Y. Orikasa, arXiv:1504.01204 [hep-ph].
  * (56) C. Q. Geng, D. Huang and L. H. Tsai, Phys. Lett. B 745, 56 (2015) [arXiv:1504.05468 [hep-ph]].
  * (57) A. Ahriche, K. L. McDonald, S. Nasri and T. Toma, Phys. Lett. B 746, 430 (2015) [arXiv:1504.05755 [hep-ph]].
  * (58) D. Restrepo, A. Rivera, M. Sánchez-Peláez, O. Zapata and W. Tangarife, arXiv:1504.07892 [hep-ph].
  * (59) S. Kashiwase, H. Okada, Y. Orikasa and T. Toma, arXiv:1505.04665 [hep-ph].
  * (60) Y. H. Ahn and H. Okada, Phys. Rev. D 85, 073010 (2012) [arXiv:1201.4436 [hep-ph]].
  * (61) E. Ma, A. Natale and A. Rashed, Int. J. Mod. Phys. A 27, 1250134 (2012) [arXiv:1206.1570 [hep-ph]].
  * (62) Y. Kajiyama, H. Okada and K. Yagyu, JHEP 10, 196 (2013) arXiv:1307.0480 [hep-ph].
  * (63) A. E. Carcamo Hernandez, I. d. M. Varzielas, S. G. Kovalenko, H. Päs and I. Schmidt, Phys. Rev. D 88, 076014 (2013) [arXiv:1307.6499 [hep-ph]].
  * (64) E. Ma and A. Natale, Phys. Lett. B 723, 403 (2014) [arXiv:1403.6772 [hep-ph]].
  * (65) M. Aoki and T. Toma, JCAP 1409, 016 (2014) [arXiv:1405.5870 [hep-ph]].
  * (66) E. Ma, Phys. Lett. B 741, 202 (2015) [arXiv:1411.6679 [hep-ph]].
  * (67) E. Ma, arXiv:1504.02086 [hep-ph].
  * (68) A. E. Carcamo Hernandez, RMartinez and F. Ochoa, arXiv:1309.6567 [hep-ph].
  * (69) E. Ma, Phys. Rev. Lett. 112, 091801 (2014) [arXiv:1311.3213 [hep-ph]].
  * (70) H. Okada and K. Yagyu, Phys. Rev. D 89, 053008 (2014) [arXiv:1311.4360 [hep-ph]].
  * (71) S. Baek, H. Okada and T. Toma, Phys. Lett. B 732, 85 (2014) [arXiv:1401.6921 [hep-ph]].
  * (72) H. Okada and K. Yagyu, Phys. Rev. D 90, no. 3, 035019 (2014) [arXiv:1405.2368 [hep-ph]].
  * (73) F. Bonnet, M. Hirsch, T. Ota and W. Winter, JHEP 1207, 153 (2012) [arXiv:1204.5862 [hep-ph]].
  * (74) D. Aristizabal Sierra, A. Degee, L. Dorame and M. Hirsch, JHEP 1503, 040 (2015) [arXiv:1411.7038 [hep-ph]].
  * (75) H. Davoudiasl and I. M. Lewis, Phys. Rev. D 90, no. 3, 033003 (2014) [arXiv:1404.6260 [hep-ph]].
  * (76) M. Lindner, S. Schmidt and J. Smirnov, arXiv:1405.6204 [hep-ph].
  * (77) H. Okada and Y. Orikasa, arXiv:1412.3616 [hep-ph].
  * (78) K. S. Babu and C. Macesanu, Phys. Rev. D 67, 073010 (2003) [hep-ph/0212058].
  * (79) D. Aristizabal Sierra and D. Restrepo, JHEP 0608, 036 (2006) [hep-ph/0604012].
  * (80) D. Aristizabal Sierra and M. Hirsch, JHEP 0612, 052 (2006) [hep-ph/0609307].
  * (81) M. Nebot, J. F. Oliver, D. Palao and A. Santamaria, Phys. Rev. D 77, 093013 (2008) [arXiv:0711.0483 [hep-ph]].
  * (82) R. Bouchand and A. Merle, JHEP 1207, 084 (2012) [arXiv:1205.0008 [hep-ph]].
  * (83) Y. Kajiyama, H. Okada and K. Yagyu, Nucl. Phys. B 887, 358 (2014) [arXiv:1309.6234 [hep-ph]].
  * (84) K. L. McDonald, JHEP 1311, 131 (2013) [arXiv:1310.0609 [hep-ph]].
  * (85) E. Ma, Phys. Lett. B 732, 167 (2014) [arXiv:1401.3284 [hep-ph]].
  * (86) D. Schmidt, T. Schwetz and H. Zhang, Nucl. Phys. B 885, 524 (2014) [arXiv:1402.2251 [hep-ph]].
  * (87) J. Herrero-Garcia, M. Nebot, N. Rius and A. Santamaria, Nucl. Phys. B 885, 542 (2014) [arXiv:1402.4491 [hep-ph]].
  * (88) A. Ahriche, S. Nasri and R. Soualah, Phys. Rev. D 89, 095010 (2014) [arXiv:1403.5694 [hep-ph]].
  * (89) H. N. Long and V. V. Vien, Int. J. Mod. Phys. A 29, no. 13, 1450072 (2014) [arXiv:1405.1622 [hep-ph]].
  * (90) V. Van Vien, H. N. Long and P. N. Thu, arXiv:1407.8286 [hep-ph].
  * (91) S. Weinberg, Phys. Rev. Lett. 110, no. 24, 241301 (2013) [arXiv:1305.1971 [astro-ph.CO]].
  * (92) Y. Chikashige, R. N. Mohapatra and R. D. Peccei, Phys. Lett. B 98, 265 (1981); A. Adulpravitchai, P. H. Gu and M. Lindner, Phys. Rev. D 82, 073013 (2010).
  * (93) D. Cadamuro and J. Redondo, JCAP 1202, 032 (2012) [arXiv:1110.2895 [hep-ph]].
  * (94) A. Pilaftsis, Phys. Rev. D 49, 2398 (1994) [hep-ph/9308258].
  * (95) K. Cheung, W. Y. Keung and T. C. Yuan, Phys. Rev. D 89, no. 1, 015007 (2014) [arXiv:1308.4235 [hep-ph]].
  * (96) R. Battye and A. Moss, Phys. Rev. D 82, 023521 (2010) [arXiv:1005.0479 [astro-ph.CO]].
  * (97) Z. Maki, M. Nakagawa and S. Sakata, Prog. Theor. Phys. 28, 870 (1962).
  * (98) B. Pontecorvo, Sov. Phys. JETP 26, 984 (1968) [Zh. Eksp. Teor. Fiz. 53, 1717 (1967)].
  * (99) M. C. Gonzalez-Garcia, M. Maltoni and T. Schwetz, JHEP 1411, 052 (2014) [arXiv:1409.5439 [hep-ph]].
  * (100) E. Arganda and M. J. Herrero, Phys. Rev. D 73, 055003 (2006) [hep-ph/0510405].
  * (101) T. Toma and A. Vicente, JHEP 1401, 160 (2014) [arXiv:1312.2840 [hep-ph]].
  * (102) T. A. Chowdhury and S. Nasri, arXiv:1506.00261 [hep-ph].
  * (103) G. Aad et al. [ATLAS Collaboration], JHEP 1503, 041 (2015) [arXiv:1412.0237 [hep-ex]].
  * (104) P. P. Giardino, K. Kannike, M. Raidal and A. Strumia, Phys. Lett. B 718, 469 (2012) [arXiv:1207.1347 [hep-ph]].
  * (105) D. Carmi, A. Falkowski, E. Kuflik, T. Volansky and J. Zupan, JHEP 1210, 196 (2012) [arXiv:1207.1718 [hep-ph]].
  * (106) J. Ellis and T. You, JHEP 1209, 123 (2012) [arXiv:1207.1693 [hep-ph]].
  * (107) G. Aad et al. [ATLAS Collaboration], Phys. Rev. D 90, no. 11, 112015 (2014) [arXiv:1408.7084 [hep-ex]].
  * (108) G. Aad et al. [ATLAS Collaboration], Phys. Rev. D 91, no. 1, 012006 (2015) [arXiv:1408.5191 [hep-ex]].
  * (109) G. Aad et al. [ATLAS Collaboration], arXiv:1412.2641 [hep-ex].
  * (110) G. Aad et al. [ATLAS Collaboration], JHEP 1501, 069 (2015) [arXiv:1409.6212 [hep-ex]].
  * (111) The ATLAS collaboration, ATLAS-CONF-2014-061, ATLAS-COM-CONF-2014-080.
  * (112) V. Khachatryan et al. [CMS Collaboration], Eur. Phys. J. C 75, no. 5, 212 (2015) [arXiv:1412.8662 [hep-ex]].
  * (113) S. Heinemeyer et al. [LHC Higgs Cross Section Working Group Collaboration], arXiv:1307.1347 [hep-ph].
  * (114) J. R. Ellis, M. K. Gaillard and D. V. Nanopoulos, Nucl. Phys. B 106, 292 (1976).
  * (115) M. A. Shifman, A. I. Vainshtein, M. B. Voloshin and V. I. Zakharov, Sov. J. Nucl. Phys. 30, 711 (1979) [Yad. Fiz. 30, 1368 (1979)].
  * (116) A. Djouadi, Phys. Rept. 457, 1 (2008) [hep-ph/0503172].
  * (117) M. Carena, I. Low and C. E. M. Wagner, JHEP 1208, 060 (2012) [arXiv:1206.1082 [hep-ph]].
  * (118) C. S. Chen, C. Q. Geng, D. Huang and L. H. Tsai, Phys. Rev. D 87, 075019 (2013) [arXiv:1301.4694 [hep-ph]].
  * (119) K. Griest and D. Seckel, Phys. Rev. D 43, 3191 (1991).
  * (120) P. A. R. Ade et al. [Planck Collaboration], Astron. Astrophys. (2014) [arXiv:1303.5076 [astro-ph.CO]].
  * (121) P. Gondolo and G. Gelmini, Nucl. Phys. B 360, 145 (1991).
  * (122) J. Edsjo and P. Gondolo, Phys. Rev. D 56, 1879 (1997) [hep-ph/9704361].
  * (123) J. Hisano, K. Ishiwata and N. Nagata, Phys. Rev. D 82, 115007 (2010) [arXiv:1007.2601 [hep-ph]].
  * (124) T. Abe and R. Sato, JHEP 1503, 109 (2015) [arXiv:1501.04161 [hep-ph]].
  * (125) G. Belanger, F. Boudjema, A. Pukhov and A. Semenov, Comput. Phys. Commun. 185, 960 (2014) [arXiv:1305.0237 [hep-ph]].
  * (126) D. S. Akerib et al. [LUX Collaboration], Phys. Rev. Lett. 112, 091303 (2014) [arXiv:1310.8214 [astro-ph.CO]].
  * (127) S. Baek, P. Ko and W. I. Park, JHEP 1202, 047 (2012) [arXiv:1112.1847 [hep-ph]].
  * (128) S. Baek, P. Ko, W. I. Park and E. Senaha, JHEP 1305, 036 (2013) [arXiv:1212.2131 [hep-ph]].
  * (129) J. M. Cline and K. Kainulainen, JCAP 1301, 012 (2013) [arXiv:1210.4196 [hep-ph]].
  * (130) J. M. Cline, K. Kainulainen, P. Scott and C. Weniger, Phys. Rev. D 88, 055025 (2013) [arXiv:1306.4710 [hep-ph]].
  * (131) K. Cheung and O. Seto, Phys. Rev. D 69, 113009 (2004) [hep-ph/0403003].
  * (132) The ATLAS collaboration [ATLAS Collaboration], ATLAS-CONF-2015-004, ATLAS-COM-CONF-2015-004.