#  Heavy singly charged Higgs bosons and inverse seesaw neutrinos as origins
of large $(g-2)_{e,\mu}$ in two Higgs doublet models

L. T. Hue<EMAIL_ADDRESS>Subatomic Physics Research Group, Science and
Technology Advanced Institute, Van Lang University, Ho Chi Minh City 70000,
Vietnam  Faculty of Applied Technology, School of Engineering and Technology,
Van Lang University, Ho Chi Minh City 70000, Vietnam  A. E. Cárcamo Hernández
<EMAIL_ADDRESS>Universidad Técnica Federico Santa María, Casilla
110-V, Valparaíso, Chile  Centro Científico-Tecnológico de Valparaíso, Casilla
110-V, Valparaíso, Chile  Millennium Institute for Subatomic Physics at High-
Energy Frontier (SAPHIR), Fernández Concha 700, Santiago, Chile  H. N. Long
<EMAIL_ADDRESS>Institute of Physics, Vietnam Academy of Science and
Technology, 10 Dao Tan, Ba Dinh, Hanoi 10000, Vietnam  T. T. Hong
111Corresponding author<EMAIL_ADDRESS>An Giang University, VNU - HCM, Ung
Van Khiem Street, Long Xuyen, An Giang 88000, Vietnam

###### Abstract

We show that simple extensions of two Higgs doublet models consisting of new
heavy neutrinos and a singly charged Higgs boson singlet can successfully
explain the experimental data on muon and electron anomalous magnetic moments
thanks to large chirally-enhanced one-loop level contributions. These
contributions arise from the large couplings of inverse seesaw neutrinos with
singly charged Higgs bosons and right-handed charged leptons. The regions of
parameter space satisfying the experimental data on $(g-2)_{e,\mu}$ anomalies
allow heavy singly charged Higgs boson masses above the TeV scale, provided
that heavy neutrino masses are above few hundred GeV, the non-unitary part of
the active neutrino mixing matrix must be large enough, two singly charged
Higgs bosons are non degenerate, and the mixing between singly charged Higgs
bosons must be non-zero.

##  I Introduction

The latest experimental measurement for the anomalous magnetic moment (AMM) of
the muon $a_{\mu}\equiv(g-2)_{\mu}/2$ has been reported from Fermilab
Abi:2021gix  and is in agreement with the previous experimental result
measured by Brookhaven National Laboratory (BNL) E82 Muong-2:2006rrc . A
combination of these results in the new average value of
$a^{\mathrm{exp}}_{\mu}=116592061(41)\times 10^{-11}$, which leads to the
improved standard deviation of 4.2 $\sigma$ from the Standard Model (SM)
prediction, namely

$\Delta a^{\mathrm{NP}}_{\mu}\equiv
a^{\mathrm{exp}}_{\mu}-a^{\mathrm{SM}}_{\mu}=\left(2.51\pm 0.59\right)\times
10^{-9},$ (1)

where $a^{\mathrm{SM}}_{\mu}=116591810(43)\times 10^{-11}$ is the SM
prediction Aoyama:2020ynm  combined from various different contributions
Davier:2010nc ; Davier:2017zfy ; Keshavarzi:2018mgv ; Colangelo:2018mtw ;
Hoferichter:2019mqg ; Davier:2019can ; Keshavarzi:2019abf ; Kurz:2014wya ;
Melnikov:2003xd ; Masjuan:2017tvw ; Colangelo:2017fiz ; Hoferichter:2018kwz ;
Gerardin:2019vio ; Bijnens:2019ghy ; Colangelo:2019uex ; Colangelo:2014qya ;
Blum:2019ugy ; Aoyama:2012wk ; Aoyama:2019ryr ; Czarnecki:2002nt ;
Gnendiger:2013pva . On the other hand, the recent experimental $a_{e}$ data
was reported from different groups Hanneke:2008tm ; Parker:2018vye ;
Morel:2020dww , leading to the latest deviation between experiment and the SM
prediction Aoyama:2012wj ; Aoyama:2012wk ; Laporta:2017okg ; Aoyama:2017uqe ;
Terazawa:2018pdc ; Volkov:2019phy ; Gerardin:2020gpp  as follows:

$\Delta a^{\mathrm{NP}}_{e}\equiv
a^{\mathrm{exp}}_{e}-a^{\mathrm{SM}}_{e}=\left(4.8\pm 3.0\right)\times
10^{-13}.$ (2)

Many Beyond Standard Model (BSM) theories have been constructed to explain the
$(g-2)_{\mu,e}$ anomalies. Such theories rely on the inclusion of vector-like
lepton multiplets Dermisek:2013gta ; Crivellin:2018qmi ; Escribano:2021css ;
Hernandez:2021tii ; Crivellin:2021rbq ; Dermisek:2021ajd ; Chun:2020uzw ;
Frank:2020smf ; Endo:2020tkb ; Chen:2020tfr ; Hati:2020fzp ; Ferreira:2021gke
; Borah:2021khc ; Bharadwaj:2021tgp ; DeJesus:2020yqx ; Cogollo:2020nrc ;
Saez:2021qta , leptoquarks Bigaran:2020jil ; Crivellin:2020tsz ; Zhang:2021dgl
; Keung:2021rps , both neutral and charged Higss bosons as $SU(2)_{L}$
singlets Mondal:2021vou ; Kang:2021jmi . Some two Higgs doublet models (THDM)
can provide sizeable two-loop contributions to $\Delta a_{\mu}$ arising from
new $SU(2)_{L}$ Higgs doublets CarcamoHernandez:2021iat ; Li:2020dbg ;
DelleRose:2020oaa ; Botella:2020xzf ; Han:2021gfu ; Chen:2021jok ; Han:2018znu
; Rose:2021cav , where some of them require rather light masses of new neutral
and/or charged Higgs bosons at few hundred GeV.

There are many types of different versions of THDM, such as for instance: type
I, type II, type X and type Y, which are discussed in Ref. Aoki:2009ha , where
collider phenomenology will give rise to different physical results. They are
distinguished among themselves by the ways in which the two different
$SU(2)_{L}$ Higgs doublets $\Phi_{1,2}$ generate the masses of up and down
quarks as well as of charged leptons. These models arise from suitable choices
of $Z_{2}$ charge assignments. The Yukawa couplings of Higgs bosons with heavy
fermions will be constrained by the perturbative limits, for example
$0.4\leq\tan\beta\leq 91$ corresponding to the definitions of Yukawa couplings
of the top and bottom quarks
$\left|Y_{t}\right|=(\sqrt{2}/v)m_{t}\cot\beta<\sqrt{\pi}$ and
$\left|Y_{b}\right|=(\sqrt{2}/v)m_{b}\tan\beta<\sqrt{\pi}$, where $\tan\beta$
is the ratio between the two vacuum expectations values (vevs) of the two
neutral Higgs components. Discussions on the original THDM have shown that a
new $SU(2)_{L}$ Higgs doublet can be used to accommodate the $(g-2)_{e,\mu}$
anomalies in a rather narrow allowed region of parameter space at the price of
imposing dangerous requirements that can be experimentally tested in the near
future Li:2020dbg ; Jueid:2021avn . The conclusion is also true for the
Minimal Supersymmetric Model (MSSM) Badziak:2019gaf ; Li:2021koa ;
Athron:2021iuf . Another solution for MSSM requires large SUSY threshold
corrections enough to change both sizes and signs of the SUSY electron and
muon Yukawa couplings Endo:2019bcj , hence allow regions satisfying both
$(g-2)_{e,\mu}$ anomalies with large TeV scale of slepton masses, but rather
large $t_{\beta}\geq 70$ is needed. As usual solution, there is a variety of
extensions of THDM that successfully explain and accommodate the
$(g-2)_{e,\mu}$ anomalies by adding new vector-like fermions as $SU(2)_{L}$
multiplets Chun:2020uzw ; Frank:2020smf ; Hernandez:2021tii ; Chen:2020tfr ;
Dermisek:2021ajd ; Ferreira:2021gke  as well as $SU(2)_{L}$ singlets of
neutral and charged Higss bosons CarcamoHernandez:2019xkb ;
CarcamoHernandez:2019cbd ; Mondal:2021vou ; CarcamoHernandez:2021iat ;
Hernandez:2021xet ; Adhikari:2021yvx . Only few THDM models with scalar
singlets and fermions can successfully explain the experimental data of both
muon and electron anomalous magnetic moments Chen:2020tfr ; Dutta:2020scq ;
Hernandez:2021xet . Our work here will introduce a more general solution that
a wide class of the THDM by adding heavy right handed (RH) neutrino singlets
and one singly electrically charged Higgs boson can give sizeable one-loop
contributions enough to explain the experimental ranges of values of both muon
and electron anomalous magnetic moments in a wide region of parameter space.

This paper is organized as follows. In Sec. II, the THDM extension with new RH
neutrinos and singly charged Higgs boson (THDM$N_{R}h^{\pm}$) will be
introduced, where we pay attention to the leptons, gauge bosons, and Higgs
sectors, giving all physical states as well the couplings that may give large
one-loop contributions to AMM. In Sec. III, the simple inverse seesaw (ISS)
and analytic formulas for one-loop contributions to AMM are constructed.
Numerical discussion will be shown in details. Our main results will be
collected in Sec. IV.

##  II  Review of the THDM$N_{R}h^{\pm}$

In this section we will focus on the analysis of the lepton sector of the
THDM$N_{R}h^{\pm}$. For details of the quark sector of different types of the
THDM, the reader is referred to Refs. Davidson:2005cw ; Herrero-Garcia:2017xdu
. The leptonic, quark and scalar spectrum with their assignments under the
$SU(2)_{L}\times U(1)_{Y}$ gauge group are given by Davidson:2005cw ; Herrero-
Garcia:2017xdu :

$\displaystyle L_{a}$ $\displaystyle=\begin{pmatrix}\nu_{aL}\\\
e_{aL}\end{pmatrix}\sim(2,-1);\;\quad Q_{aL}=\begin{pmatrix}u_{aL}\\\
d_{aL}\end{pmatrix}\sim\left(2,\dfrac{1}{3}\right),$ $\displaystyle e_{aR}$
$\displaystyle\sim(1,-2);\quad u_{aR}\sim\left(1,\dfrac{4}{3}\right);\quad
d_{aR}\sim\left(1,-\dfrac{2}{3}\right);\quad
N_{IR}\sim(1,0),\;a=1,2,3,\;I=1,2,...,6;$
$\displaystyle\Phi_{i}=\begin{pmatrix}\phi^{+}_{i}\\\
\phi^{0}_{i}\end{pmatrix}\sim(2,1),\;\langle\Phi_{i}\rangle=\begin{pmatrix}0\\\
\frac{v_{i}}{\sqrt{2}}\end{pmatrix},\;i=1,2;$ $\displaystyle h^{\pm}$
$\displaystyle\sim(1,\pm 2).$ (3)

Note that we use the old convention for the defining the hypercharge via the
original Gell-Mann- Nishijima relation:

$Q=T_{3}+\frac{1}{2}Y.$ (4)

Here we add six new neutrino singlets to the leptonic spectrum of the model in
order to implement the standard ISS mechanism. Furthermore, the model scalar
sector is augmented by the inclusion of singly charged Higgs bosons, which
allow to generate the new Yukawa interaction
$\overline{(N_{IR})^{c}}e_{aR}h^{+}$ thus resulting in sizeable one-loop
contributions to AMM.

One of the most important parameters introduced in the THDM is

$t_{\beta}\equiv\tan\beta=\dfrac{v_{2}}{v_{1}},\;s_{\beta}\equiv\sin\beta,\;c_{\beta}\equiv\cos\beta.$
(5)

We will also follow the notations $s_{x}\equiv\sin x$, $c_{x}\equiv\cos x$,
and $t_{x}=s_{x}/c_{x}$ for any parameter $x$. Because the Yukawa couplings of
up type quarks and $\Phi_{2}$ are fixed, we always have a lower bound
$t_{\beta}\geq 0.4$. In the THDM type I, all fermions couple with only
$\Phi_{2}$, hence no upper bounds for $t_{\beta}$ arise. In the THDM type II
and Y, all down type quarks and charged leptons couple with $\Phi_{1}$, hence
$t_{\beta}\leq 91$. In the type-X, the Yukawa coupling of tau with $\Phi_{1}$
gives $t_{\beta}\leq 348$.

The THDM type-II and type-X where lepton doublets couple with only $\Phi_{1}$
are called THDM type-A. The corresponding lepton Yukawa terms are:

$\displaystyle-\mathcal{L}^{\text{Y}}_{\ell}=$ $\displaystyle
Y^{e}_{ab}\overline{e_{aR}}\Phi_{1}^{\dagger}L_{b}+Y^{N}_{k,Ia}\overline{N_{IR}}\left(i\sigma_{2}\Phi_{k}\right)^{T}L_{a}+\frac{1}{2}M_{N,IJ}\overline{N_{IR}}(N_{JR})^{c}$
$\displaystyle+Y^{h}_{Ia}\overline{(N_{IR})^{c}}e_{aR}h^{+}+\mathrm{h.c.},$
(6)

where $a,b=1,2,3$ are family indices; $k=1,2$; and $I=1,2,3,...6$ are the
number of new neutral lepton singlets. A rough condition for the pertubative
limit of $Y^{h}$ is $|Y^{h}_{Ia}|<\sqrt{4\pi}$, and should be smaller for
trusty values Allwicher:2021rtd . The $Z_{2}$ charges of $N_{IR}$ and
$h^{\pm}$ are always chosen to guarantee the $Z_{2}$ symmetry of Lagrangian
(II). Using the $Z_{2}$ charge assignments discussed in Ref. Aoki:2009ha , two
Higgs triplets $\Phi_{1,2}$ have different $Z_{2}$ charges, therefore $N_{IR}$
couple with only one of them. In addition, we can consider the inclusion of
new $Z_{3}$ discrete symmetry in the model in order to forbid the quartic
scalar interaction
$\left[\left(\Phi_{1}^{\dagger}\Phi_{2}\right)^{2}+\mathrm{H.c.}\right]$ that
generates active neutrino masses at the one-loop level Ma:2006km . This
guarantees that the Yukawa couplings $Y^{h}_{Ia}$ do not have suppressed
constraints arising from the neutrino oscillation experimental data. The
detailed assignments of the $Z_{2}$ and $Z_{3}$ charges are shown in table 1.

Table 1:  $Z_{2}\times Z_{3}$ charge assignments of the model type-A for the case I (II) with heavy neutrinos couplings with $\Phi_{1}$ ($\Phi_{2}$), $w=e^{i2\pi/3}$ | $Q_{aL}$ | $u_{aR}$ | $d_{aR}$ | $L_{a}$ | $e_{aR}$ | $N_{1,2,3R}$ | $N_{4,5,6R}$ | $\Phi_{1}$ | $\Phi_{2}$ | $h^{+}$
---|---|---|---|---|---|---|---|---|---|---
$Z_{2}$ | - | - | - | 1 | 1 | 1(-1) | -1(1) | 1 | -1 | -1(1)
$Z_{3}$ | $w$ | $w^{2}$ | $w^{2}$ | $w$ | $w^{2}$ | 1 | $w$ | $w^{2}$ | $w^{2}$ | $w$

The $Z_{3}$ charges of leptons are chosen for the implementation of the ISS
mechanism in the total neutrino mass matrix, namely: i) the charged lepton
mass term always consists of $\Phi_{1}$, ii) only $N_{aR}$ ($a=1,2,3$)
generate non-zero Dirac mass term of neutrinos, iii) and only $N_{(a+3)R}$
couple with $h^{\pm}$. The $Z_{2}$ charges of quarks depend on which types of
well-known THDMs, hence we do not list here. The $Z_{3}$ charges of all quarks
have the same values as those of lepton doublets and charged lepton singlets.
This guarantees that the quark Yukawa terms are consistent with those
previously discussed.

Different from the well-known THDM, the model under considerations have new
kind of Yukawa couplings with singly charged Higgs boson $h^{\pm}$ and right-
handed neutrinos $N_{I}$. In the SM, the vacuum stability around the Plank
scale requires the quartic coupling in the Higgs potential close to zero, and
the upper bound of the SM Higgs boson mass is very close to the experimental
value, $M_{h}<126$ GeV Degrassi:2012ry . This problem can be solved in the
model by adding new scalars such as in the THDM Jangid:2020qgo , and even with
the models by adding neutrino gauge singlets Jangid:2020dqh ; Coriano:2015sea
; DelleRose:2015bms . The reason is that additional scalar couplings to Higgs
doublets generating SM-like Higgs boson mass give positive one-loop
contributions to the $\beta$-functions of the quartic couplings, hence the
negative values of these couplings implying the vacuum instability will be
relaxed at higher energy scale. This consequence still hold with the
appearance of large Yukawa couplings of new fermions with Higgs doublets,
which give negative one-loop contributions to the $\beta$-functions of the
quartic couplings. On the other hand, the upper bounds of these Yukawa
couplings are required besides their bounds from the perturbative and
unitarity limits, previously discussed for THDM in refs. Bandyopadhyay:2020djh
, for example. This situation also applies to the THDM$N_{R}h^{\pm}$ model
under consideration. On the other hand, the appearance of the singly charged
Higgs boson $h^{\pm}$ resulting on the Yukawa coupling matrix $Y^{h}$, which
is irrelevant with any Higgs doublets. Hence $Y^{h}$ is not constrained by the
requirement of vacuum stability. Furthermore, the quartic couplings of
$h^{\pm}$ with Higgs doublets in the Higgs potential will result in the higher
energy scale satisfying the vacuum stability.

In the THDM type-I and type-Y, where charged lepton couple with $\Phi_{2}$,
the first term of Eq. (II) should replace $\Phi_{1}$ with $\Phi_{2}$. We call
them the THDM type-B and will discuss later that the qualitative results for
AMM do not change. This model type including the model introduced in Ref.
Chen:2020tfr .

We also emphasize that the models under consideration including the Yukawa
part discussed on Ref. Mondal:2021vou  focusing on one particular case of
$\Phi_{k}=\Phi_{2}$. As we will show that in our numerical analysis, the
regions of the parameter space predicting $\Delta a_{e,\mu}$ consistent with
the experimental data, favor small values of $t_{\beta}$, including the range
$t_{\beta}\in[0.4,\;10]$ which was not mentioned in Ref. Mondal:2021vou .

Defining $\nu^{\prime}_{L}=(\nu_{1},\nu_{2},\nu_{3})_{L}^{T}$ and
$N_{R}=(N_{1},N_{2},...,N_{6})_{R}^{T}$, we find the following leptonic mass
terms:

$\displaystyle-\mathcal{L}^{\text{mass}}_{\mathrm{lepton}}$
$\displaystyle=\frac{Y^{e}_{ab}v_{1}}{\sqrt{2}}\overline{e^{\prime}_{aR}}e^{\prime}_{bL}+\frac{1}{2}\begin{pmatrix}\overline{(\nu^{\prime}_{L})^{c}}&\overline{N_{R}}\end{pmatrix}\mathcal{M}^{\nu}\begin{pmatrix}\nu^{\prime}_{L}\\\
(N_{R})^{c}\end{pmatrix}+\mathrm{h.c.},\quad\mathcal{M}^{\nu}=\begin{pmatrix}0_{3\times
3}&M^{T}_{D}\\\ M_{D}&M_{N}\end{pmatrix}\;,$ (7)

where

$(M_{D})_{Ia}\equiv M_{D,Ia}=\sum_{k=1}^{2}\frac{Y^{N}_{k,Ia}v_{k}}{\sqrt{2}}$
(8)

are the components of a $6\times 3$ Dirac neutrino mass matrix and $M_{N}$ is
a $6\times 6$ symmetric Majorana mass matrix. Here we choose basis where the
charged lepton mass matrix is diagonal, which implies
$Y^{e}_{ab}=\delta_{ab}\sqrt{2}m_{e_{a}}/v_{1}$, i.e., the flavor and mass
states of charged leptons are the same. The total unitary mixing matrix is
defined as

$\displaystyle U^{\nu T}\mathcal{M}^{\nu}U^{\nu}$
$\displaystyle=\hat{\mathcal{M}}^{\nu}=\mathrm{diag}(m_{n_{1}},\;m_{n_{2}},\;...,m_{n_{9}})\equiv\mathrm{diag}(\hat{m}_{\nu},\;\hat{M}_{N}),$
$\displaystyle\begin{pmatrix}\nu^{\prime}_{L}\\\ (N_{R})^{c}\end{pmatrix}$
$\displaystyle=U^{\nu}n_{L},\;\begin{pmatrix}(\nu^{\prime}_{L})^{c}\\\
N_{R}\end{pmatrix}=U^{\nu*}n_{R}=U^{\nu*}(n_{L})^{c},\;$ (9)

where $n_{L,R}=(n_{1},n_{2},...,n_{9})_{L,R}$ are the Majorana neutrino mass
eigenstates satisfying $n_{iL,R}=(n_{iR,L})^{c}$, and the four-component forms
are $n_{i}=(n_{iL},n_{iR})^{T}$.

Given that we are working in the basis where the charged lepton mass matrix is
diagonal, the leptonic mixing entirely arises from the neutrino sector. This
implies that, in order to successfully reproduce the neutrino oscillation
experimental data, the neutrino mixing matrix is parameterised in the
following form Ibarra:2010xw :

$U^{\nu}=\left(\begin{array}[]{cc}\left(I_{3}-\frac{1}{2}RR^{\dagger}\right)U_{\mathrm{PMNS}}&RV\\\
-R^{\dagger}U_{\mathrm{PMNS}}&\left(I_{6}-\frac{1}{2}R^{\dagger}R\right)V\\\
\end{array}\right)+\mathcal{O}(R^{3}),$ (10)

where $V$ is a $6\times 6$ unitary matrix; $R$ is a $3\times 6$ matrix
satisfying $|R_{aI}|<1$ for all $a=1,2,3$, and $I=1,2,...,6$. The $3\times 3$
unitary matrix $U_{\mathrm{PMNS}}$ is the Pontecorvo-Maki-Nakagawa-Sakata
(PMNS) matrix ParticleDataGroup:2020ssz . The condition $|R_{aI}|<1$ needed to
guarantee that the $U^{\nu}$ as a power series of $R$ given in Ref.
Ibarra:2010xw  is convergent and the approximation up to the order $R^{2}$
results in seesaw relations well-known in literature. Namely, the total
neutrino mass matrix given in Eq. (7) will result in three light masses for
active neutrinos and new heavy neutrinos, $m_{n_{1,2,3}}\ll
m_{n_{4,5,\dots}}$. This form of $U^{\nu}$ is consistent with many other well-
known parameterizations Schechter:1981cv ; Korner:1992zk ; Grimus:2000vj .

The successful implementation of the ISS mechanism requires to construct the
Dirac and Majorana mass matrices in terms of $3\times 3$ submatrices as
follows Arganda:2014dta ; Ibarra:2010xw

$\displaystyle M^{T}_{D}=(m^{T}_{D},\;0_{3\times 3}),\hskip
14.22636ptM_{N}=\left(\begin{array}[]{cc}0_{3\times 3}&M_{R}\\\
M^{T}_{R}&\mu_{X}\\\ \end{array}\right),$ (13)

where $0_{3\times 3}$ is the $3\times 3$ null matrix. Using a new notation
$M=M_{R}\mu_{X}^{-1}M_{R}^{T}$, we have the following ISS relations:

$\displaystyle R$
$\displaystyle=M^{\dagger}_{D}{M^{*}_{N}}^{-1}=\left(-m_{D}^{\dagger}M^{*-1},\hskip
14.22636ptm^{\dagger}_{D}\left(M^{\dagger}_{R}\right)^{-1}\right),$
$\displaystyle U^{*}_{\mathrm{PMNS}}\hat{m}_{\nu}U^{\dagger}_{\mathrm{PMNS}}$
$\displaystyle=m_{\nu}=-M^{T}_{D}M_{N}^{-1}M_{D}=m_{D}^{T}\left(M_{R}^{T}\right)^{-1}\mu_{X}M_{R}^{-1}m_{D},$
$\displaystyle V^{*}\hat{M}_{N}V^{\dagger}$ $\displaystyle\simeq
M_{N}+\frac{1}{2}R^{T}R^{*}M_{N}+\frac{1}{2}M_{N}R^{\dagger}R.$ (14)

Because $m_{D}$ is parameterized in terms of many free parameters, hence it is
enough to choose $\mu_{X}=\mu_{0}I_{3}$ with $\mu_{0}>0$,
$M_{R}=\hat{M}_{R}=M_{0}I_{3}$. We choose a simple form
$m_{D}=\frac{M_{0}}{\sqrt{\mu_{0}}}\sqrt{\hat{m}_{\nu}}U^{\dagger}_{\mathrm{PMNS}}$
Casas:2001sr ; Arganda:2014dta ; Ibarra:2010xw . The ISS condition
$|\hat{m}_{\nu}|\ll|\mu_{X}|\ll|m_{D}|\ll M_{0}$ gives
$\frac{\sqrt{\mu_{0}\hat{m}_{\nu}}}{M_{0}}\simeq 0$ and

$\displaystyle\hat{M}_{N}=\left(\begin{matrix}\hat{M}_{R}&0_{3\times 3}\\\
0_{3\times 3}&\hat{M}_{R}\end{matrix}\right)\simeq
M_{0}I_{6},\;V\simeq\dfrac{1}{\sqrt{2}}\left(\begin{matrix}-iI_{3}&I_{3}\\\
iI_{3}&I_{3}\end{matrix}\right),$ (15)

meaning that all of six heavy neutrinos have approximately the same mass
$m_{n_{i}}\simeq M_{0}$ for all $i>3$. Hence, the above simple forms of
$M_{R}$, and $\mu_{X}$ will result in degenerate heavy neutrino masses,
leading to small rates of lepton flavor violating decays of charged leptons
(cLFV) that satisfy the current experimental constraints MEG:2016leq ;
BaBar:2009hkt . Relations in (14) reduce to the following simple forms:

$\displaystyle m_{D}$
$\displaystyle=M_{0}\hat{x}_{\nu}^{1/2}U^{\dagger}_{\mathrm{PMNS}},\quad
R=\left(-U_{\mathrm{PMNS}}\frac{\sqrt{\mu_{0}\hat{m}_{\nu}}}{M_{0}},\;U_{\mathrm{PMNS}}\hat{x}_{\nu}^{1/2}\right)\simeq\left(0_{3\times
3},\;U_{\mathrm{PMNS}}\hat{x}_{\nu}^{1/2}\right),$ (16)

where

$\hat{x}_{\nu}\equiv\frac{\hat{m}_{\nu}}{\mu_{0}}\equiv\mathrm{diag}\left(\hat{x}_{\nu
1},\hat{x}_{\nu 2},\hat{x}_{\nu 3}\right),\;\hat{x}_{\nu
a}\equiv\frac{m_{n_{a}}}{\mu_{0}},\;a=1,2,3$ (17)

satisfying max$[|\hat{x}_{\nu a}|]\ll 1$ for all $a=1,2,3$.

In our numerical analysis, we will use the best-fit values of the neutrino
oscillation data ParticleDataGroup:2020ssz  corresponding to the normal order
(NO) scheme with $m_{n_{1}}<m_{n_{2}}<m_{n_{3}}$, namely

$\displaystyle
s^{2}_{12}=0.32,\;s^{2}_{23}=0.547,\;s^{2}_{13}=0.0216,\;\delta=218\;[\mathrm{Deg}],$
$\displaystyle\Delta m^{2}_{21}=7.55\times 10^{-5}[\mathrm{eV}^{2}],\;\Delta
m^{2}_{32}=2.424\times 10^{-3}[\mathrm{eV}^{2}].$ (18)

In our numerical analysis, we have used the following relations

$\displaystyle\hat{x}_{\nu}$
$\displaystyle=\mu_{0}^{-1}\mathrm{diag}\left(m_{n_{1}},\;\sqrt{m_{n_{1}}^{2}+\Delta
m^{2}_{21}},\;\sqrt{m_{n_{1}}^{2}+\Delta m^{2}_{21}+\Delta
m^{2}_{32}}\right),$ $\displaystyle U_{\mathrm{PMNS}}$
$\displaystyle=\left(\begin{array}[]{ccc}c_{12}c_{13}&c_{13}s_{12}&s_{13}e^{-i\delta}\\\
-c_{23}s_{12}-c_{12}s_{13}s_{23}e^{i\delta}&c_{12}c_{23}-s_{12}s_{13}s_{23}e^{i\delta}&c_{13}s_{23}\\\
s_{12}s_{23}-c_{12}c_{23}s_{13}e^{i\delta}&-c_{23}s_{12}s_{13}e^{i\delta}-c_{12}s_{23}&c_{13}c_{23}\\\
\end{array}\right)$ (22)
$\displaystyle\simeq\left(\begin{array}[]{ccc}0.816&0.56&0.147e^{-i\delta}\\\
-0.381-0.09e^{i\delta}&0.555-0.062e^{i\delta}&0.732\\\
0.418-0.082e^{i\delta}&-0.61-0.056e^{i\delta}&0.666\\\ \end{array}\right).$
(26)

These numerical values of neutrino masses satisfy the cosmological constraint
arising from the Planck 2018 experimental data Planck:2018vyg :
$\sum_{i=a}^{3}m_{n_{a}}\leq 0.12\;\mathrm{eV}$. In order to simplify our
numerical analysis, we assume
$m_{n_{1}}=0.01\;\mathrm{eV}<m_{n_{2}}<m_{n_{3}}$.

The other well-known numerical parameters are given in Ref.
ParticleDataGroup:2020ssz , namely

$\displaystyle g$
$\displaystyle=0.652,\;\alpha_{e}=\frac{1}{137}=\frac{e^{2}}{4\pi},\;s^{2}_{W}=0.231,$
$\displaystyle m_{e}$ $\displaystyle=5\times
10^{-4}\;\mathrm{GeV},\;m_{\mu}=0.105\;\mathrm{GeV},\;m_{W}=80.385\;\mathrm{GeV}.$
(27)

Also the inverted order (IO) scheme with $m_{n_{3}}<m_{n_{1}}<m_{n_{2}}$ can
be considered in the similar way, but the qualitative results are the same
with those from NO scheme, so we will not present here.

The non-unitary of the active neutrino mixing matrix
$\left(I_{3}-\frac{1}{2}RR^{\dagger}\right)U_{\mathrm{PMNS}}$ is constrained
by other phenomenological aspects such as, for instance, electroweak precision
tests, cLFV decays Fernandez-Martinez:2016lgt ; Pinheiro:2021mps ;
Agostinho:2017wfs , thus leading to the following constraints

$\displaystyle\eta\equiv\frac{1}{2}\left|RR^{\dagger}\right|<\eta_{0}=\begin{pmatrix}2\times
10^{-3}&3.5\times 10^{-5}&8.\times 10^{-3}\\\ 3.5\times 10^{-5}&8\times
10^{-4}&5.1\times 10^{-3}\\\ 8\times 10^{-3}&5.1\times 10^{-3}&2.7\times
10^{-3}\end{pmatrix}.$ (28)

This constraint is consistent with the data popularly used in recent works
discussion on the ISS framework Bandyopadhyay:2012px ; Dao:2021vqp ;
Mondal:2021vou . The constraint on $\eta$ may be more strict, depending on
particular models. For example in the type III and inverse seesaw models, one
has the constraint $|\eta_{aa}|\leq\mathcal{O}(10^{-4})$ Biggio:2019eeo ;
Escribano:2021css . In our numerical analysis, we will choose the values
satisfying $|\eta_{33}|\leq 10^{-3}$, which are also consistent with the
updated constraints on neutrino couplings discussed in Ref. Coutinho:2019aiy ;
Manzari:2020eum , including the constraint from lepton universality previously
discussed deBlas:2013gla .

In the next section, we will discuss the model scalar potential as well as the
one-loop contributions to $\Delta a_{\mu,e}$.

##  III  Higgs bosons and one-loop contributions to $\Delta a_{e_{a}}$

The Higgs potential satisfying the symmetry mentioned in this work is

$\displaystyle V$
$\displaystyle=\mu_{1}^{2}\Phi_{1}^{\dagger}\Phi_{1}+\mu_{2}^{2}\Phi_{2}^{\dagger}\Phi_{2}-\left(\mu_{3}^{2}\Phi^{\dagger}_{2}\Phi_{1}+\mathrm{H.c.}\right)$
$\displaystyle+\dfrac{1}{2}\lambda_{1}\left(\Phi_{1}^{{\dagger}}\Phi_{1}\right)^{2}+\dfrac{1}{2}\lambda_{2}\left(\Phi_{2}^{{\dagger}}\Phi_{2}\right)^{2}+\lambda_{3}\left(\Phi_{1}^{{\dagger}}\Phi_{1}\right)\left(\Phi_{2}^{{\dagger}}\Phi_{2}\right)+\lambda_{4}\left(\Phi_{1}^{{\dagger}}\Phi_{2}\right)\left(\Phi_{2}^{{\dagger}}\Phi_{1}\right)$
$\displaystyle+\lambda_{h}\left|{h^{+}}\right|^{4}+\left|{h^{+}}\right|^{2}\left[\mu_{h}^{2}+\lambda_{8}\Phi_{1}^{{\dagger}}\Phi_{1}+\lambda_{9}\Phi_{2}^{{\dagger}}\Phi_{2}\right]+\left(\mu\epsilon_{\alpha\beta}\Phi_{1}^{\alpha}\Phi_{2}^{\beta}h^{-}+\mathrm{H.c.}\right),$

where the $Z_{2}$ soft breaking term
$\left(\mu_{3}^{2}\Phi^{\dagger}_{2}\Phi_{1}\mathrm{+H.c.}\right)$ is kept in
order to generate non-zero mass for the CP-odd neutral Higgs predicted in this
model. The quartic term $\left(\Phi^{\dagger}_{2}\Phi_{1}\right)^{2}$ vanishes
because of the $Z_{3}$ symmetry, whose assignments for the particle spectrum
are given in Table 1. Consequently, the one-loop contributions to active
neutrino masses discussed in Ref. Ma:2006km  do not appear in this case,
implying that all entries of $Y^{h}$ are not constrained by this condition. In
general, $\mu_{3}$, and $\mu$ can be complex, while the remaining parameters
$\mu_{1,2}$, $\mu_{h}$, and $\lambda_{1,2,3,4,8,9,h}$ are real Herrero-
Garcia:2017xdu ; Davidson:2005cw . In this work, all of these parameters and
vacuum expectation values (vev) are assumed to be real, which corresponds to
CP conservation.

It is easy to find the two minimization conditions of the Higgs potential,
then inserting them into the Higgs potential in order to find the physical
electrically charged scalars as follows:

$\displaystyle\begin{pmatrix}G^{\pm}\\\ H^{\pm}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\ c_{\beta}&s_{\beta}\\\
-s_{\beta}&c_{\beta}\end{pmatrix}\begin{pmatrix}\phi_{1}^{\pm}\\\
\phi_{2}^{\pm}\end{pmatrix},$ $\displaystyle\begin{pmatrix}h_{1}^{\pm}\\\
h_{2}^{\pm}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}s_{\varphi}&c_{\varphi}\\\
c_{\varphi}&-s_{\varphi}\end{pmatrix}\begin{pmatrix}h^{\pm}\\\
H^{\pm}\end{pmatrix},\quad
s_{2\varphi}=\dfrac{\sqrt{2}v\mu}{m^{2}_{h_{2}^{+}}-m^{2}_{h_{1}^{+}}},$
$\displaystyle m^{2}_{{h_{1}^{+}},{h_{2}^{+}}}$
$\displaystyle\equiv\dfrac{1}{2}\left[M^{2}_{H^{+}}+M^{2}_{33}\mp\sqrt{\left(M^{2}_{H^{+}}-M^{2}_{33}\right)^{2}+2v^{2}\mu^{2}}\right],$
(29)

where $t_{\beta}\equiv\tan\beta=\dfrac{v_{2}}{v_{1}}$, $s_{x}\equiv\sin x$,
$c_{x}\equiv\cos x$,
$M_{H^{+}}^{2}=\frac{\mu_{3}^{2}}{s_{\beta}c_{\beta}}-\dfrac{1}{2}v^{2}\lambda_{4}$,
$M_{33}^{2}=\mu_{h}^{2}+\frac{v^{2}}{2}\left(c_{\beta}^{2}\lambda_{8}+s_{\beta}^{2}\lambda_{9}\right)$.
This result is consistent with the one obtained from the Higgs potential used
in Ref Herrero-Garcia:2017xdu  after some transformations into a new Higgs
basis and parameters. Namely, the two Higgs doublets $\Phi_{i}$ are changed
into a new basis as follows:

$\begin{pmatrix}H_{1}\\\ H_{2}\end{pmatrix}=\begin{pmatrix}\
c_{\beta}&s_{\beta}\\\
-s_{\beta}&c_{\beta}\end{pmatrix}\begin{pmatrix}\Phi_{1}\\\
\Phi_{2}\end{pmatrix},$ (30)

which can be expanded around the minimum as shown below:

$\displaystyle H_{1}$ $\displaystyle=\begin{pmatrix}G^{+}\\\
\dfrac{v+\varphi^{0}_{1}+iG^{0}}{\sqrt{2}}\end{pmatrix},\;H_{2}=\begin{pmatrix}H^{+}\\\
\dfrac{\varphi^{0}_{2}+iA}{\sqrt{2}}\end{pmatrix},$ (31)

where $v_{1}^{2}+v_{2}^{2}=v^{2}=(246\;\mathrm{GeV})^{2}$. In addition, the
free parameters of the Higgs potential are transformed as follows
$\mu_{1}^{2}c_{\beta}^{2}+\mu_{2}^{2}s_{\beta}^{2}-\mu_{3}^{2}2s_{\beta}c_{\beta}\to\mu_{1}^{2}$,
$\mu_{1}^{2}s_{\beta}^{2}+\mu_{2}^{2}c_{\beta}^{2}+\mu_{3}^{2}2s_{\beta}c_{\beta}\to\mu_{2}^{2}$,
$(\mu_{1}^{2}-\mu_{2}^{2})s_{\beta}c_{\beta}+\mu_{3}^{2}(c^{2}_{\beta}-s^{2}_{\beta})\to\mu_{3}^{2}$,
$\lambda_{8}c^{2}_{\beta}+\lambda_{9}s^{2}_{\beta}\to\lambda_{8}$,
$\lambda_{8}s^{2}_{\beta}+\lambda_{9}c^{2}_{\beta}\to\lambda_{9}$,
$s_{\beta}c_{\beta}(\lambda_{9}-\lambda_{8})\to\lambda_{10}$,
$\mu\to\mu,\;\dots$ The Higgs potential after the transformation (30) takes
the form Herrero-Garcia:2017xdu , where the physical charged scalar states and
their masses are consistent with those ones given in Eq. (III).

In our numerical calculation, we will use $m^{2}_{{h_{k}^{+}}}$ and the mixing
angle $\varphi$ as free parameters. Three Higgs mass parameters $\mu_{2}^{2}$,
$\mu_{h}^{2}$, and $\mu$ are functions of the remaining parameters. Thus, no
perturbative limits on the Higgs selfcouplings are necessary to constrain the
dependent functions chosen here.

From the above information we obtain all vertices providing one-loop
contributions to the $e_{b}\to e_{a}\gamma$ decay rates as well as to $\Delta
a_{e_{a}}$. They are collected from the lepton Yukawa terms given in Eq. (II)
respecting all symmetries given in Table (1), namely

$\displaystyle\mathcal{L}^{\text{yuk}}_{\mathrm{lepton}}$
$\displaystyle=-\frac{\sqrt{2}m_{e_{a}}}{v_{1}}\overline{L_{a}}\Phi^{\dagger}_{1}e_{aR}-Y^{N}_{k,ba}\overline{N_{bR}}\left(i\sigma_{2}\Phi_{k}\right)^{T}L_{a}-Y^{h}_{(b+3)a}\overline{(N_{(b+3)R})^{c}}e_{aR}h^{+}+\mathrm{h.c.},$

In order to derive the total neutrino mass matrix in the general ISS form, the
leptonic Yukawa terms are rewritten as follows

$\displaystyle\mathcal{L}^{\text{yuk}}_{\mathrm{lepton}}$
$\displaystyle=-\frac{gm_{e_{a}}(-t_{\beta})}{\sqrt{2}m_{W}}\overline{L_{a}}H_{2}e_{aR}-Y^{N}_{0,Ia}\overline{N_{IR}}\left(i\sigma_{2}H_{2}\right)^{T}L_{a}-Y^{h}_{Ia}\overline{(N_{IR})^{c}}e_{aR}h^{+}+\mathrm{h.c.},$
(32)

where we denote $Y^{N}_{0}\equiv\left(Y^{N}_{0},\;0_{3\times
3}\right)^{T}=-s_{\beta}Y^{N}_{1}+c_{\beta}Y^{N}_{2}$ for the ISS mechanism
discussed in this work. It is interesting to link this matrix with $M_{D}$
given in Eq. (8), where depending on the $Z_{2}$ charges of $N_{IR}$ there are
two cases where non-zero couplings with only $\Phi_{1}$ or $\Phi_{2}$
correspond to $Y^{N}_{2}=0$ or $Y^{N}_{1}=0$, respectively. Namely

$\displaystyle
M_{D}^{T}=\left(\frac{vf_{H}^{-1}}{\sqrt{2}}Y^{N}_{0},\;0_{3\times
3}\right)^{T},\;f_{H}=\left[\begin{array}[]{cc}t_{\beta}^{-1},&Y^{N}_{1}=0_{6\times
3},\;Y^{N}_{2}\neq 0;\\\ -t_{\beta},&Y^{N}_{2}=0_{6\times 3},\;Y^{N}_{1}\neq
0\end{array}\right..$ (35)

With this new notation, the Yukawa Lagrangian for the THDM type-A is

$\displaystyle\mathcal{L}^{\text{yuk}}_{\mathrm{lepton}}=$
$\displaystyle\frac{gm_{e_{a}}t_{\beta}}{\sqrt{2}m_{W}}U^{\nu*}_{ai}\overline{n_{iL}}e_{aR}\left(c_{\varphi}h^{+}_{1}-s_{\varphi}h^{+}_{2}\right)+\frac{gf_{H}}{\sqrt{2}m_{W}}U^{\nu}_{(I+3)i}M_{D,ia}\overline{n_{iR}}e_{aL}\left(c_{\varphi}h^{+}_{1}-s_{\varphi}h^{+}_{2}\right)$
$\displaystyle-Y^{h}_{Ib}U^{\nu*}_{(I+3)i}\overline{n_{iL}}e_{aR}\left(s_{\varphi}h^{+}_{1}+c_{\varphi}h^{+}_{2}\right)+\mathrm{h.c.}.$
(36)

Then, all relevant couplings are given in the following lepton Yukawa
interaction Lagrangian

$\displaystyle\mathcal{L}=\sum_{a=1}^{3}\sum_{i=1}^{9}\left[\frac{g}{\sqrt{2}m_{W}}\sum_{k=1}^{2}\overline{n_{i}}\left(\lambda^{L,k}_{ia}P_{L}+\lambda^{R,k}_{ia}P_{R}\right)e_{a}h^{+}_{k}+\frac{g}{\sqrt{2}}U^{\nu*}_{ai}\overline{n_{i}}\gamma^{\mu}P_{L}e_{a}W^{+}_{\mu}\right]+\mathrm{h.c.},$
(37)

where

$\displaystyle\lambda^{L,1}_{ia}$
$\displaystyle=\sum_{I=1}^{6}f_{H}M_{D,Ia}c_{\varphi}U^{\nu}_{(I+3)i}\simeq
f_{H}c_{\varphi}\times\left[\begin{array}[]{cc}0,&\quad i\leq 3\\\
\left(M^{T}_{D}V\right)_{a(i-3)},&\quad i>3\end{array}\right.,$ (40)
$\displaystyle\lambda^{L,2}_{ia}$
$\displaystyle\simeq-\lambda^{L,1}_{ia}t_{\varphi}=\lambda^{L,1}_{ia}\left[c_{\varphi}\to-
s_{\varphi}\right],$ $\displaystyle\lambda^{R,1}_{ia}$
$\displaystyle=m_{e_{a}}t_{\beta}c_{\varphi}U^{\nu*}_{ai}-\sum_{I=1}^{6}\frac{v}{\sqrt{2}}s_{\varphi}Y^{h}_{Ia}U^{\nu*}_{(I+3)i}$
$\displaystyle\simeq\left[\begin{array}[]{cc}m_{e_{a}}t_{\beta}c_{\varphi}\left[U^{*}_{\mathrm{PMNS}}\left(I_{3}-\frac{1}{2}\hat{x}_{\nu}\right)\right]_{ai}+\frac{vs_{\varphi}}{\sqrt{2}}\left(Y^{hT}R^{T}U^{*}_{\mathrm{PMNS}}\right)_{ai},&\quad
i\leq 3\\\
m_{e_{a}}t_{\beta}c_{\varphi}\left(RV\right)^{*}_{a(i-3)}-\frac{vs_{\varphi}}{\sqrt{2}}\left[Y^{hT}\left(I_{6}-\frac{1}{2}R^{T}R^{*}\right)V^{*}\right]_{a(i-3)}&\quad
i>3\end{array}\right.,$ (43) $\displaystyle\lambda^{R,2}_{ia}$
$\displaystyle=-m_{e_{a}}t_{\beta}s_{\varphi}U^{\nu*}_{ai}-\sum_{I=1}^{6}\frac{v}{\sqrt{2}}c_{\varphi}Y^{h}_{Ia}U^{\nu*}_{(I+3)i}$
$\displaystyle\simeq\left[\begin{array}[]{cc}-m_{e_{a}}t_{\beta}s_{\varphi}\left[U^{*}_{\mathrm{PMNS}}\left(I_{3}-\frac{1}{2}\hat{x}_{\nu}\right)\right]_{ai}+\frac{vc_{\varphi}}{\sqrt{2}}\left(Y^{hT}R^{T}U^{*}_{\mathrm{PMNS}}\right)_{ai},&\quad
i\leq 3\\\
-m_{e_{a}}t_{\beta}s_{\varphi}\left(RV\right)^{*}_{a(i-3)}-\frac{vc_{\varphi}}{\sqrt{2}}\left[Y^{hT}\left(I_{6}-\frac{1}{2}R^{T}R^{*}\right)V^{*}\right]_{a(i-3)}&\quad
i>3\end{array}\right.$ (46)
$\displaystyle=\lambda^{R,1}_{ia}\left[s_{\varphi}\to
c_{\varphi},\;c_{\varphi}\to-s_{\varphi}\right].$ (47)

The Feynman diagrams giving one-loop contributions to $(g-2)_{e_{a}}$
corresponding to the Lagrangian (37) are shown in Fig. 1.

Figure 1: One-loop contributions of $W^{\pm}$ and $h^{\pm}_{k}$ to
$(g-2)_{e_{a}}$ in the unitary gauge, where $k=1,2$.

We do not list here the couplings of neutral gauge and Higgs bosons because
they give suppressed contributions to $a^{\mathrm{NP}}_{e_{a}}$. In
particular, the relevant couplings are only the ones with usual charged
leptons $s^{0}\overline{e_{a}}e_{a}$ and
$Z_{\mu}\overline{e_{a}}\gamma^{\mu}e_{a}$. The one-loop level contribution
arising from the $Z$ gauge boson exchange is the same as the predicted by the
SM. The contributions arising from new neutral Higgs bosons are not larger
than the one coming from the SM-like Higgs boson since they are suppressed by
a factor of the order of $\mathcal{O}(10^{-14})$, because we assume here their
masses are at the TeV scale.

We will use the approximation that $m^{2}_{n_{i}}/m_{W}^{2}=0$ with $i\leq 3$
and $m^{2}_{n_{i}}/m_{W}^{2}=M_{0}^{2}/m_{W}^{2}=x_{W}$ with $i>3$. Then, the
contribution to $a_{e_{a}}$ arising from the $W$ exchange has the form:

$\displaystyle a_{e_{a}}(W)$
$\displaystyle=-\frac{g^{2}m^{2}_{e_{a}}}{8\pi^{2}m^{2}_{W}}\left[-\frac{5}{12}+\left(R^{*}R^{T}\right)_{aa}\times\left(\tilde{f}_{V}\left(x_{W}\right)+\frac{5}{12}\right)\right],$
(48)

where

$\displaystyle\tilde{f}_{V}(x)$
$\displaystyle=\frac{-4x^{4}+49x^{3}-78x^{2}+43x-10-18x^{3}\ln
x}{24(x-1)^{4}},$
$\displaystyle\tilde{f}_{V}(0)=-\frac{5}{12}\leq\tilde{f}_{V}(x)\leq\tilde{f}_{V}(\infty)=-\frac{1}{6}.$
(49)

Because $|\tilde{f}_{V}\left(x_{W}\right)+\frac{5}{12}|\leq\frac{5}{12}$ and
$\left(R^{*}R^{T}\right)_{aa}\leq\mathcal{O}(10^{-3})\ll 1$, Eq. (48) equals
to the one-loop level contribution predicted by the SM, see example in Ref.
Jegerlehner:2009ry :

$a^{(1)\mathrm{SM}}_{\mu}(W)=\frac{g^{2}m^{2}_{\mu}}{8\pi^{2}m^{2}_{W}}\times\frac{5}{12}\simeq
383\times 10^{-11},\;\frac{g^{2}m^{2}_{\mu}}{8\pi^{2}m^{2}_{W}}\simeq
9.19\times 10^{-9}.$ (50)

The one-loop level contribution to $\Delta a_{e_{a}}$ arising from the
exchange of the electrically charged scalar singlet $h^{\pm}_{k}$ is given by
Crivellin:2018qmi :

$\displaystyle a_{e_{a}}(h^{\pm}_{k})$
$\displaystyle=-\frac{g^{2}m_{e_{a}}\;}{8\pi^{2}m^{2}_{W}}\sum_{i=1}^{9}\frac{\lambda^{L,k*}_{ia}\lambda^{R,k}_{ia}m_{n_{i}}f_{\Phi}(x_{i,k})+m_{e_{a}}\left(\lambda^{L,k*}_{ia}\lambda^{L,k}_{ia}+\lambda^{R,k*}_{ia}\lambda^{R,k}_{ia}\right)\tilde{f}_{\Phi}(x_{i,k})}{m^{2}_{h^{\pm}_{k}}},$
(51)

where $x_{i,k}\equiv m^{2}_{n_{i}}/m^{2}_{h^{\pm}_{k}}$ and the loop functions
appearing in Eq. (51) have the forms:

$\displaystyle f_{\Phi}(x)$ $\displaystyle=\frac{x^{2}-1-2x\ln
x}{4(x-1)^{3}},\quad\tilde{f}_{\Phi}(x)=\frac{2x^{3}+3x^{2}-6x+1-6x^{2}\ln
x}{24(x-1)^{4}},$ $\displaystyle f_{\Phi}(\infty)=0\leq f_{\Phi}(x)\leq
f_{\Phi}(0)=\frac{1}{4},\;\tilde{f}_{\Phi}(\infty)=0\leq\tilde{f}_{\Phi}(x)\leq\tilde{f}_{\Phi}0)=\frac{1}{24}.$
(52)

And the deviation from the SM is defined as follows:

$\displaystyle\Delta a_{e_{a}}=\sum_{k=1}^{2}a_{e_{a}}(h^{\pm}_{k})+\Delta
a_{e_{a}}(W),\;\Delta a_{e_{a}}(W)\equiv
a_{e_{a}}(W)-a^{(1)\mathrm{SM}}_{e_{a}}(W),$ (53)

where $a^{(1)\mathrm{SM}}_{\mu}(W)\simeq 3.83\times 10^{-9}$
Jegerlehner:2009ry .

Using the approximations $m^{2}_{n_{i}}/m^{2}_{h^{\pm}_{k}}\simeq 0$ for
$i\leq 3$ and $m^{2}_{n_{i}}/m^{2}_{h^{\pm}_{k}}\simeq
M_{0}^{2}/m^{2}_{h^{\pm}_{k}}=x_{k}$, we have $f_{\Phi}(x_{i,k})\simeq
f_{\Phi}(0)$ for $i\leq 3$ and $f_{\Phi}(x_{i,k})\simeq f_{\Phi}(x_{k})$ for
$i>3$. Following Eqs. (15) and (16), we obtain that the one loop level
contribution to $a_{e_{a}}$ due to the exchange of $h^{\pm}_{1}$ is given by

$\displaystyle a_{e_{a}}(h^{\pm}_{1})$ $\displaystyle=-9.19\times
10^{-9}\times\frac{m^{2}_{e_{a}}}{m^{2}_{\mu}}$
$\displaystyle\times\mathrm{Re}\left\\{\sum_{c=1}^{3}f_{H}\left[c^{2}_{\varphi}t_{\beta}\left|U_{\mathrm{PMNS},ac}\right|^{2}\frac{m_{n_{c}}}{\mu_{0}}-\frac{vc_{\varphi}s_{\varphi}}{\sqrt{2}m_{e_{a}}}U_{\mathrm{PMNS},ac}\left(\frac{m_{n_{c}}}{\mu_{0}}\right)^{1/2}Y^{h}_{(c+3)a}\right]x_{1}f_{\Phi}(x_{1})\right.$
$\displaystyle\quad+\sum_{c=1}^{3}\left[\left|U_{\mathrm{PMNS},ac}\right|^{2}\left(f_{H}^{2}c^{2}_{\varphi}\frac{m_{n_{c}}}{\mu_{0}}\right)\right]x_{1}\tilde{f}_{\Phi}(x_{1})$
$\displaystyle\quad+\frac{m^{2}_{e_{a}}t_{\beta}^{2}c_{\varphi}^{2}}{m^{2}_{H^{\pm}_{1}}}\left[\frac{1}{24}-\sum_{c=1}^{3}\left[\left|U_{\mathrm{PMNS},ac}\right|^{2}\frac{m_{n_{c}}}{\mu_{0}}\right]\left(\frac{1}{24}-\tilde{f}_{\Phi}(x_{1})\right)\right]$
$\displaystyle\quad+\frac{v^{2}s_{\varphi}^{2}}{2m^{2}_{H^{\pm}_{1}}}\left[\sum_{c=1}^{3}\left[\left|Y^{h}_{(c+3)a}\right|^{2}\frac{m_{n_{c}}}{\mu_{0}}\right]\left(\frac{1}{24}-\tilde{f}_{\Phi}(x_{1})\right)+\left(Y^{h\dagger}Y^{h}\right)_{aa}\tilde{f}_{\Phi}(x_{1})\right]$
$\displaystyle\quad-\left.\frac{vm_{e_{a}}t_{\beta}s_{2\varphi}}{\sqrt{2}m^{2}_{H^{\pm}_{1}}}\left(\frac{1}{24}-\tilde{f}_{\Phi}(x_{1})\right)\sum_{c=1}^{3}\left[U_{\mathrm{PMNS},ac}Y^{h}_{(c+3)a}\left(\frac{m_{n_{c}}}{\mu_{0}}\right)^{1/2}\right]\right\\},$
(54) $\displaystyle a_{e_{a}}(h^{\pm}_{2})$
$\displaystyle=a_{e_{a}}(h^{\pm}_{1})\left[x_{1}\to\;x_{2},\;s_{\varphi}\to
c_{\varphi},\;c_{\varphi}\to-s_{\varphi}\right].$ (55)

In the real part of Eq. (54), the first line corresponds to the chirally-
enhanced part proportional to $\lambda^{L,1*}\lambda^{R,1}$ whereas the second
and remaining lines are the parts proportional to
$\lambda^{L,1*}\lambda^{L,1}$ and $\lambda^{R,1*}\lambda^{R,1}$, respectively.

Now we compare our results given in Eq. (54) with the one-loop contribution
due to the exchange of singly electrically charged Higgs bosons in the
original versions Aoki:2009ha  without $N_{IR}$ and the singlet $h^{\pm}$. Now
we assume $h^{\pm}_{1}\equiv H^{\pm}$ in Eq. (III), corresponding to
$c_{\varphi}=1,s_{\varphi}=0$. The absence of $N_{IR}$ can be conveniently
derived from $f_{H}=0$ and $f(x_{1})=0$, thus implying that the only one-loop
contribution from $h^{\pm}_{1}$ only consists of the first term in the third
line of the real part given in Eq. (54), which is proportional to $-9.19\times
10^{-9}\times
m^{4}_{e_{a}}t_{\beta}^{2}c_{\varphi}^{2}/(24m^{2}_{\mu}m^{2}_{H^{\pm}_{1}})$,
which yields a small and negative contribution to $\Delta
a^{\mathrm{NP}}_{\mu}$ Jueid:2021avn . The dominant contributions to AMM arise
from two-loop Barr-Zee type diagrams. The similar conclusion for the THDM
type-B where $t^{2}_{\beta}\to t^{-2}_{\beta}$, hence has suppressed two-loop
contributions to AMM for $t_{\beta}\geq 0.4$.

If the mixing between two singly charged Higgs bosons vanishes, namely
$s_{\alpha}c_{\alpha}=0$, then all of the remaining terms in both Eqs. (54)
and (55) are negative, thus not allowing to accommodate the experimental data
on muon and electron anomalous magnetic moments.

Using the constraint (28) for
$RR^{\dagger}=U_{\mathrm{PMNS}}\hat{x}_{\nu}U_{\mathrm{PMNS}}^{\dagger}$ we
have $\hat{x}_{\nu}<\mathcal{O}(10^{-3})$. Therefore, we will choose a safe
upper bound as follows

$\hat{x}_{\nu
3}\equiv\left(\hat{x}_{\nu}\right)_{33}=\frac{m_{n_{3}}}{\mu_{0}}\leq
10^{-3}.$ (56)

To avoid unnecessary independent parameters of $Y^{h}_{Ia}$ without any
qualitative AMM results discussed on this work and in order to cancel large
one-loop contributions from these Higgs bosons to the cLFV decays $e_{b}\to
e_{a}\gamma$, we assume that

$\displaystyle\sum_{c=1}^{3}U_{\mathrm{PMNS},ac}\left(\frac{m_{n_{c}}}{m_{n_{3}}}\right)^{1/2}Y^{h}_{(c+3)b}=Y^{d}_{a}\delta_{ab}\,.$
(57)

The total one-loop level contribution arising from the exchange of two singly
charged Higgs bosons is written as

$\displaystyle a_{e_{a}}(h^{\pm})$ $\displaystyle\equiv
a_{e_{a}}(h^{\pm}_{1})+a_{e_{a}}(h^{\pm}_{2})=a_{e_{a},0}(h^{\pm})+\dots,$
(58) $\displaystyle a_{e_{a},0}(h^{\pm})$ $\displaystyle=9.19\times
10^{-9}f_{H}\times\frac{m_{e_{a}}^{2}}{m^{2}_{\mu}}\times\mathrm{Re}\left\\{\frac{vs_{2\varphi}}{2\sqrt{2}m_{e_{a}}}\hat{x}_{\nu
3}^{1/2}Y^{d}_{a}\left[x_{1}f_{\Phi}(x_{1})-x_{2}f_{\Phi}(x_{2})\right]\right\\},$
(59)

where $a_{e_{a},0}(h^{\pm})$ denotes the dominant term of the chirally-
enhanced part coming from the second one in the first line of the real part
given in Eq. (54) where $x_{k}$ is the part relating with the contribution
from $h^{\pm}_{k}$ exchange. This conclusion can be qualitatively understood
from the property of the large factor $v/(2\sqrt{2}m_{e_{a}})$ as well as
large free Yukawa couplings up to the perturbative limit
max[$|Y^{h}_{Ia}|]\leq\sqrt{4\pi}$. In addition, the sign of this term can be
the same as $\Delta a^{\mathrm{NP}}_{e_{a}}$ depending on the sign of
Re[$Y^{d}_{a}$] when all other factors are fixed. As a result, this term can
easily explain both signs of $\Delta a^{\mathrm{NP}}_{\mu}$ and $\Delta
a^{\mathrm{NP}}_{e}$ that are still in conflict between different experimental
results and needed to be confirmed in the future. Our numerical investigation
showed that the term in Eq. (59) is the dominant one, and the sum of all the
remaining terms is suppressed in the allowed regions of parameter space. More
detailed estimations confirming that the remaining contributions are
suppressed were given in Ref. Hue:2021zyw . Finally, $a_{e_{a},0}(h^{\pm})\neq
0$ only when the mixing between two $SU(2)_{L}$ doublets and singlets is non-
zero $s_{\varphi}c_{\varphi}\neq 0$, and their masses are non degenerate
$m_{h^{\pm}_{1}}\neq m_{h^{\pm}_{2}}$.

We comment here an important property that $a_{e_{a},0}(h^{\pm})$ in Eq. (59)
keeps the same form for both types of THDM A and B, because these models
control only the first term of $\lambda^{R}_{ia}$ in Eq. (40), which is
proportional to $t_{\beta}$ or $t^{-1}_{\beta}$. This key factor controls loop
contributions to AMM, which are proportional to suppressed power of
$t_{\beta}^{-1}$ with large $t_{\beta}$ in THDM of type-B, including the model
type-I. Hence, it is impossible to accommodate the AMM data in the original
version of the THDM type I. On the other hand, the THDM type-A has one and two
loop contributions to AMM consisting of much enhanced factors of
$t^{2}_{\beta}$ and $t^{4}_{\beta}$, respectively. Therefore, original
versions can predict large loop contributions to AMM, provided that the new
Higgs bosons are light having low masses of about few hundred GeV. These might
be excluded by future collider experiments. Then, the presence of
$a_{e_{a},0}(h^{\pm})$ is an alternative way to explain the AMM data.

Now we consider the case of $f_{H}=t_{\beta}^{-1}$, which corresponds to the
models where the RH neutrino singlets couple with $\Phi_{2}$, whereas the
charged leptons couple with $\Phi_{1}$ Li:2020dbg ; Hue:2021zyw . Now,
increasing values of $a_{e_{a},0}(h^{\pm})$ in Eq. (59) require small
$t_{\beta}$ values, thus implying that the scanning range of $t_{\beta}$
should be chosen from the lower bound $t_{\beta}\geq 0.4$. Numerical
illustrations of $a_{e_{a},0}(h^{\pm})$ are shown in Fig. 2 with fixed
$\hat{x}_{\nu 3}=10^{-3}$, i.e., $0\neq|(Y^{N}_{2})_{ab}|<\sqrt{4\pi}$.

|
---|---
Figure 2:  Contour plots of $\Delta a_{\mu}(h^{\pm})\times 10^{9}$ and $\Delta
a_{e}(h^{\pm})\times 10^{13}$ as functions of $x_{1}$ and $x_{2}$, where
$\hat{x}_{\nu 3}=10^{-3}$ and $f_{H}=t_{\beta}^{-1}$. The green (orange)
regions are excluded by the 1$\sigma$ data of $\Delta a^{\mathrm{NP}}_{\mu}$
($\Delta a^{\mathrm{NP}}_{e}$). The black and dashed-black curves show the
constant values of $\Delta a_{e}(h^{\pm})\times 10^{13}$ and $\Delta
a_{\mu}(h^{\pm})\times 10^{9}$, respectively.

In addition, there are different fixed values of
$s_{\varphi},t_{\beta},Y^{d}_{1}$, and $Y^{d}_{2}$ shown in the respective
panels, namely small $t_{\beta}=0.5$ and large $t_{\beta}=10$. We have checked
that $a_{e_{a},0}(h^{\pm})\simeq a_{e_{a}}(h^{\pm})$. Now we estimate the
allowed ranges of $M_{0}$, $m_{h^{\pm}_{{1,2}}}$, which are affected from the
perturbative limit of the Yukawa coupling matrix $Y^{N}_{0}$ defined in Eq.
(35) and related with $M_{0}$ through Eqs. (13), (8), and (16). We have two
different constraints corresponding to $f_{H}=t^{-1}_{\beta}$ and
$(-t_{\beta})$. The constraint for $f_{H}=t^{-1}_{\beta}$ is

$\displaystyle
M_{0}\left|\left(\hat{x}_{\nu}^{1/2}U^{\dagger}_{\mathrm{PMNS}}\right)_{ab}\right|=\frac{vt_{\beta}}{\sqrt{2}}\left|c_{\beta}(Y^{N}_{2})_{(a+3)b}\right|<\frac{vs_{\beta}}{\sqrt{2}}\sqrt{4\pi}.$
(60)

We can choose a more strict upper bound of $M_{0}\leq\hat{x}_{\nu
3}^{-1/2}s_{\beta}v\sqrt{\pi/2}=9.7s_{\beta}$ TeV for $f_{H}=t_{\beta}^{-1}$
and $\hat{x}_{\nu 3}=10^{-3}$. Therefore for $\hat{x}_{\nu 3}\leq 10^{-3}$ and
$t_{\beta}\geq 0.4$, values of $M_{0}\leq 6.4$ TeV are always acceptable.
Applying this constraint, we consider a benchmark where $M_{0}=3$ TeV in order
to estimate the allowed values of $m_{h^{\pm}_{1,2}}$. We can see that in the
left-panel of figure 2, the range $0.2\leq
x_{2}=M_{0}^{2}/m^{2}_{h^{\pm}_{2}}\leq 1.5$ with
$x_{1}=M_{0}^{2}/m^{2}_{h^{\pm}_{1}}=10$ is allowed with respect to
$2.5\;\mathrm{TeV}\leq m_{h^{\pm}_{2}}\leq 5$ TeV and $m_{h^{\pm}_{1}}\simeq
0.95$ TeV. Similarly for the right panel of figure 2, we can choose $M_{0}=$ 5
TeV and larger $m_{h^{\pm}_{2}}=M_{0}/\sqrt{2.5}$ and $m_{h^{\pm}_{1}}\geq
5/\sqrt{15}>1$ TeV.

In the last discussion we will focus on the allowed regions consisting of
masses of heavy neutrinos and singly charged Higgs bosons below few TeV so
that they can be detected by future colliders. The allowed regions are defined
as they result in the two values of $\Delta a_{\mu}(h^{\pm})$ and $\Delta
a_{e}(h^{\pm})$ both satisfying the experimental data of the muon and electron
anomalous magnetic moments within the 1 $\sigma$ level, and all perturbative
limits of the Yukawa couplings $Y^{h}_{Ia}$ and $Y^{N}_{Ia}$ are satisfied,
namely $|Y^{h}_{Ia}|,|Y^{N}_{2,Ia}|<\sqrt{4\pi}$. The region of parameter
space used to scan is chosen as follows:

$\displaystyle 0.8\;\mathrm{TeV}\leq m_{h^{\pm}_{1}},m_{h^{\pm}_{2}}\leq
5\;\mathrm{TeV};10\;\mathrm{GeV}\leq M_{0}\leq\;5\;\mathrm{TeV};$
$\displaystyle 0.3\leq t_{\beta}\leq 50;\;\left|s_{\varphi}\right|\leq
1;\;|Y^{d}_{ab}|\leq 3;\;10^{-7}\leq\;\hat{x}_{\nu
3}\equiv\left(\hat{x}_{\nu}\right)_{33}=\frac{m_{n_{3}}}{\mu_{X}}\leq
10^{-3}.$ (61)

Here we fix $m_{n_{1}}=0.01$ eV corresponding to the NO scheme used in our
numerical analysis. Smaller values of $m_{n_{1}}$ will result in small allowed
ranges of $Y^{d}_{1,2}$ because of the perturbative limit affecting the
relation given in Eq. (57). The scanned range of $\hat{x}_{\nu 3}$ satisfies
the non-unitary constraint given in Eq. (28). The numerical results confirm
that $\left|a_{e_{a},0}(h^{\pm})/a_{e_{a}}(h^{\pm})\right|\simeq 1$, namely
$0.995<\left|a_{e,0}(h^{\pm})/a_{e}(h^{\pm})\right|<1.005$ and
$1\leq\left|a_{\mu,0}(h^{\pm})/a_{\mu}(h^{\pm})\right|\leq 1.03$. Hence the
discussion about correlations between different contributions in Eq. (55) will
not be shown. The correlations between important free parameters vs. $\Delta
a_{\mu}(h^{\pm})$ are shown in Fig. 3.

|
---|---
Figure 3:  The correlations of $\Delta a_{\mu}(h^{\pm})$ vs $t_{\beta}$ and
$Y^{d}_{1,2}$ with $f_{H}=t_{\beta}^{-1}$.

As mentioned above, large $t_{\beta}$ favors small $\Delta a_{\mu}(h^{\pm})$,
leading to an upper bound on $t_{\beta}$ in the allowed regions, see the left-
panel of Fig. 3. The scanned ranges in Eq. (III) allow all experimental ranges
of $\Delta a_{e,\mu}$. In addition, the dependence between $\Delta a_{e}$ and
$\Delta a_{\mu}$ is not interesting. We know instead that $\Delta
a_{e}(h^{\pm})\simeq\Delta a_{e,0}(h^{\pm})$ is a function of the Yukawa
coupling $Y^{d}_{1}$. Hence, the dependence of $\Delta a_{e}(h^{\pm})$ on
$\Delta a_{\mu}(h^{\pm})$ can be seen from the dependence of the Yukawa
coupling $Y^{d}_{1}$ in the right panel, where it is bounded in a more
restrictive range than the one given in Eq. (III), see Table 2, where other
allowed ranges are also listed.

| $t_{\beta}$ | $|s_{\varphi}|$ |  $M_{0}$ [GeV] | $\left|Y_{1}^{d}\right|$ | $\left|Y_{2}^{d}\right|$ | $\hat{x}_{\nu 3}$
---|---|---|---|---|---|---
Min | 0.4 | 0.03 | 282 | 0.003 | 0.171 | $2.92\times 10^{-7}$
Max | 33.7 | 0.999 | 5000 | 0.203 | 2.529 | $10^{-3}$
Table 2: Allowed ranges corresponding to the scanned region given in Eq.
(III), considered in case of $f_{H}=t_{\beta}^{-1}$.

As we mentioned above, the dominant contributions to $\Delta a_{e_{a}}$ are
$\Delta a_{e_{a},0}$ given in Eq. (59). This property can be seen in Fig. 4,
showing the dependence of the ratio $Y^{d}_{1}/Y^{d}_{2}$ and $Y^{d}_{1,2}$ on
$\Delta a_{e}$.

|
---|---
Figure 4:  $\frac{Y^{d}_{1}}{Y^{d}_{d2}}$ (left panel) and $Y^{d}_{1}$ (right
panel) as functions of $\Delta a_{e}$ with $f_{H}=t_{\beta}^{-1}$.

The allowed region of this ratio $Y^{d}_{1}/Y^{d}_{2}\sim\Delta
a_{e,0}(h^{\pm})/\Delta a_{\mu,0}(h^{\pm})$ also linearly depends on $\Delta
a_{e}(h^{\pm})$. The vertical width of the allowed region is controlled by
both $1\sigma$ ranges of $\Delta a_{e}(h^{\pm})$ and $\Delta
a_{\mu}(h^{\pm})$. Also, the right panel shows the linear dependence of the
allowed region of $Y^{d}_{1}$ on $\Delta a_{e}(h^{\pm})$. The linear behavior
of $Y^{d}_{2}$ is less clear than the one of $Y^{d}_{1}$ because $Y^{d}_{2}$
is also affected by the perturbative condition.

An interesting property is that excepting $t_{\beta}$, all other parameters
like $M_{0}$, the non-unitary parameter $\hat{x}_{\nu 3}$, $Y^{d}_{1,2}$, and
$s_{\varphi}$ must have lower bounds. The allowed ranges of heavy neutrino
masses $M_{0}\geq\mathcal{O}(100)$ GeV might be confirmed by recent
experimental searches in colliders such as LHC and ILC Das:2012ze ;
Das:2014jxa ; Das:2015toa ; Das:2016hof ; Das:2018usr . Because of the
sizeable mixing angle $\sim\sqrt{x_{\nu,3}}$ between ISS and active neutrinos
$\nu_{aL}$, the main production channel of heavy neutrinos $n_{I}$ (I=4,…,9)
with mass $M_{0}$ at the LHC is via the Drell Yan anihilation process
$u\bar{d}\to n_{I}e^{+}_{a}$ mediated by the $W$ gauge boson in the $s$
channel. Then the decay channel of $n_{I}$ can be $n_{I}\to
e_{a}^{-}W^{+},n_{a}Z,\;n_{a}h$, where $h$ is the standard model-like Higgs
boson. The ILC can produce heavy neutrino in the processes
$e^{+}e^{-}\to\bar{n}_{a}n_{I}$ through the exchange of virtual $W$ and $Z$
bosons in the $t$ and $s$-channels, respectively. The model under
consideration also predicts the production channel of a heavy neutrino pair
$e^{+}e^{-}\to\bar{n}_{I}n_{I}$ through the virtual exchange of $h^{\pm}_{k}$.
In addition, the singly charged Higgs bosons in the model under consideration
can be searched in a proton-proton collider through the processes
$pp\to\gamma^{*}/Z^{*}\to
h^{+}_{k}h^{-}_{l}\to(e^{+}_{c}n_{a})(\bar{n}_{b}e^{-}_{d})$ (with
$a,b,c,d=1,2,3$), where the Yukawa couplings $Y^{h}$ give an important
contribution Calle:2021tez . The ILC can produce two singly charged Higgs
bosons $e^{+}e^{-}\to h^{+}_{k}h^{-}_{l}$ through the $n_{I}$ exchange in the
$t$-channel. Studying these processes are beyond the scope of this work, but
will be investigated in more detail elsewhere. Because of the non-vanishing
mixing between $h^{\pm}$ and the singly charged components of the Higgs
doublets, another decay into a CP-odd neutral Higgs boson $A$, such as
$h^{\pm}_{k}\to AW^{\pm}$, can occur Rose:2021cav .

The correlations relating the masses with $\Delta a_{\mu}(h^{\pm})$ are shown
in Fig. 5.

|
---|---
Figure 5:  The correlations of $\Delta a_{\mu}(h^{\pm})$ and $m_{h^{\pm}_{1}}$
vs. masses $M_{0}$ and $m_{h^{\pm}_{1,2}}$ with $f_{H}=t_{\beta}^{-1}$.

We can see that $m^{2}_{h^{\pm}_{1}}$ must be different than
$m^{2}_{h^{\pm}_{2}}$ and our numerical analysis indicates
$\left|m^{2}_{h^{\pm}_{2}}-m^{2}_{h^{\pm}_{1}}\right|\geq 100$ GeV. The
allowed regions of large $m_{h^{\pm}_{1,2}}$ at TeV scale corresponding to
$(g-2)_{\mu}$ experimental data may be tested indirectly through the process
$\mu^{+}\mu^{-}\to\gamma^{*}\to h\gamma$ at multi-TeV muon colliders
Yin:2020afe .

Finally, the correlations showing significant dependence of free parameters
and $t_{\beta}$ are given in Fig. 6,

|
---|---
Figure 6:  The correlations of $t_{\beta}$ vs other free parameters in the
allowed regions with $f_{H}=t_{\beta}^{-1}$.

where the allowed regions with large $t_{\beta}\geq 20$ require both
conditions of large mixing $|s_{2\varphi}|=1$ and large $\hat{x}_{\nu 3}$. We
obtain a small allowed range of $t_{\beta}\leq 10$ that was missed in Ref.
Mondal:2021vou . Our result is consistent with the discussion corresponding to
3-3-1 models given in Ref. Hue:2021xap , where the THDM is embedded.

For the case of $f_{H}=(-t_{\beta})$, i.e., RH neutrino and charged leptons
singlets couple with the same $SU(2)_{L}$ Higgs doublet $\Phi_{1}$ in Yukawa
Lagrangian (II), for example Ref. Hue:2021xap . Then large
$a_{e_{a}}(h^{\pm})$ in Eq. (59) support large $t_{\beta}$. In addition, even
when $a_{e_{a}}(h^{\pm})\sim s_{2\varphi}\neq 0$, the first terms in the first
lines of Eq. (54) or (55) may be large enough consistent with $\Delta
a^{\mathrm{NP}}_{\mu,e}$ in both sign and amplitude. But the simple
assumptions of the couplings and the total neutrino mass matrix in this work
are not enough to explain both $(g-2)_{e,\mu}$ experimental data. The
perturbative constraint gives an upper bound on $M_{0}$, namely

$\displaystyle
M_{0}\left|\left(\hat{x}_{\nu}^{1/2}U^{\dagger}_{\mathrm{PMNS}}\right)_{ab}\right|=\frac{vc_{\beta}}{\sqrt{2}}\left|(Y^{N}_{1})_{(a+3)b}\right|\Rightarrow
M_{0}<\hat{x}_{\nu,3}^{-1/2}\left|U_{\mathrm{PMNS},23}\right|^{-1}vc_{\beta}\sqrt{2\pi},$
(62)

therefore $M_{0}$ may be small with small $c_{\beta}$ equivalently with large
$t_{\beta}$. We always have $M_{0}<1.6$ TeV for $t_{\beta}=0.4$ and
$\hat{x}_{\nu,3}=10^{-3}$. Although $t_{\beta}\geq 0.4$ is always kept,
smaller $\hat{x}_{\nu,3}\leq 10^{-3}$ may allow large $M_{0}$ which also allow
large $m_{h^{\pm}_{1,2}}$. On the other hand, from Eq. (59), $|a_{e_{a}}|\sim
t_{\beta}\left(x_{\nu 3}\right)^{1/2}|s_{2\varphi}Y^{d}_{a}|$, hence the
allowed regions are easily satisfied for small $t_{\beta}$ and large values of
other parameters. Our numerical analysis shows that large $x_{\nu 3}\geq
10^{-4}$ allows all $t_{\beta}\leq 50$ and all singly charged Higgs masses at
the TeV scale. Contour plots of $\Delta a_{\mu}(h^{\pm})\times 10^{9}$ and
$\Delta a_{e}(h^{\pm})\times 10^{13}$ for $\hat{x}_{\nu,3}=5\times 10^{-4}$
with small $t_{\beta}=0.5$ and large $t_{\beta}=20$ are shown in Fig. 7.

|
---|---
Figure 7:  Contour plots of $\Delta a_{\mu}(h^{\pm})\times 10^{9}$ and $\Delta
a_{e}(h^{\pm})\times 10^{13}$ as functions of $x_{1}$ and $x_{2}$, where
$\hat{x}_{\nu 3}=5\times 10^{-4}$ and $f_{H}=-t_{\beta}$. The green (orange)
regions are excluded by the 1$\sigma$ data of $\Delta a^{\mathrm{NP}}_{\mu}$
($\Delta a^{\mathrm{NP}}_{e}$). The black and dashed-black curves show the
constant values of $\Delta a_{e}(h^{\pm})\times 10^{13}$ and $\Delta
a_{\mu}(h^{\pm})\times 10^{9}$, respectively.

Regarding the numerical analysis in the scanning range (III) with
$f_{H}=-t_{\beta}$, the allowed ranges are tighter than the scanning regions
as shown below

$\displaystyle|s_{\varphi}|>10^{-3},\;0.21\geq|Y^{d}_{1}|>0.002,\;2.6>|Y^{d}_{2}|>0.02,\;\hat{x}_{\nu
3}>3\times 10^{-10}.$ (63)

In addition small $M_{0}\leq 30$ GeV is also allowed. The numerical results of
the correlations between $t_{\beta}$ vs. $Y^{d}_{1}$, $Y^{d}_{2}$, and
$\hat{x}_{\nu 3}$ are shown in Fig. 8,

|
---|---
Figure 8:  The correlations of $t_{\beta}$ vs other free parameters in the
allowed regions with $f_{H}=-t_{\beta}$.

where the allowed ranges of $Y^{d}_{2}$ correspond to a rather narrow allowed
range of $Y^{d}_{1}$, see Eq. (63). This implies that the phenomenology of the
singly charged Higgs boson at colliders related with these two couplings will
have some certain relations that should be experimentally verified.

We have checked numerically that although lower bounds for allowed ranges of
$s_{\varphi}$ and $\hat{x}_{\nu 3}$ are tiny, but never vanishes. In addition,
small allowed values of $\hat{x}_{\nu 3}$ near lower bounds require both large
$t_{\beta}\to 50$ and $|s_{2\varphi}|=2|c_{\varphi}s_{\varphi}|\to 1$,
implying the existence of $\Delta a_{e_{a,0}}(h^{\pm})$. Also, small allowed
values of $|s_{2\varphi}|$ near the lower bound require both large
$\hat{x}_{\nu 3}\to 10^{-3}$ and $t_{\beta}$. Illustrations for these comments
are shown in Fig. 9.

|
---|---
|
Figure 9:  The correlations of $t_{\beta}$ vs $s_{2\varphi}$ (upper panels)
and $\hat{x}_{\nu 3}$ (lower panels) near lower bounds of their allowed ranges
with $f_{H}=-t_{\beta}$.

We see also that, the relation shown in Fig. 4 is also true with
$f_{H}=-t_{\beta}$.

Finally, we have comments regarding the $(g-2)_{e}$ data Parker:2018vye :
$\Delta a^{\mathrm{NP}}_{e}=-(8.7\pm 3.6)\times 10^{-13}$. The allowed regions
of parameter space corresponding to this data can be derived from the above
described numerical analysis. Namely, excepting $Y^{d}_{1}$ all allowed ranges
of free parameters are kept unchanged to guarantee consistency with the
experimental $(g-2)_{\mu}$ data. On the other hand, $Y^{d}_{1}$ is changed
into new values such that $Y^{d}_{1}\to\frac{(-8.7\pm 3.6)}{(4.8\pm 3)}\times
Y^{d}_{2}$ and exclude too large values violating the perturbative limit. The
two models under consideration predict allowed regions of parameter space that
are different from those discussed in Ref. DelleRose:2020oaa ; Rose:2021cav
for the THDM model where fermions couple to two Higgs doublets with the
aligned assumption of the two respective Yukawa couplings
$Y^{\prime}_{2f}=\zeta_{f}Y^{\prime}_{1f}$ with $f=u,d,\ell,\nu$, where $\nu$
denote new right handed neutrinos $\nu_{aR}$. The large contributions to
$(g-2)_{\mu}$ comes from the two-loop Barr-Zee contributions with the
necessary condition of very light neutral CP-odd Higgs boson mass, $m_{A}\leq
60$ GeV. On the other hand, large and negative sign of $(g-2)_{e}$ satisfying
the $1\sigma$ experimental data comes from one-loop contribution with the
condition that $\zeta_{\ell}\zeta_{\nu}<0$. The model in Ref.
DelleRose:2020oaa  does not include the case $Y^{\prime}_{1\nu}=0$ and
$Y^{\prime}_{2\nu}\neq 0$, corresponding to $f_{H}=t_{\beta}^{-1}$ mentioned
in our work. In addition, the region of parameter space with
$\zeta_{\nu}=\zeta_{\ell}=0$ corresponding to $f_{H}=-t_{\beta}$ is excluded
by the $1\sigma$ range of $(g-2)_{e}$ data. In contrast, our models always
assume that $Y^{\prime}_{1u}=\zeta_{\ell}=0$, and only one of the two Yukawa
coupling matrices $Y^{\prime}_{1\nu}$ or $Y^{\prime}_{2\nu}$ being non-zero.
The Yukawa couplings $Y^{d}_{a}$ between only gauge singlets
$N_{a(b+3)}e_{aR}h^{+}$ give main one-loop contributions to both $\Delta
a_{e,\mu}$, leading to nearly linear relations of these two quantities.
Finally, our models predict regions of parameter space that successfully
accommodate the experimental data on both $(g-2)_{e,\mu}$ anomalies without
the requirement of small $m_{A}$ and rather light
$m_{h^{\pm}_{k}}\sim\mathcal{O}(10^{2})$ GeV. Therefore, our models will be a
another solution for the $(g-2)_{e,\mu}$ anomalies if a light CP odd scalar
$A$ is excluded by future experiments.

##  IV  Conclusion

In this work we have shown that the appearance of heavy ISS neutrinos and
singly charged Higgs bosons is a very promoting solution to explain the
experimental data on both $(g-2)_{\mu,e}$ anomalies in many different types of
THDM, and in regions of parameter space allowing heavy singly charged Higgs
boson masses up to the TeV scale and small values of the $\tan\beta$ parameter
satisfying $\tan\beta\geq 0.4$. In particularly, the most important terms
$\Delta a_{e_{a},0}$ given in Eq. (59) are enough to explain successfully the
experimental AMM data of both $a_{e}$ and $a_{\mu}$ in THDM, including the
model type I, where other loop contributions to $\Delta a_{e_{a}}$ caused by
power of factor $t^{-1}_{\beta}$ are suppressed. In other types of THDM
needing large $t_{\beta}$ giving sizeable loop contributions to $\Delta
a_{e_{a}}$, light masses of new Higgs bosons around few hundred GeV in the
loop are also necessary to successfully accommodate the recent experimental
AMM data. The presence of $\Delta a_{e_{a},0}$ is an alternative way to
explain the AMM data if light mass ranges are excluded by future collider
searches. This solution will enlarge the allowed regions of parameter spaces
of the THDM which can simultaneously explain the $(g-2)_{e,\mu}$ data thanks
to the one loop exchange of electrically charged Higgs bosons with masses
within the LHC reach. The existence of $\Delta a_{e_{a},0}$ also yieds the
following consequences: i) the non-zero mixing $s_{2\varphi}\neq 0$, ii) the
non-unitary parameter $\hat{x}_{\nu,3}$ has a lower bound $\hat{x}_{\nu,3}\geq
4\times 10^{-7}$ for $f_{H}=t_{\beta}^{-1}$ and
$\hat{x}_{\nu,3}\geq\mathcal{O}(10^{-10})$ for $f_{H}=-t_{\beta}$, iii) and
lower bounds of new heavy neutrino masses are of the order of
$\mathcal{O}(100)$ GeV. Tiny but non vanishing values of $\hat{x}_{\nu,3}\sim
10^{-10}$ require very large $t_{\beta}$ and $|s_{2\varphi}|\to 1$. Despite
the large number of parameters, the model is economical with a small amount of
BSM fields, much lower than the corresponding to several models considered in
the literature. Besides that, apart from explaining the $(g-2)$ anomalies, the
model considered in this paper can feature interesting collider signatures
mainly related with electrically charged scalar and heavy neutrino production
at the LHC, that can be useful to test that theory at colliders.

## Acknowledgments

We thank Prof. Ray Volkas, Dr. Claudio Andrea, Dr. Lei Wang, and Dr. Wen Yin
for useful comments. L. T. Hue is thankful to Van Lang University. This
research is funded by the Vietnam National Foundation for Science and
Technology Development (NAFOSTED) under the grant number 103.01-2019.387 as
well as by ANID-Chile FONDECYT 1210378, ANID PIA/APOYO AFB180002, and Milenio-
ANID-ICN2019_044.

## References

  * (1) B. Abi et al. [Muon g-2], Phys. Rev. Lett. 126, 141801 (2021) [arXiv:2104.03281 [hep-ex]].
  * (2) G. W. Bennett et al. [Muon g-2], Phys. Rev. D 73, 072003 (2006) [arXiv:hep-ex/0602035 [hep-ex]].
  * (3) T. Aoyama, N. Asmussen, M. Benayoun, J. Bijnens, T. Blum, M. Bruno, I. Caprini, C. M. Carloni Calame, M. Cè and G. Colangelo, et al. Phys. Rept. 887, 1 (2020) [arXiv:2006.04822 [hep-ph]].
  * (4) A. Keshavarzi, D. Nomura and T. Teubner, Phys. Rev. D 97, no.11, 114025 (2018) [arXiv:1802.02995 [hep-ph]].
  * (5) G. Colangelo, M. Hoferichter and P. Stoffer, JHEP 02, 006 (2019) [arXiv:1810.00007 [hep-ph]].
  * (6) M. Hoferichter, B. L. Hoid and B. Kubis, JHEP 08, 137 (2019) [arXiv:1907.01556 [hep-ph]].
  * (7) M. Davier, A. Hoecker, B. Malaescu and Z. Zhang, Eur. Phys. J. C 80, no.3, 241 (2020) [erratum: Eur. Phys. J. C 80, no.5, 410 (2020)] [arXiv:1908.00921 [hep-ph]].
  * (8) A. Keshavarzi, D. Nomura and T. Teubner, Phys. Rev. D 101, no.1, 014029 (2020) [arXiv:1911.00367 [hep-ph]].
  * (9) A. Kurz, T. Liu, P. Marquard and M. Steinhauser, Phys. Lett. B 734, 144-147 (2014) [arXiv:1403.6400 [hep-ph]].
  * (10) K. Melnikov and A. Vainshtein, Phys. Rev. D 70, 113006 (2004) [arXiv:hep-ph/0312226 [hep-ph]].
  * (11) P. Masjuan and P. Sanchez-Puertas, Phys. Rev. D 95, no.5, 054026 (2017) [arXiv:1701.05829 [hep-ph]].
  * (12) G. Colangelo, M. Hoferichter, M. Procura and P. Stoffer, JHEP 04, 161 (2017) [arXiv:1702.07347 [hep-ph]].
  * (13) M. Hoferichter, B. L. Hoid, B. Kubis, S. Leupold and S. P. Schneider, JHEP 10, 141 (2018) [arXiv:1808.04823 [hep-ph]].
  * (14) A. Gérardin, H. B. Meyer and A. Nyffeler, Phys. Rev. D 100, no.3, 034520 (2019) [arXiv:1903.09471 [hep-lat]].
  * (15) J. Bijnens, N. Hermansson-Truedsson and A. Rodríguez-Sánchez, Phys. Lett. B 798, 134994 (2019) [arXiv:1908.03331 [hep-ph]].
  * (16) G. Colangelo, F. Hagelstein, M. Hoferichter, L. Laub and P. Stoffer, JHEP 03, 101 (2020) [arXiv:1910.13432 [hep-ph]].
  * (17) G. Colangelo, M. Hoferichter, A. Nyffeler, M. Passera and P. Stoffer, Phys. Lett. B 735, 90-91 (2014) [arXiv:1403.7512 [hep-ph]].
  * (18) T. Blum, N. Christ, M. Hayakawa, T. Izubuchi, L. Jin, C. Jung and C. Lehner, Phys. Rev. Lett. 124, no.13, 132002 (2020) [arXiv:1911.08123 [hep-lat]].
  * (19) T. Aoyama, M. Hayakawa, T. Kinoshita and M. Nio, Phys. Rev. Lett. 109, 111808 (2012) [arXiv:1205.5370 [hep-ph]].
  * (20) T. Aoyama, T. Kinoshita and M. Nio, Atoms 7, no.1, 28 (2019)
  * (21) A. Czarnecki, W. J. Marciano and A. Vainshtein, Phys. Rev. D 67, 073006 (2003) [erratum: Phys. Rev. D 73, 119901 (2006)] [arXiv:hep-ph/0212229 [hep-ph]].
  * (22) C. Gnendiger, D. Stöckinger and H. Stöckinger-Kim, Phys. Rev. D 88, 053005 (2013) [arXiv:1306.5546 [hep-ph]].
  * (23) M. Davier, A. Hoecker, B. Malaescu and Z. Zhang, Eur. Phys. J. C 77, no.12, 827 (2017) [arXiv:1706.09436 [hep-ph]].
  * (24) M. Davier, A. Hoecker, B. Malaescu and Z. Zhang, Eur. Phys. J. C 71, 1515 (2011) [erratum: Eur. Phys. J. C 72, 1874 (2012)] [arXiv:1010.4180 [hep-ph]].
  * (25) D. Hanneke, S. Fogwell and G. Gabrielse, Phys. Rev. Lett. 100, 120801 (2008) [arXiv:0801.1134 [physics.atom-ph]].
  * (26) R. H. Parker, C. Yu, W. Zhong, B. Estey and H. Müller, Science 360, 191 (2018) [arXiv:1812.04130 [physics.atom-ph]].
  * (27) L. Morel, Z. Yao, P. Cladé and S. Guellati-Khélifa, Nature 588, no.7836, 61-65 (2020)
  * (28) T. Aoyama, M. Hayakawa, T. Kinoshita and M. Nio, Phys. Rev. Lett. 109, 111807 (2012) [arXiv:1205.5368 [hep-ph]].
  * (29) S. Laporta, Phys. Lett. B 772, 232-238 (2017) [arXiv:1704.06996 [hep-ph]].
  * (30) T. Aoyama, T. Kinoshita and M. Nio, Phys. Rev. D 97, no.3, 036001 (2018) [arXiv:1712.06060 [hep-ph]].
  * (31) H. Terazawa, Nonlin. Phenom. Complex Syst. 21, no.3, 268-272 (2018)
  * (32) S. Volkov, Phys. Rev. D 100, no.9, 096004 (2019) [arXiv:1909.08015 [hep-ph]].
  * (33) A. Gérardin, Eur. Phys. J. A 57, no.4, 116 (2021) [arXiv:2012.03931 [hep-lat]].
  * (34) R. Dermisek and A. Raval, Phys. Rev. D 88, 013017 (2013) [arXiv:1305.3522 [hep-ph]].
  * (35) A. Crivellin, M. Hoferichter and P. Schmidt-Wellenburg, Phys. Rev. D 98, no.11, 113002 (2018) [arXiv:1807.11484 [hep-ph]].
  * (36) P. Escribano, J. Terol-Calvo and A. Vicente, Phys. Rev. D 103, no.11, 115018 (2021) [arXiv:2104.03705 [hep-ph]].
  * (37) A. Crivellin and M. Hoferichter, JHEP 07, 135 (2021) [arXiv:2104.03202 [hep-ph]].
  * (38) R. Dermisek, K. Hermanek and N. McGinnis, Phys. Rev. D 104, no.5, 055033 (2021) [arXiv:2103.05645 [hep-ph]].
  * (39) A. E. C. Hernández, S. F. King and H. Lee, Phys. Rev. D 103, no.11, 115024 (2021) [arXiv:2101.05819 [hep-ph]].
  * (40) E. J. Chun and T. Mondal, JHEP 11, 077 (2020) [arXiv:2009.08314 [hep-ph]].
  * (41) M. Frank and I. Saha, Phys. Rev. D 102, no.11, 115034 (2020) [arXiv:2008.11909 [hep-ph]].
  * (42) K. F. Chen, C. W. Chiang and K. Yagyu, JHEP 09, 119 (2020) [arXiv:2006.07929 [hep-ph]].
  * (43) P. M. Ferreira, B. L. Gonçalves, F. R. Joaquim and M. Sher, Phys. Rev. D 104, no.5, 053008 (2021) [arXiv:2104.03367 [hep-ph]].
  * (44) M. Endo and S. Mishima, JHEP 08, no.08, 004 (2020) [arXiv:2005.03933 [hep-ph]].
  * (45) C. Hati, J. Kriewald, J. Orloff and A. M. Teixeira, JHEP 07, 235 (2020) [arXiv:2005.00028 [hep-ph]].
  * (46) D. Borah, M. Dutta, S. Mahapatra and N. Sahu, Phys. Rev. D 105, no.1, 015029 (2022) [arXiv:2109.02699 [hep-ph]].
  * (47) H. Bharadwaj, S. Dutta and A. Goyal, JHEP 11, 056 (2021) [arXiv:2109.02586 [hep-ph]].
  * (48) A. S. De Jesus, S. Kovalenko, F. S. Queiroz, C. Siqueira and K. Sinha, Phys. Rev. D 102, no.3, 035004 (2020) [arXiv:2004.01200 [hep-ph]].
  * (49) D. Cogollo, Y. M. Oviedo-Torres and Y. S. Villamizar, Int. J. Mod. Phys. A 35, no.23, 2050126 (2020) [arXiv:2004.14792 [hep-ph]].
  * (50) B. D. Sáez and K. Ghorbani, Phys. Lett. B 823, 136750 (2021) [arXiv:2107.08945 [hep-ph]].
  * (51) I. Bigaran and R. R. Volkas, Phys. Rev. D 102, 075037 (2020) [arXiv:2002.12544 [hep-ph]].
  * (52) A. Crivellin, D. Mueller and F. Saturnino, Phys. Rev. Lett. 127, no.2, 021801 (2021) [arXiv:2008.02643 [hep-ph]].
  * (53) D. Zhang, JHEP 07, 069 (2021) [arXiv:2105.08670 [hep-ph]].
  * (54) W. Y. Keung, D. Marfatia and P. Y. Tseng, LHEP 2021, 209 (2021) [arXiv:2104.03341 [hep-ph]].
  * (55) T. Mondal and H. Okada, Nucl. Phys. B 976, 115716 (2022) [arXiv:2103.13149 [hep-ph]].
  * (56) D. W. Kang, J. Kim and H. Okada, Phys. Lett. B 822, 136666 (2021) [arXiv:2107.09960 [hep-ph]].
  * (57) A. E. Cárcamo Hernández, C. Espinoza, J. Carlos Gómez-Izquierdo and M. Mondragón, [arXiv:2104.02730 [hep-ph]].
  * (58) S. P. Li, X. Q. Li, Y. Y. Li, Y. D. Yang and X. Zhang, JHEP 01, 034 (2021) [arXiv:2010.02799 [hep-ph]].
  * (59) L. Delle Rose, S. Khalil and S. Moretti, Phys. Lett. B 816, 136216 (2021) [arXiv:2012.06911 [hep-ph]].
  * (60) F. J. Botella, F. Cornet-Gomez and M. Nebot, Phys. Rev. D 102, no.3, 035023 (2020) [arXiv:2006.01934 [hep-ph]].
  * (61) X. F. Han, T. Li, H. X. Wang, L. Wang and Y. Zhang, Phys. Rev. D 104, no.11, 115001 (2021) [arXiv:2104.03227 [hep-ph]].
  * (62) C. H. Chen, C. W. Chiang and T. Nomura, Phys. Rev. D 104, no.5, 055011 (2021) [arXiv:2104.03275 [hep-ph]].
  * (63) X. F. Han, T. Li, L. Wang and Y. Zhang, Phys. Rev. D 99, no.9, 095034 (2019) [arXiv:1812.02449 [hep-ph]].
  * (64) M. Aoki, S. Kanemura, K. Tsumura and K. Yagyu, Phys. Rev. D 80, 015017 (2009) [arXiv:0902.4665 [hep-ph]].
  * (65) A. Jueid, J. Kim, S. Lee and J. Song, Phys. Rev. D 104, no.9, 095008 (2021) [arXiv:2104.10175 [hep-ph]].
  * (66) P. Athron, C. Balázs, D. H. Jacob, W. Kotlarski, D. Stöckinger and H. Stöckinger-Kim, JHEP 09, 080 (2021) [arXiv:2104.03691 [hep-ph]].
  * (67) M. Badziak and K. Sakurai, JHEP 10 (2019), 024 [arXiv:1908.03607 [hep-ph]].
  * (68) S. Li, Y. Xiao and J. M. Yang, Eur. Phys. J. C 82, no.3, 276 (2022) [arXiv:2107.04962 [hep-ph]].
  * (69) M. Endo and W. Yin, JHEP 08, 122 (2019) [arXiv:1906.08768 [hep-ph]].
  * (70) B. Dutta, S. Ghosh and T. Li, Phys. Rev. D 102, no.5, 055017 (2020) [arXiv:2006.01319 [hep-ph]].
  * (71) A. E. Cárcamo Hernández, S. Kovalenko, R. Pasechnik and I. Schmidt, JHEP 06, 056 (2019) [arXiv:1901.02764 [hep-ph]].
  * (72) A. E. C. Hernández, D. T. Huong and I. Schmidt, Eur. Phys. J. C 82, no.1, 63 (2022) [arXiv:2109.12118 [hep-ph]].
  * (73) A. E. Cárcamo Hernández, S. Kovalenko, R. Pasechnik and I. Schmidt, Eur. Phys. J. C 79, no.7, 610 (2019) [arXiv:1901.09552 [hep-ph]].
  * (74) R. Adhikari, I. A. Bhat, D. Borah, E. Ma and D. Nanda, Phys. Rev. D 105, no.3, 035006 (2022) [arXiv:2109.05417 [hep-ph]]..
  * (75) J. Herrero-García, T. Ohlsson, S. Riad and J. Wirén, JHEP 04, 130 (2017) [arXiv:1701.05345 [hep-ph]].
  * (76) S. Davidson and H. E. Haber, Phys. Rev. D 72, 035004 (2005) [erratum: Phys. Rev. D 72, 099902 (2005)] [arXiv:hep-ph/0504050 [hep-ph]].
  * (77) L. Allwicher, P. Arnan, D. Barducci and M. Nardecchia, JHEP 10, 129 (2021) [arXiv:2108.00013 [hep-ph]].
  * (78) A. Ibarra, E. Molinaro and S. T. Petcov, JHEP 09, 108 (2010) [arXiv:1007.2378 [hep-ph]].
  * (79) P. A. Zyla et al. [Particle Data Group], PTEP 2020, 083C01 (2020)
  * (80) E. Arganda, M. J. Herrero, X. Marcano and C. Weiland, Phys. Rev. D 91, 015001 (2015) [arXiv:1405.4300 [hep-ph]].
  * (81) J. A. Casas and A. Ibarra, Nucl. Phys. B 618, 171 (2001) [arXiv:hep-ph/0103065 [hep-ph]].
  * (82) A. M. Baldini et al. [MEG], Eur. Phys. J. C 76, no.8, 434 (2016) [arXiv:1605.05081 [hep-ex]].
  * (83) B. Aubert et al. [BaBar], Phys. Rev. Lett. 104, 021802 (2010) [arXiv:0908.2381 [hep-ex]].
  * (84) N. Aghanim et al. [Planck], Astron. Astrophys. 641, A6 (2020) [erratum: Astron. Astrophys. 652, C4 (2021)] [arXiv:1807.06209 [astro-ph.CO]]
  * (85) J. P. Pinheiro, C. A. de S. Pires, F. S. Queiroz and Y. S. Villamizar, Phys. Lett. B 823, 136764 (2021) [arXiv:2107.01315 [hep-ph]].
  * (86) E. Fernandez-Martinez, J. Hernandez-Garcia and J. Lopez-Pavon, JHEP 08, 033 (2016) [arXiv:1605.08774 [hep-ph]].
  * (87) N. R. Agostinho, G. C. Branco, P. M. F. Pereira, M. N. Rebelo and J. I. Silva-Marcos, Eur. Phys. J. C 78, no.11, 895 (2018) [arXiv:1711.06229 [hep-ph]].
  * (88) T. N. Dao, M. Mühlleitner and A. V. Phan, Eur. Phys. J. C 82, no.8, 667 (2022) [arXiv:2108.10088 [hep-ph]].
  * (89) C. Biggio, E. Fernandez-Martinez, M. Filaci, J. Hernandez-Garcia and J. Lopez-Pavon, JHEP 05, 022 (2020) [arXiv:1911.11790 [hep-ph]].
  * (90) A. M. Coutinho, A. Crivellin and C. A. Manzari, Phys. Rev. Lett. 125, no.7, 071802 (2020) [arXiv:1912.08823 [hep-ph]].
  * (91) C. A. Manzari, A. M. Coutinho and A. Crivellin, PoS LHCP2020, 242 (2021) [arXiv:2009.03877 [hep-ph]].
  * (92) F. Jegerlehner and A. Nyffeler, Phys. Rept. 477, 1-110 (2009) [arXiv:0902.3360 [hep-ph]].
  * (93) A. Nepomuceno and B. Meirose, Phys. Rev. D 101, 035017 (2020) [arXiv:1911.12783 [hep-ph]].
  * (94) A. Das and N. Okada, Phys. Rev. D 88, 113001 (2013) [arXiv:1207.3734 [hep-ph]].
  * (95) A. Das, P. S. Bhupal Dev and N. Okada, Phys. Lett. B 735, 364-370 (2014) [arXiv:1405.0177 [hep-ph]].
  * (96) A. Das and N. Okada, Phys. Rev. D 93, no.3, 033003 (2016) [arXiv:1510.04790 [hep-ph]].
  * (97) A. Das, P. Konar and S. Majhi, JHEP 06, 019 (2016) [arXiv:1604.00608 [hep-ph]].
  * (98) A. Das, S. Jana, S. Mandal and S. Nandi, Phys. Rev. D 99, no.5, 055030 (2019) [arXiv:1811.04291 [hep-ph]].
  * (99) W. Yin and M. Yamaguchi, Phys. Rev. D 106 (2022) no.3, 033007 [arXiv:2012.03928 [hep-ph]].
  * (100) L. T. Hue, K. H. Phan, T. P. Nguyen, H. N. Long and H. T. Hung, Eur. Phys. J. C 82, no.8, 722 (2022) [arXiv:2109.06089 [hep-ph]].
  * (101) L. T. Hue, H. T. Hung, N. T. Tham, H. N. Long and T. P. Nguyen, Phys. Rev. D 104, 033007 (2021) [arXiv:2104.01840 [hep-ph]].
  * (102) J. Schechter and J. W. F. Valle, Phys. Rev. D 25, 774 (1982)
  * (103) J. G. Korner, A. Pilaftsis and K. Schilcher, Phys. Rev. D 47, 1080-1086 (1993) [arXiv:hep-ph/9301289 [hep-ph]].
  * (104) W. Grimus and L. Lavoura, JHEP 11, 042 (2000) [arXiv:hep-ph/0008179 [hep-ph]].
  * (105) E. Ma, Phys. Rev. D 73, 077301 (2006) [arXiv:hep-ph/0601225 [hep-ph]].
  * (106) G. Degrassi, S. Di Vita, J. Elias-Miro, J. R. Espinosa, G. F. Giudice, G. Isidori and A. Strumia, JHEP 08 (2012), 098 [arXiv:1205.6497 [hep-ph]].
  * (107) S. Jangid and P. Bandyopadhyay, Eur. Phys. J. C 80 (2020) no.8, 715 [arXiv:2003.11821 [hep-ph]].
  * (108) S. Jangid, P. Bandyopadhyay, P. S. Bhupal Dev and A. Kumar, JHEP 08 (2020), 154 [arXiv:2001.01764 [hep-ph]].
  * (109) C. Coriano, L. Delle Rose and C. Marzo, JHEP 02 (2016), 135 [arXiv:1510.02379 [hep-ph]].
  * (110) L. Delle Rose, C. Marzo and A. Urbano, JHEP 12 (2015), 050 [arXiv:1506.03360 [hep-ph]].
  * (111) P. Bandyopadhyay, S. Jangid and M. Mitra, JHEP 02 (2021), 075 [arXiv:2008.11956 [hep-ph]].
  * (112) L. D. Rose, S. Khalil and S. Moretti, [arXiv:2111.12185 [hep-ph]].
  * (113) J. de Blas, EPJ Web Conf. 60 (2013), 19008 [arXiv:1307.6173 [hep-ph]].
  * (114) P. Bandyopadhyay, E. J. Chun, H. Okada and J. C. Park, JHEP 01, 079 (2013) [arXiv:1209.4803 [hep-ph]].
  * (115) J. Calle, D. Restrepo and Ó. Zapata, Phys. Rev. D 104, no.1, 1 (2021) [arXiv:2103.15328 [hep-ph]].