#  Large $(g-2)_{\mu}$ and signals of decays $e_{b}\rightarrow e_{a}\gamma$ in
a 3-3-1 model with inverse seesaw neutrinos

L. T. Hue<EMAIL_ADDRESS>Institute for Research and Development, Duy
Tan University, Da Nang City 50000, Vietnam  H. T. Hung<EMAIL_ADDRESS>Department of Physics, Hanoi Pedagogical University 2, Phuc Yen, Vinh Phuc
15000, Vietnam  N. T. Tham<EMAIL_ADDRESS>Department of Physics,
Hanoi Pedagogical University 2, Phuc Yen, Vinh Phuc 15000, Vietnam  H. N. Long
<EMAIL_ADDRESS>Institute of Physics, Vietnam Academy of Science and
Technology, 10 Dao Tan, Ba Dinh, 10000 Hanoi, Vietnam  T.Phong Nguyen 111
Corresponding author<EMAIL_ADDRESS>Department of Physics, Can Tho
University, 3/2 Street, Ninh Kieu, Can Tho City 94000, Vietnam

###### Abstract

We show that under current experimental bounds of the decays $e_{a}\rightarrow
e_{b}\gamma$, the recent experimental data of the muon anomalous magnetic
dipole moment $(g-2)_{\mu}$ can be explained in the framework of the 3-3-1
model with right-handed neutrinos. In addition, all of these branching ratios
can reach closely the recent experimental upper bounds.

###### pacs:

12.60.Fr, 13.15.+g, 14.60.St, 14.80.Bn

##  I  Introduction

At present, the experimental data on the anomalous dipole magnetic moments of
electron and muon $a_{e,\mu}=(g_{e,\mu}-2)/2$ show significant deviations from
their values predicted by the Standard Model (SM) Hagiwara:2011af ;
Davier:2017zfy ; Parker:2018vye ; Zyla:2020zbs . From the combination of
various different contributions 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 , the recent improved value of $a_{\mu}$
predicted by the SM is accepted widely as follows Aoyama:2020ynm :
$a^{\mathrm{SM}}_{\mu}=116591810(43)\times 10^{-11}$. The latest experimental
measurement has been reported from Fermi National Accelerator Laboratory
Abi:2021gix , $a^{\mathrm{exp}}_{\mu}=116592061(41)\times 10^{-11}$, leading
to the improved standard deviation of 4.2 $\sigma$ from the SM prediction,
namely

$\Delta a_{\mu}\equiv a^{\mathrm{exp}}_{\mu}-a^{\mathrm{SM}}_{\mu}=251\times
10^{-11}\pm 59\times 10^{-11}.$ (1)

On the other hand, the recent constraints on the charged lepton flavor
violating (cLFV) decays, $e_{b}\rightarrow e_{a}\gamma$ are TheMEG:2016wtm ;
Aubert:2009ag :

$\displaystyle\mathrm{Br}(\tau\rightarrow\mu\gamma)$ $\displaystyle<4.4\times
10^{-8},$ $\displaystyle\mathrm{Br}(\tau\rightarrow e\gamma)$
$\displaystyle<3.3\times 10^{-8},$ $\displaystyle\mathrm{Br}(\mu\rightarrow
e\gamma)$ $\displaystyle<4.2\times 10^{-13}.$ (2)

Many recent versions of the 3-3-1 models were indicated that they are
difficult to explain simultaneously all of these experimental constraints
Kelso:2014qka ; Binh:2015jfz ; Lindner:2016bgg ; deJesus:2020upp ;
deJesus:2020ngn ; Ky:2000ku  with the very large TeV values of the $SU(3)_{L}$
symmetry scale. Namely, the discussion on Ref. Kelso:2014qka  needs the cLFV
constraints from experimental data to rule out large $\Delta a_{\mu}$. The
remaining models rule out large $\Delta a_{\mu}$ for large $SU(3)_{L}$
symmetry scale with order of $\mathcal{O}(1)$ TeV, if no new $SU(3)_{L}$ Higgs
triplet or vectorlike charged lepton are added. This result can be explained
qualitatively from a consequence that a one-loop contribution from a heavy
gauge boson $V$ is different from that of the $W^{\pm}$ boson by a small
factor $m^{2}_{W}/m_{V}^{2}\geq 10^{-3}$. Similarly, one-loop contributions
from heavy Higgs boson $S$ have a suppressed factor $m^{2}_{h}/m^{2}_{S}$,
where $m_{h}$ is the mass of the standard model (SM-like) Higgs boson. In
addition, these Higgs contributions are constrained strictly by the small
upper bound of Br$(\mu\rightarrow e\gamma)$, leading to a strict constraint on
the doubly Higgs mass for the 3-3-1 models adding a $SU(3)_{L}$ Higgs sextet
to explain the experimental neutrino oscillation data. Adding new particles as
Higgs triplets or vectorlike charged leptons into the original 3-3-1 models to
generate new couplings contributing to $\Delta a_{\mu}$ is a popular way to
explain successful the experimental data of $a_{\mu}$ deJesus:2020upp ;
deJesus:2020ngn , but there seems irrelevant with neutrino oscillation data.
Some recent extensions of 3-3-1 models with discrete symmetries
CarcamoHernandez:2019lhv ; CarcamoHernandez:2020pxw  need a large number of
new leptons and Higgs bosons for the explanation of large $\Delta a_{\mu,e}$
consistent with experiments. On the other hand, a recent note indicated that a
version of the 3-3-1 model with right-handed neutrino (331RN) with heavy
neutral fermions assigned as $SU(3)_{L}$ gauge singlets (called the 331ISS
model for short) can predict large one-loop contributions from singly charged
Higgs bosons and inverse seesaw (ISS) neutrinos enough to explain the recent
$(g-2)_{\mu}$ data Dinh:2020pqn . More interesting, the model contains two
singly charged Higgs bosons, which may result in a special possibility that
two one-loop contributions to $\Delta a_{\mu}$ are large and constructive,
while those relate with cLFV decay amplitudes are strongly destructive. In
this work, we will pay attention to this possibility, namely we will try to
answer a question whether there exist any allowed regions of the parameter
space that the destructive properties of the Higgs contributions are enough to
satisfy the cLFV experimental constraints given in Eq. (I), and explain
successfully the recent data given in Eq. (1). We will use the 3-3-1 model
with the general Higgs potential given in Ref. Chang:2006aa ; Sanchez-
Vega:2016dwe . The 3-3-1 models explaining active neutrino data based on the
ISS mechanism has been discussed widely previously Catano:2012kw ; Dias:2012xp
; Pires:2018kaj ; Boucenna:2015zwa , but the interesting regions of the
parameter space allowing large $\Delta a_{\mu}$ data and consistent with
recent cLFV experimental constraints were not shown. In addition, the
Br$(\tau\rightarrow\mu\gamma,e\gamma)$ were predicted to be smaller than
Br$(\mu\rightarrow e\gamma)$, which is very suppressed with the recent and
upcoming experimental sensitivities of the order $\mathcal{O}(10^{-9})$
Baldini:2013ke ; Aushev:2010bq . Many other models beyond the SM with the ISS
mechanism can explain consistently the experimental data of $\Delta a_{\mu}$
and cLFV constraints Cao:2019evo ; Cao:2021lmj ; Nomura:2021adf ;
Mondal:2021vou . Here we analyze predictions of the 3-3-1 model with right-
handed neutrinos for the above observables.

Our work is arranged as follows. We will review the 331ISS model in Sec. II,
summarize the gauge, Higgs bosons and the lepton sectors. In Sec. III, we
introduce the analytic formulas to calculate the muon magnetic dipole moment
and the cLFV branching ratios. In Sec. IV, we discuss on the effect of a new
singly charged Higgs boson that can give one-loop contributions to $\Delta
a_{\mu}$ and cLFV amplitudes enough to explain successful all the experimental
data under consideration. In Sec. V, illustrations for numerical results are
given to indicate the existence of the allowed regions satisfying the
experimental data mentioned in this work. The conclusion is presented in the
last Sec. VI, where important results will be summarized.

##  II  Review the 3-3-1ISS model

###  II.1 Gauge bosons and fermions

The particle content of the 331ISS model was introduced in Refs.
Boucenna:2015zwa ; Nguyen:2018rlb  where active neutrino masses and
oscillations are originated from the ISS mechanism. The quark sector and
$SU(3)_{C}$ representations are irrelevant in this work, and hence they are
omitted here. We refer Ref. Boucenna:2015zwa  for a quark discussion. The
electric charge operator corresponding to the gauge group $SU(3)_{L}\times
U(1)_{X}$ is $Q=T_{3}-\frac{1}{\sqrt{3}}T_{8}+X$, where $T_{3,8}$ are the
diagonal $SU(3)_{L}$ generators. Each lepton family consists of a $SU(3)_{L}$
triplet $\psi_{aL}=(\nu_{a},~{}e_{a},N_{a})_{L}^{T}\sim(3,-\frac{1}{3})$ and a
right-handed charged lepton $e_{aR}\sim(1,-1)$ with $a=1,2,3$. Each left-
handed neutrino $N_{aL}=(N_{aR})^{c}$ is equivalent with a new right-handed
neutrinos defined in previous 331RN models Foot:1994ym . The only difference
between the two models 331RN and 331ISS is that, the 331ISS model contains
three more right-handed neutrinos transforming as gauge singlets,
$X_{aR}\sim(1,0)$, $a=1,2,3$. They couple with the $SU(3)_{L}$ Higgs triplets
to generate the neutrino mass term relating with the ISS mechanism. The three
Higgs triplets
$\rho=(\rho^{+}_{1},~{}\rho^{0},~{}\rho^{+}_{2})^{T}\sim(3,\frac{2}{3})$,
$\eta=(\eta_{1}^{0},~{}\eta^{-},\eta^{0}_{2})^{T}\sim(3,-\frac{1}{3})$, and
$\chi=(\chi_{1}^{0},~{}\chi^{-},\chi^{0}_{2})^{T}\sim(3,-\frac{1}{3})$ have
the following necessary vacuum expectation values for generating all tree-
level quark masses and leptons:
$\langle\rho\rangle=(0,\,\frac{v_{1}}{\sqrt{2}},\,0)^{T}$,
$\langle\eta\rangle=(\frac{v_{2}}{\sqrt{2}},\,0,\,0)^{T}$ and
$\langle\chi\rangle=(0,\,0,\,\frac{w}{\sqrt{2}})^{T}$.

The gauge bosons get masses through the covariant kinetic term of the Higgs
triplets,
$\mathcal{L}^{H}=\sum_{H=\chi,\eta,\rho}\left(D_{\mu}H\right)^{\dagger}\left(D^{\mu}H\right)$,
where the covariant derivative for the electroweak symmetry is
$D_{\mu}=\partial_{\mu}-ig{W}_{\mu}^{a}{T^{a}}-i{g_{X}}{T^{9}}X{X_{\mu}}$,
$a=1,2,..,8$. Note that $T^{9}\equiv\frac{I_{3}}{\sqrt{6}}$ and
$\frac{1}{\sqrt{6}}$ for (anti)triplets and singlets Buras:2012dp . Matching
with the SM gives $e=g\,s_{W}$ and
$\frac{g_{X}}{g}=\frac{3\sqrt{2}s_{W}}{\sqrt{3-4s^{2}_{W}}}$, where $e$ and
$s_{W}$ are respective the electric charge and sine of the Weinberg angle,
$s^{2}_{W}\simeq 0.231$. The relation $\frac{g_{X}}{g}$ is the same for both
choices of triplet or antitriplets representations of the left-handed leptons
Buras:2014yna ; Hue:2018dqf . The derivation of this relation is summarized as
follows. The 3-3-1 models have two spontaneous breaking steps:
$SU(3)_{L}\times U(1)_{X}\overset{w}{\longrightarrow}SU(2)_{L}\times
U(1)_{Y}\overset{v_{1},v_{2}}{\longrightarrow}U(1)_{Q}$. The first breaking
step with $w\neq 0$ generates masses for heavy particles predicted by the
$SU(3)_{L}$ symmetry. The neutral gauge bosons will change into the basis
containing the SM ones $W^{3}_{\mu}$ and $B_{\mu}$:
$(W^{3}_{\mu},\;W^{8}_{\mu},X_{\mu})\overset{w\neq
0,v_{1}=v_{2}=0}{\longrightarrow}(W^{3}_{\mu},Z^{\prime}_{\mu},B_{\mu})$.
Diagonalizing the squared mass matrix of these neutral gauge bosons will get a
massive eigenstate $Z^{\prime}$ with $m^{2}_{Z^{\prime}}\sim w^{2}$ and two SM
massless states $W^{3}_{\mu}$ and $B_{\mu}$. The relations between the two
bases before and after the first breaking step are $W^{8}_{\mu}=\frac{\beta
t}{\sqrt{6+\beta^{2}t^{2}}}B_{\mu}-\frac{\sqrt{6}}{\sqrt{6+\beta^{2}t^{2}}}Z^{\prime}_{\mu}$,
and $X_{\mu}=\frac{\sqrt{6}}{\sqrt{6+\beta^{2}t^{2}}}B_{\mu}+\frac{\beta
t}{\sqrt{6+\beta^{2}t^{2}}}Z^{\prime}_{\mu}$, with $t\equiv g_{X}/g$.
Inserting these relations to the covariant derivation of the $3-3-1$ gauge
group and keeping the part used to identify with the SM one, we have

$\displaystyle D^{3-3-1}_{\mu}\to
D^{\mathrm{SM}}_{\mu}=\partial_{\mu}-igT^{3}W^{3}_{\mu}-i\frac{gt}{\sqrt{6+\beta^{2}t^{2}}}\left(\beta
T^{8}+\sqrt{6}T^{9}X\right)B_{\mu},$

which results in the consequences that $g$ and
$\frac{gt}{\sqrt{6+\beta^{2}t^{2}}}=gt_{W}$ are the gauge couplings of the SM,
and the $U(1)_{Y}$ charge of the SM is $Y/2=\beta T^{8}+\mathbb{I}X$.

Like the 331RN model, the 331ISS model includes two pairs of singly charged
gauge bosons with the following physical states $W^{\pm}$ and $Y^{\pm}$ and
masses

$\displaystyle W^{\pm}_{\mu}$ $\displaystyle=\frac{W^{1}_{\mu}\mp
iW^{2}_{\mu}}{\sqrt{2}},\;Y^{\pm}_{\mu}=\frac{W^{6}_{\mu}\pm
iW^{7}_{\mu}}{\sqrt{2}},\;m_{W}^{2}=\frac{g^{2}}{4}\left(v_{1}^{2}+v_{2}^{2}\right),\;m_{Y}^{2}=\frac{g^{2}}{4}\left(w^{2}+v_{1}^{2}\right).$
(3)

The bosons $W^{\pm}$ are identified with the SM ones, leading to the
consequence that

$\displaystyle v_{1}^{2}+v_{2}^{2}\equiv v^{2}=(246\mathrm{GeV})^{2}.$ (4)

The general Higgs potential relating with the 331RN model will be applied in
our work with $v_{1}\neq v_{2}$. We will use the following parameters for this
general case.

$t_{\beta}\equiv\tan\beta=\frac{v_{2}}{v_{1}},\quad v_{1}=vc_{\beta},\quad
v_{2}=vs_{\beta}.$ (5)

The parameter $t_{\beta}$ plays a similar role known in the well-known models
with two Higgs doublet and the minimal supersymmetric Standard Model. This is
different from Ref. Nguyen:2018rlb , where $v_{1}=v_{2}$ was assumed so that
the Higgs potential given in Ref. Hue:2015fbb  was used to find the exact
physical state of the SM-like Higgs boson. This simple condition was also used
in previous discussions in 3-3-1 models addressed with anomalous magnetic
dipole moments deJesus:2020upp ; deJesus:2020ngn . As we will show below,
large $t_{\beta}\neq 1$ is one of the key condition for predicting large
$(g-2)_{\mu}$ consistent with experiments. The reason is that the physical
states of the charged Higgs bosons are determined analytically from this Higgs
potential, and only these Higgs bosons contribute significantly to one-loop
corrections to the $(g-2)_{\mu}$.

The Yukawa Lagrangian for generating lepton masses is:

$\displaystyle\mathcal{L}^{\mathrm{Y}}_{l}=-h^{e}_{ab}\overline{\psi_{aL}}\rho
e_{bR}+h^{\nu}_{ab}\epsilon^{ijk}\overline{(\psi_{aL})_{i}}(\psi_{bL})^{c}_{j}\rho^{*}_{k}-Y_{ab}\overline{\psi_{aL}}\,\chi
X_{bR}-\frac{1}{2}(\mu_{X})^{*}_{ba}\overline{(X_{aR})^{c}}X_{bR}+\mathrm{H.c.}.$
(6)

Here we assumed that the model under consideration respects a new lepton
number symmetry $\mathcal{L}$ discussed in Ref. Chang:2006aa  so that the term
$\overline{\psi_{aL}}\,\eta X_{bR}$ is not allowed in the above Yukawa
Lagrangian, while the soft-breaking term
$(\mu_{X})^{*}_{ba}\overline{(X_{aR})^{c}}X_{bR}$ is allowed with small
$(\mu_{X})_{ba}$. The new lepton number $\mathcal{L}$ called by generalized
lepton number CarcamoHernandez:2017cwi  is defined as
$L=\frac{4}{\sqrt{3}}T^{8}+\mathcal{L}\mathbb{I}$, where $L$ is the normal
lepton number. The specific assignment of $\mathcal{L}$ is
$\mathcal{L}(\rho)=-1/3$, $\mathcal{L}(\eta)=-2/3$, $\mathcal{L}(\chi)=4/3$,
$\mathcal{L}(\psi_{aL})=1/3$, which guarantees the consistence for the well-
known definition of $L$, namely $L(\ell)=1$ for $\ell=e_{aL,R},\nu_{aL}$,
$L(\ell)=-1$ for $\ell=N_{aL},X_{aR}$, and $L(q)=0$ for all SM quarks
Chang:2006aa .

The first term in Lagrangian (6) generates charged lepton masses
$m_{e_{a}}\equiv\frac{h^{e}_{ab}v_{1}}{\sqrt{2}}\delta_{ab}$, i.e, the mass
matrix of the charged leptons is assumed to be diagonal, hence the flavor
states of the charged leptons are also the physical ones. In the basis
$\nu^{\prime}_{L}=(\nu_{L},N_{L},(X_{R})^{c})^{T}$ and
$(\nu^{\prime}_{L})^{c}=((\nu_{L})^{c},(N_{L})^{c},X_{R})^{T}$ of the neutral
leptons, Lagrangian (6) gives a neutrino mass term corresponding to a block
form of the mass matrix Nguyen:2018rlb , namely

$\displaystyle-{\mathcal{L}}^{\nu}_{\mathrm{mass}}=\frac{1}{2}\overline{\nu^{\prime}_{L}}M^{\nu\dagger}(\nu^{\prime}_{L})^{c}+\mathrm{H.c.},\,\mathrm{where}\quad
M^{\nu\dagger}=\begin{pmatrix}0&m_{D}&0\\\ m_{D}^{T}&0&M_{R}\\\
0&M_{R}^{T}&\mu_{X}^{\dagger}\end{pmatrix},$ (7)

where $M_{R}$ is a $3\times 3$ matrix $(M_{R})_{ab}\equiv
Y_{ab}\frac{w}{\sqrt{2}}$, $(m_{D})_{ab}\equiv\sqrt{2}h^{\nu}_{ab}v_{1}$ with
$a,b=1,2,3$. Neutrino subbases are denoted as
$\nu_{R}=((\nu_{1L})^{c},(\nu_{2L})^{c},(\nu_{3L})^{c})^{T}$,
$N_{R}=((N_{1L})^{c},(N_{2L})^{c},(N_{3L})^{c})^{T}$, and
$X_{L}=((X_{1R})^{c},(X_{2R})^{c},(X_{3R})^{c})^{T}$. The mass matrix $M_{R}$
does not appear in the 331RN. The Dirac neutrino mass matrix $m_{D}$ must be
antisymmetric. The matrix $\mu_{X}$ defined in Eq. (6) is symmetric and it can
be diagonalized by a transformation $U_{X}$:

$U_{X}^{T}\mu_{X}U_{X}=\mathrm{diag}\left(\mu_{X,1},\mu_{X,2},\mu_{X,3}\right).$
(8)

The matrix $U_{X}$ will be absorbed by redefinition the states $X_{a}$,
therefore $\mu_{X}$ will be set as the diagonal matrix given in the right hand
side of Eq. (8).

The mass matrix $M^{\nu}$ is diagonalized by a $9\times 9$ unitary matrix
$U^{\nu}$,

$\displaystyle U^{\nu
T}M^{\nu}U^{\nu}=\hat{M}^{\nu}=\mathrm{diag}(m_{n_{1}},m_{n_{2}},...,m_{n_{9}})=\mathrm{diag}(\hat{m}_{\nu},\hat{M}_{N}),$
(9)

where $m_{n_{i}}$ ($i=1,2,...,9$) are masses of the nine physical neutrino
states $n_{iL}$. They consist of three active neutrinos $n_{aL}$ ($a=1,2,3$)
corresponding to the mass submatrix
$\hat{m}_{\nu}=\mathrm{diag}(m_{n_{1}},\;m_{n_{2}},\;m_{n_{3}})$, and the six
extra neutrinos $n_{IL}$ ($I=4,5,..,9$) with
$\hat{M}_{N}=\mathrm{diag}(m_{n_{4}},\;m_{n_{5}},...,\;m_{n_{9}})$. The ISS
mechanism leads to the following approximation solution of $U^{\nu}$,

$\displaystyle
U^{\nu}=\Omega\left(\begin{array}[]{cc}U_{\mathrm{PMNS}}&\mathbf{O}\\\
\mathbf{O}&V\\\
\end{array}\right),\;\;\Omega=\exp\left(\begin{array}[]{cc}\mathbf{O}&R\\\
-R^{\dagger}&\mathbf{O}\\\
\end{array}\right)=\left(\begin{array}[]{cc}1-\frac{1}{2}RR^{\dagger}&R\\\
-R^{\dagger}&1-\frac{1}{2}R^{\dagger}R\\\
\end{array}\right)+\mathcal{O}(R^{3}),$ (16)

where

$\displaystyle R^{*}$ $\displaystyle\simeq\left(-m^{*}_{D}M^{-1},\quad
m^{*}_{D}(M_{R}^{\dagger})^{-1}\right),\quad M\equiv
M^{*}_{R}\mu_{X}^{-1}M_{R}^{\dagger},$ (17) $\displaystyle
m^{*}_{D}M^{-1}m^{\dagger}_{D}$ $\displaystyle\simeq m_{\nu}\equiv
U^{*}_{\mathrm{PMNS}}\hat{m}_{\nu}U^{\dagger}_{\mathrm{PMNS}},$ (18)
$\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,\;M_{N}\equiv\begin{pmatrix}0&M_{R}^{*}\\\
M_{R}^{\dagger}&\mu_{X}\end{pmatrix}.$ (19)

The relations between the flavor and mass eigenstates are

$\nu^{\prime}_{L}=U^{\nu}n_{L},\quad\mathrm{and}\;(\nu^{\prime}_{L})^{c}=U^{\nu*}(n_{L})^{c},$
(20)

where $n_{L}\equiv(n_{1L},n_{2L},...,n_{9L})^{T}$ and
$(n_{L})^{c}\equiv((n_{1L})^{c},(n_{2L})^{c},...,(n_{9L})^{c})^{T}$. The
standard form of the lepton mixing matrix $U_{\mathrm{PMNS}}$ is the function
of three angles $\theta_{ij}$, one Dirac phase $\delta$ and two Majorana
phases $\alpha_{1}$, and $\alpha_{2}$ Tanabashi:2018oca , namely

$\displaystyle U^{\mathrm{PDG}}_{\mathrm{PMNS}}$
$\displaystyle=\begin{pmatrix}1&0&0\\\ 0&c_{23}&s_{23}\\\
0&-s_{23}&c_{23}\end{pmatrix}\,\begin{pmatrix}c_{13}&0&s_{13}e^{-i\delta}\\\
0&1&0\\\
-s_{13}e^{i\delta}&0&c_{13}\end{pmatrix}\,\begin{pmatrix}c_{12}&s_{12}&0\\\
-s_{12}&c_{12}&0\\\
0&0&1\end{pmatrix}\mathrm{diag}\left(1,e^{i\alpha_{1}},\,e^{i\alpha_{2}}\right)$
$\displaystyle=U^{0}_{\mathrm{PMNS}}\;\mathrm{diag}\left(1,e^{i\alpha_{1}},\,e^{i\alpha_{2}}\right),$
(21)

where $s_{ij}\equiv\sin\theta_{ij}$,
$c_{ij}\equiv\cos\theta_{ij}=\sqrt{1-s^{2}_{ij}}$, $i,j=1,2,3$ ($i<j$),
$0\leq\theta_{ij}<90\;[\mathrm{Deg.}]$ and $0<\delta\leq
720\;[\mathrm{Deg.}]$. The Majorana phases are chosen in the range
$-180\leq\alpha_{i}\leq 180$ [Deg.]

In this paper, we will work on the normal ordered scheme (NO) of the active
neutrino masses, which allows $\delta=\pi$ using in this work. The respective
best fit and the confidence level of 3$\sigma$ of the neutrino oscillation
experimental data is given as Zyla:2020zbs

$\displaystyle s^{2}_{12}=0.32,\;0.273\leq s^{2}_{12}\leq 0.379;$
$\displaystyle s^{2}_{23}=0.547,\;0.445\leq s^{2}_{23}\leq 0.599;$
$\displaystyle s^{2}_{13}=0.0216,\;0.0196\leq s^{2}_{13}\leq 0.0241;$
$\displaystyle\delta=218\;[\mathrm{Deg}],\;157\;[\mathrm{Deg}]\leq\delta\leq
349\;[\mathrm{Deg}];$ $\displaystyle\Delta m^{2}_{21}=7.55\times
10^{-5}[\mathrm{eV}^{2}],\quad 7.05\times 10^{-5}[\mathrm{eV}^{2}]\leq\Delta
m^{2}_{21}\leq 8.24\times 10^{-5}[\mathrm{eV}^{2}];$ $\displaystyle\Delta
m^{2}_{32}=2.424\times 10^{-3}[\mathrm{eV}^{2}],\quad 2.334\times
10^{-3}[\mathrm{eV}^{2}]\leq\Delta m^{2}_{32}\leq 2.524\times
10^{-3}[\mathrm{eV}^{2}].$ (22)

The above $CP$ phase is consistent with the updated one given in Ref.
Abe:2019vii , where the allowed range corresponding to $3\sigma$ confidence
level are $-3.41\leq\delta\leq-0.03$ ($164.6\leq\delta\leq 358.3$ [Deg.]) for
the NO scheme. The lepton mixing matrix defined in Eq. (18) relates with the
experimental parameters appearing in Eq.(II.1) are Tanabashi:2018oca

$\displaystyle s^{2}_{12}$
$\displaystyle=\frac{|\left(U_{\mathrm{PMNS}}\right)_{12}|^{2}}{1-|\left(U_{\mathrm{PMNS}}\right)_{13}|^{2}},\;s^{2}_{13}=|\left(U_{\mathrm{PMNS}}\right)_{13}|^{2},\;s^{2}_{23}=\frac{|\left(U_{\mathrm{PMNS}}\right)_{23}|^{2}}{1-|\left(U_{\mathrm{PMNS}}\right)_{13}|^{2}}.$
(23)

Additionally, it is easily to derive that

$\displaystyle e^{i\delta}$
$\displaystyle=\frac{c_{23}(c^{2}_{12}+ys^{2}_{12})}{s_{13}s_{23}s_{12}c_{12}(1-y)},\;y=\frac{\left(U_{\mathrm{PMNS}}\right)_{22}\left(U_{\mathrm{PMNS}}\right)_{11}}{\left(U_{\mathrm{PMNS}}\right)_{12}\left(U_{\mathrm{PMNS}}\right)_{21}},$
$\displaystyle e^{i\alpha_{1}}$
$\displaystyle=\frac{\left(U_{\mathrm{PMNS}}\right)_{12}c_{12}}{|\left(U_{\mathrm{PMNS}}\right)_{11}|s_{12}},\;e^{i\left(\alpha_{2}-\delta\right)}=\frac{\left(U_{\mathrm{PMNS}}\right)_{13}c_{13}c_{12}}{|\left(U_{\mathrm{PMNS}}\right)_{11}|s_{13}}.$
(24)

The detailed calculation shown in Ref. Nguyen:2018rlb , using the ISS
relations, yields

$m_{D}=zc_{\beta}\times\tilde{m}_{D},\;\tilde{m}_{D}=\begin{pmatrix}0&x_{12}&x_{13}\\\
-x_{12}&0&1\\\ -x_{13}&-1&0\end{pmatrix},$ (25)

where $z=\sqrt{2}v\,h^{\nu}_{23}$ is assumed to be positive and real,

$\displaystyle x^{*}_{12}$
$\displaystyle=\frac{(m_{\nu})_{11}(m_{\nu})_{23}-(m_{\nu})_{13}(m_{\nu})_{12}}{(m_{\nu})_{12}(m_{\nu})_{33}-(m_{\nu})_{13}(m_{\nu})_{23}},\;x^{*}_{13}=\frac{(m_{\nu})_{11}(m_{\nu})_{33}-(m_{\nu})^{2}_{13}}{(m_{\nu})_{12}(m_{\nu})_{33}-(m_{\nu})_{13}(m_{\nu})_{23}}.$
(26)

We note that the lightest active neutrino mass is zero at the tree level, but
can be nonzero when loop-corrections are included Chang:2006aa . Also, the
quantum effects can be considered for the charged lepton masses, so that the
regions predicting large $\Delta a_{\mu}$ may be larger Yin:2021yqy ;
Baker:2021yli  than the ones discussed in this work. The perturbative limit
requires that $h^{\nu}_{23}<\sqrt{4\pi}$, leading to the following upper bound
of $z$,

$z<1233\;[\mathrm{GeV}].$ (27)

The two formulas in Eq. (26) were found in the general symmetric from of
$M^{-1}$, namely they are be found by using Eq. (18) for off-diagonal entries
of $m_{\nu}$ to determine $\left(M^{-1}\right)_{ij}$, then insert them into
the diagonal ones. The off-diagonal elements of $M^{-1}$ are determined as
follows:

$\displaystyle\left(M^{-1}\right)_{12}$
$\displaystyle=\frac{1}{2}\left[x^{*}_{13}\left(M^{-1}\right)_{11}-\frac{\left(M^{-1}\right)_{22}}{x^{*}_{13}}-\frac{\left(m_{\nu}\right)_{13}+x^{*}_{13}\left(m_{\nu}\right)_{23}}{x^{*}_{12}x^{*}_{13}z^{2}}\right],$
$\displaystyle\left(M^{-1}\right)_{13}$
$\displaystyle=\frac{1}{2}\left[x^{*}_{12}\left(M^{-1}\right)_{11}+\frac{\left(M^{-1}\right)_{33}}{x^{*}_{12}}-\frac{\left(m_{\nu}\right)_{12}+x^{*}_{12}\left(m_{\nu}\right)_{23}}{x^{*}_{12}x^{*}_{13}z^{2}}\right],$
$\displaystyle\left(M^{-1}\right)_{23}$
$\displaystyle=\frac{1}{2}\left[\frac{x^{*2}_{12}\left(M^{-1}\right)_{22}-x^{*2}_{13}\left(M^{-1}\right)_{33}}{x^{*}_{12}x^{*}_{13}}+\frac{x^{*}_{13}\left(m_{\nu}\right)_{12}-x^{*}_{12}\left(m_{\nu}\right)_{13}}{x^{*}_{12}x^{*}_{13}z^{2}}\right].$
(28)

Hence all elements of the matrix $M^{-1}$ depend on only three complex
parameters $\left(M^{-1}\right)_{ii}$ with $i=1,2,3$. When identifying with
$M^{-1}=\left(M_{R}^{\dagger}\right)^{-1}\mu_{X}\left(M_{R}^{*}\right)^{-1}$
given in Eq. (17), six parameters $\mu_{X,i}$ and $\left(M^{-1}\right)_{ii}$
are determined as functions of elements of $M_{R}$. In this work, we will
consider all elements of $M_{R}$ are free parameters, namely

$\left(M_{R}\right)_{ij}=zc_{\beta}\times\left(\widetilde{M}_{R}\right)_{ij},\;\left(\widetilde{M}_{R}\right)_{ij}\equiv
k_{ij},$ (29)

where all $k_{ij}$ are assumed to be real for simplicity. The ISS relations
are valid with at least some $|k_{ij}|\gg 1$ and det$M\neq 0$. In the
numerical investigation, $m_{\nu}$ is determined from the $3\sigma$ neutrino
oscillation data through Eq. (18). The Dirac matrix $m_{D}$ is then determined
by Eq. (25). The free parameters $k_{ij}$ and $z$ are assumed to be real, and
$z$ is positive. The three elements of the matrix $\mu_{X}$ are determined as
functions of these free parameters. The respective formulas are lengthy hence
they are not written down explicitly here. In our work, we only consider the
case max$|\mu_{X,i}|\ll z$ hence all $(\mu_{X})_{i}$ gives suppressed mixing
elements in the total lepton mixing matrix $U^{\nu}$. This condition will
always be checked numerically to derive the final results.

In the numerical investigation, the free parameters $z$ and $k_{ij}$ will be
scanned in the valid ranges to construct the total neutrino mass matrix
defined in Eq. (7). After that, the mass eigenstates and the total mixing
matrix are calculated numerically with at least 30 digits of precision. Using
the relations listed in Eqs. (23) and (II.1), we reproduce all of the
oscillation parameters $\Delta m^{2}_{ij}$ and $s^{2}_{ij}$ then force them
satisfying the $3\sigma$ allowed data. This will help us to collect the
allowed values of $z$ and $k_{ij}$ in evaluating the cLFV branching ratios and
$(g-2)_{\mu}$ data. We emphasize that the regions of the parameter space in
our numerical investigation are more general than those mentioned in Refs.
Nguyen:2018rlb ; Dinh:2020pqn .

The Lagrangian for quark masses was discussed previously Chang:2006aa . Here,
we just remind the reader that the Yukawa couplings of the top quark must
satisfy the perturbative limit $h^{u}_{33}<\sqrt{4\pi}$, leading to a lower
bound $v_{2}>\frac{\sqrt{2}m_{t}}{\sqrt{4\pi}}$. Combining this with the
relations in Eqs. (4) and (5) gives a lower bound $t_{\beta}\geq 0.3$, which
will be used in the numerical discussion.

###  II.2 Higgs bosons

The Higgs potential used here respect the new lepton number defined in Ref.
Chang:2006aa , namely

$\displaystyle V_{h}$
$\displaystyle=\sum_{S=\eta,\rho,\chi}\left[\mu_{S}^{2}S^{\dagger}S+\lambda_{S}\left(S^{\dagger}S\right)^{2}\right]+\lambda_{12}(\eta^{\dagger}\eta)(\rho^{\dagger}\rho)+\lambda_{13}(\eta^{\dagger}\eta)(\chi^{\dagger}\chi)+\lambda_{23}(\rho^{\dagger}\rho)(\chi^{\dagger}\chi)$
$\displaystyle+\tilde{\lambda}_{12}(\eta^{\dagger}\rho)(\rho^{\dagger}\eta)+\tilde{\lambda}_{13}(\eta^{\dagger}\chi)(\chi^{\dagger}\eta)+\tilde{\lambda}_{23}(\rho^{\dagger}\chi)(\chi^{\dagger}\rho)+\sqrt{2}\omega
f\left(\epsilon_{ijk}\eta^{i}\rho^{j}\chi^{k}+\mathrm{h.c.}\right),$ (30)

where $f$ is a dimensionless parameter, which $f\omega$ is the same as that
used in previous works. The minimum conditions of the Higgs potential as well
as the identification of the SM-like Higgs were discussed in detailed
previously Ninh:2005su ; Hue:2015fbb . The model always contains a light $CP$
even neutral Higgs boson identified with the SM-like Higgs boson confirmed
experimentally. This Higgs boson gives suppressed contributions to
$(g-2)_{\mu}$ hence we will ignore it from now on. The model contains two
pairs of singly charged Higgs bosons $H^{\pm}_{1,2}$ and Goldstone bosons of
the gauge bosons $W^{\pm}$ and $Y^{\pm}$, which are denoted as $G^{\pm}_{W}$
and $G^{\pm}_{Y}$, respectively. The masses of all charged Higgs bosons are
Buras:2012dp ; Hue:2017lak ; Ninh:2005su
$m^{2}_{H^{\pm}_{1}}=\left(\frac{\tilde{\lambda}_{12}v^{2}}{2}+\frac{fw^{2}}{s_{\beta}c_{\beta}}\right)$,
$m^{2}_{H^{\pm}_{2}}=(v^{2}c_{\beta}^{2}+w^{2})\left(\frac{\tilde{\lambda}_{23}}{2}+ft_{\beta}\right)$,
and $m^{2}_{G^{\pm}_{W}}=m^{2}_{G^{\pm}_{Y}}=0$. The relations between the
original and mass eigenstates of the charged Higgs bosons are Ninh:2005su

$\displaystyle\left(\begin{array}[]{c}\eta^{\pm}\\\
\rho_{1}^{\pm}\end{array}\right)=\left(\begin{array}[]{cc}-s_{\beta}&c_{\beta}\\\
c_{\beta}&s_{\beta}\\\
\end{array}\right)\left(\begin{array}[]{c}G_{W}^{\pm}\\\
H_{1}^{\pm}\end{array}\right),\quad\left(\begin{array}[]{c}\rho_{2}^{\pm}\\\
\chi^{\pm}\end{array}\right)=\left(\begin{array}[]{cc}-s_{\theta}&c_{\theta}\\\
c_{\theta}&s_{\theta}\end{array}\right)\left(\begin{array}[]{c}G_{Y}^{\pm}\\\
H_{2}^{\pm}\end{array}\right),$ (43)

where $t_{\theta}=v_{1}/w$.

The model contains five $CP$-odd neutral scalar components. Three of them are
Goldstone bosons of the neutral gauge bosons $Z,Z^{\prime}$ and $X^{0}$. The
two remaining are physical states with masses
$m^{2}_{a_{1}}=f\left(c_{\beta}s_{\beta}v^{2}+\frac{\omega^{2}}{t_{\beta}}\right)+\frac{\bar{\lambda}_{13}}{2}\left(s_{\beta}{}^{2}v^{2}+\omega^{2}\right),\quad
m^{2}_{a_{2}}=f\left(\frac{\omega^{2}}{c_{\beta}s_{\beta}}+c_{\beta}s_{\beta}v^{2}\right).$
As a consequence, the parameter $f$ must be positive. In addition, $f$ may be
small so that charged Higgs boson masses can be around 1 TeV.

##  III Analytic formulas for one loop contributions to $\Delta a_{e_{a}}$ and
cLFV decays $e_{b}\rightarrow e_{a}\gamma$

All detailed steps for calculation to derive the couplings that give large
one-loop contributions were presented in Ref. Nguyen:2018rlb . We just collect
the final results related with this work. The condition $m_{e_{b}}>m_{e_{a}}$
is always used to define the one loop form factors $c^{X}_{(ab)R}$ and
$c^{X}_{(ba)R}$ introduced in Ref. Crivellin:2018qmi , which are different
from our notations by a relative factor $m_{e_{b}}$.

The relevant Lagrangian of charged gauge bosons is

$\displaystyle\mathcal{L}^{\ell
nV}=\overline{\psi_{aL}}\gamma^{\mu}D_{\mu}\psi_{aL}$
$\displaystyle\supset\frac{g}{\sqrt{2}}\sum_{i=1}^{9}\sum_{a=1}^{3}\left[U^{\nu*}_{ai}\overline{n_{i}}\gamma^{\mu}P_{L}e_{a}W^{+}_{\mu}+U^{\nu*}_{(a+3)i}\overline{n_{i}}\gamma^{\mu}P_{L}e_{a}Y^{+}_{\mu}\right],$
(44)

corresponding to the following one-loop form factors:

$\displaystyle c^{W}_{(ab)R}$
$\displaystyle=\frac{eg^{2}}{32\pi^{2}m_{W}^{2}}\sum_{i=1}^{9}U^{\nu}_{ai}U^{\nu*}_{bi}F_{LVV}\left(\frac{m^{2}_{n_{i}}}{m^{2}_{W}}\right),$
$\displaystyle c^{W}_{(ba)R}$
$\displaystyle=\frac{eg^{2}m_{e_{a}}}{32\pi^{2}m_{W}^{2}m_{e_{b}}}\sum_{i=1}^{9}U^{\nu}_{bi}U^{\nu*}_{ai}F_{LVV}\left(\frac{m^{2}_{n_{i}}}{m^{2}_{W}}\right),$
$\displaystyle c^{Y}_{(ab)R}$
$\displaystyle=\frac{eg^{2}}{32\pi^{2}m_{W}^{2}}\sum_{i=1}^{9}U^{\nu}_{(a+3)i}U^{\nu*}_{(b+3)i}\frac{m^{2}_{W}}{m^{2}_{Y}}\times
F_{LVV}\left(\frac{m^{2}_{n_{i}}}{m^{2}_{Y}}\right),$ $\displaystyle
c^{Y}_{(ba)R}$
$\displaystyle=\frac{eg^{2}m_{e_{a}}}{32\pi^{2}m_{W}^{2}m_{e_{b}}}\sum_{i=1}^{9}U^{\nu}_{(b+3)i}U^{\nu*}_{(a+3)i}\frac{m^{2}_{W}}{m^{2}_{Y}}\times
F_{LVV}\left(\frac{m^{2}_{n_{i}}}{m^{2}_{Y}}\right),$ (45)

where

$F_{LVV}(x)=-\frac{10-43x+78x^{2}-49x^{3}+4x^{4}+18x^{3}\ln(x)}{24(x-1)^{4}},$
(46)

$e=\sqrt{4\pi\alpha_{\mathrm{em}}}$ being the electromagnetic coupling
constant, and $g=e/s_{W}$.

Lagrangian of charged Higgs bosons is

$\displaystyle\mathcal{L}^{\ell
nH}=-\frac{g}{\sqrt{2}m_{W}}\sum_{k=1}^{2}\sum_{a=1}^{3}\sum_{i=1}^{9}H_{k}^{+}\overline{n_{i}}\left(\lambda^{L,k}_{ai}P_{L}+\lambda^{R,k}_{ai}P_{R}\right)e_{a}+H.c.,$
(47)

where

$\displaystyle\lambda^{R,1}_{ai}$
$\displaystyle=m_{e_{a}}U^{\nu*}_{ai}t_{\beta},\quad\lambda^{R,2}_{ai}=\frac{m_{e_{a}}c_{\theta}U^{\nu*}_{(a+3)i}}{c_{\beta}},$
$\displaystyle\lambda^{L,1}_{ai}$
$\displaystyle=-t_{\beta}\sum_{c=1}^{3}(m_{D}^{*})_{ac}U^{\nu}_{(c+3)i}=-s_{\beta}z\sum_{c=1}^{3}(\tilde{m}_{D}^{*})_{ac}U^{\nu}_{(c+3)i},$
$\displaystyle\lambda^{L,2}_{ai}$
$\displaystyle=\sum_{c=1}^{3}\frac{c_{\theta}}{c_{\beta}}\times\left[(m_{D}^{*})_{ac}U^{\nu}_{ci}+t^{2}_{\theta}(M_{R}^{*})_{ac}U^{\nu}_{(c+6)i}\right]$
$\displaystyle=c_{\theta}z\sum_{c=1}^{3}\left[(\tilde{m}_{D}^{*})_{ac}U^{\nu}_{ci}+t^{2}_{\theta}(\widetilde{M}_{R}^{*})_{ac}U^{\nu}_{(c+6)i}\right].$
(48)

The one-loop form factors are:

$\displaystyle c^{H,k}_{(ab)R}$
$\displaystyle=\frac{eg^{2}}{32\pi^{2}m^{2}_{W}m_{e_{b}}m^{2}_{H_{k}}}\sum_{i=1}^{9}\left[\lambda^{L,k*}_{ai}\lambda^{R,k}_{bi}m_{n_{i}}F_{LHH}\left(\frac{m^{2}_{n_{i}}}{m^{2}_{H_{k}}}\right)\right.$
$\displaystyle\quad\left.+\left(m_{e_{b}}\lambda^{L,k*}_{ai}\lambda^{L,k}_{bi}+m_{e_{a}}\lambda^{R,k*}_{ai}\lambda^{R,k}_{bi}\right)\tilde{F}_{LHH}\left(\frac{m^{2}_{n_{i}}}{m^{2}_{H_{k}}}\right)\right],$
$\displaystyle c^{H,k}_{(ba)R}$
$\displaystyle=\frac{eg^{2}}{32\pi^{2}m^{2}_{W}m_{e_{b}}m^{2}_{H_{k}}}\sum_{i=1}^{9}\left[\lambda^{L,k*}_{bi}\lambda^{R,k}_{ai}m_{n_{i}}F_{LHH}\left(\frac{m^{2}_{n_{i}}}{m^{2}_{H_{k}}}\right)\right.$
$\displaystyle\quad\left.+\left(m_{e_{a}}\lambda^{L,k*}_{bi}\lambda^{L,k}_{ai}+m_{e_{b}}\lambda^{R,k*}_{bi}\lambda^{R,k}_{ai}\right)\tilde{F}_{LHH}\left(\frac{m^{2}_{n_{i}}}{m^{2}_{H_{k}}}\right)\right],$
(49)

where $b\geq a$, and

$F_{LHH}(x)=-\frac{1-x^{2}+2x\ln(x)}{4(x-1)^{3}},\quad\widetilde{F}_{LHH}(x)=-\frac{-1+6x-3x^{2}-2x^{3}+6x^{2}\ln(x)}{24(x-1)^{4}}.$
(50)

The total one-loop contribution to the cLFV and $\Delta
a^{331\mathrm{ISS}}_{\mu}$ is

$\displaystyle c_{(ab)R}$
$\displaystyle=c^{W}_{(ab)R}+c^{Y}_{(ab)R}+c^{H_{1}}_{(ab)R}+c^{H_{2}}_{(ab)R},$
$\displaystyle c_{(ba)R}$
$\displaystyle=\left\\{c_{(ab)R}\left[a\leftrightarrow
b\right]\right\\}\times\frac{m_{e_{a}}}{m_{e_{b}}}.$ (51)

The one-loop contributions from charged gauge bosons to the $a_{e_{a}}$ and
the electric dipole moment $d_{e_{a}}$ of the charged lepton $e_{a}$ are
Crivellin:2018qmi :

$\displaystyle a^{V}_{e_{a}}$
$\displaystyle=a^{W}_{e_{a}}+a^{Y}_{e_{a}}\equiv-\frac{4m^{2}_{e_{a}}}{e}\left(\mathrm{Re}[c^{W}_{(aa)R}]+\mathrm{Re}[c^{Y}_{(aa)R}]\right),$
$\displaystyle d^{V}_{e_{a}}$
$\displaystyle=d^{W}_{e_{a}}+d^{Y}_{e_{a}}\equiv-2m_{e_{a}}\left(\mathrm{Im}[c^{W}_{(aa)R}]+\mathrm{Im}[c^{Y}_{(aa)R}]\right),$
(52)

The one-loop contribution to $a_{e_{a}}$ and $d_{e_{a}}$ caused by charged
Higgs bosons is Crivellin:2018qmi :

$\displaystyle a^{H}_{e_{a}}$
$\displaystyle=\sum_{k=1}^{2}a^{H,k}_{e_{a}},\;a^{H,k}_{e_{a}}\equiv-\frac{4m^{2}_{e_{a}}}{e}\mathrm{Re}[c^{H,k}_{(aa)R}],$
$\displaystyle d^{H}_{e_{a}}$
$\displaystyle=\sum_{k=1}^{2}d^{H,k}_{e_{a}},\;d^{H,k}_{e_{a}}\equiv-2m_{e_{a}}\mathrm{Im}[c^{H,k}_{(aa)R}].$
(53)

The quantity $\Delta d_{e_{a}}=d^{V}_{e_{a}}+d^{H}_{e_{a}}$ is the new one
loop contributions predicted to the electric dipole moment of the charged
leptons. It equals to zero when our investigation is limited in the case of
the Dirac phase $\delta=\pi$. This zero value of $d_{\mu}$ satisfies the
current experimental constraint Muong-2:2008ebm  hence we will not consider
from now on.

We remind the reader that one loop contributions from neutral Higgs bosons are
very suppressed hence they are ignored here. The reason is that the 331ISS
model has no new charged leptons, hence the one-loop contributions of any
neutral Higgs bosons $H^{0}$ to $c_{(ab)R}$ must arise only from the couplings
$H^{0}\bar{e}_{a}e_{a}$ derived from the first term of the Yukawa Laragian
(6). These couplings have the same Yukawa couplings with the SM-like
$h\sim\mathrm{Re}[\rho^{0}]/\sqrt{2}$, but different mixing factors
$|c_{H^{0}}|\leq 1$ telling the contributions of $\rho^{0}$ to the physical
state $H^{0}$. Hence these contributions to $a_{\mu}$ have the same form with
the one from the SM-like Higgs boson having mass $m_{h}\simeq 125$ GeV $\gg
m_{\mu}$,
$a^{h}_{\mu}\simeq\frac{\sqrt{2}G_{\mu}m^{2}_{\mu}}{4\pi^{2}}\times\frac{m^{2}_{\mu}}{m_{h}^{2}}\ln\frac{m_{h}^{2}}{m^{2}_{\mu}}\leq\mathcal{O}(10^{-14})$
Jegerlehner:2009ry . Also, the heavy neutral Higgs will give suppressed one-
loop contributions to $\Delta a_{\mu}$. The deviation of $a_{\mu}$ between
predictions by the two models 331ISS and SM are

$\displaystyle\Delta a^{\mathrm{331ISS}}_{e_{a}}$ $\displaystyle\equiv\Delta
a_{e_{a}}=\Delta
a^{W}_{e_{a}}+a^{Y}_{e_{a}}+a^{H,1}_{e_{a}}+a^{H,2}_{e_{a}},\quad\Delta
a^{W}_{e_{a}}=a^{W}_{e_{a}}-a^{\mathrm{SM},W}_{e_{a}},$ (54)

where $a^{\mathrm{SM},W}_{\mu}=3.887\times 10^{-9}$ Jegerlehner:2009ry  is the
SM prediction for the one-loop contribution from $W$ boson
$c^{W,\mathrm{SM}}_{(22)R}$. In this work, $\Delta
a^{\mathrm{331ISS}}_{\mu}=\Delta a_{\mu}$ will be considered as new physics
predicted by the 331ISS and will be used to compare with the experimental data
in the following numerical investigation. We note that the discrepancy of
$a_{e}$ between experiments and SM is about 2.5 standard deviation
Parker:2018vye ; Aoyama:2012wj ; Aoyama:2012wk ; Laporta:2017okg ;
Terazawa:2018pdc ; Volkov:2019phy . In this work we will only pay attention to
the $\Delta a_{\mu}$ which is the very interesting result of 4.2 standard
deviation and may be a clear signal of new physics in the near future.

Based on Ref. Crivellin:2018qmi , the branching ratios of the cLFV processes
are

$\mathrm{Br}(e_{b}\rightarrow
e_{a}\gamma)\simeq\frac{48\pi^{2}}{G_{F}^{2}}\left(\left|c_{(ab)R}\right|^{2}+\left|c_{(ba)R}\right|^{2}\right)\mathrm{Br}(e_{b}\rightarrow
e_{a}\overline{\nu_{a}}\nu_{b}),$ (55)

where $G_{F}=g^{2}/(4\sqrt{2}m_{W}^{2})$. This result is consistent with the
formulas given used in Refs. Hue:2017lak ; Nguyen:2018rlb  for 3-3-1 models.

It is noted that for the gauge boson contributions, we have
$|c^{V}_{(ba)R}|/|c^{V}_{(ab)R}|=m_{e_{a}}/m_{e_{b}}\ll 1$ for
$m_{e_{b}}>m_{e_{a}}$. Similarly, we can estimate that
$|c^{H,k}_{(ba)R}|/|c^{H,k}_{(ab)R}|\ll 1$ for every particular contribution.
Anyways, in the general case we cannot ignore $c^{X}_{(ba)R}$ because of the
situation that when contributions to $c_{(ab)R}$ have the same order but some
of them have opposite signs. Then the very destructive correlations among
particular Higgs contributions in the $c_{(ab)R}$ will result in the same
order of both $|c_{(ab)R}|$ and $|c_{(ba)R}|$. This will happen in the 331ISS
model when $\Delta a^{331\mathrm{ISS}}_{\mu}=\mathcal{O}(10^{-9})$
corresponding to the order of the experimental data and Br$(\mu\rightarrow
e\gamma)<4.2\times 10^{-13}$ require both conditions of
$\mathcal{O}(10^{-9})\;[\mathrm{GeV}^{-2}]\leq|c_{(22)R}|\leq\mathcal{O}(10^{-8})\;[\mathrm{GeV}^{-2}]$
and $|c_{(21)R}|\leq\mathcal{O}(10^{-13})\;[\mathrm{GeV}^{-2}]$, respectively.
As a result, we can estimate that the one-loop contributions from two charged
Higgs bosons to Br($\mu\rightarrow e\gamma$) are strongly destructive, i.e.
$c^{H_{1}}_{(12)}\simeq-c^{H_{2}}_{(12)}$. Simultaneously,
$|c^{H_{k}}_{(12)}|\sim|c^{H_{k}}_{(22)}|$, therefore the charged Higgs
contributions to $\Delta a_{\mu}$ must be constructive and satisfy
$|c^{H_{1}}_{(22)}|\sim|c^{H_{2}}_{(22)}|\sim\mathcal{O}(10^{-9})-\mathcal{O}(10^{-8})\;[\mathrm{GeV}^{-2}]$,
or they can be destructive but $|c^{H_{i}}_{(22)}|\gg|c^{H_{j}}_{(22)}|$ with
$i\neq j$. These important properties of charged Higgs boson contributions
will be the key point in our numerical investigation to collect data points
satisfying the large values of $\Delta a_{\mu}\geq 10^{-9}$ before considering
any cLFV decay constraints. The gauge contributions are suppressed hence we do
not discuss qualitatively here, but they are also included in the numerical
investigation. We just pay attention to the two key one-loop charged Higgs
boson contributions which will affect two other cLFV decays $\tau\rightarrow
e\gamma,\mu\gamma$.

The experimental constraints of the form factors $c_{(ab)R}$ are listed in
Table 1, where the allowed values of $\Delta a_{\mu}$ are chosen in the range
of $1\sigma$ confidence level given in Eq. (1).

$192<\Delta a_{\mu}\times 10^{11}<310,$ | $-4.8\times 10^{-8}\;[\mathrm{GeV}^{-2}]<c_{(22)R}<-3.99\times 10^{-8}\;[\mathrm{GeV}^{-2}]$
---|---
Br$(\mu\rightarrow e\gamma)$ | $|c_{(21)R}|,\;|c_{(12)R}|<3.47\times 10^{-13}\;[\mathrm{GeV}^{-2}]$
Br$(\tau\rightarrow e\gamma)$ | $|c_{(31)R}|,\;|c_{(13)R}|<2.31\times 10^{-10}\;[\mathrm{GeV}^{-2}]$
Br$(\tau\rightarrow\mu\gamma)$ | $|c_{(32)R}|,\;|c_{(23)R}|<2.63\times 10^{-10}\;[\mathrm{GeV}^{-2}]$
Table 1: Constraints of $c_{(ab)R}\;[\mathrm{GeV}^{-2}]$ from experimental
data. The allowed values of $\Delta a_{\mu}$ satisfying a confidence level of
$1\sigma$ from the experimental data given in Eq. (1).

We derive that the allowed regions of the parameter space have the following
properties:

$\displaystyle\left|\frac{c_{(12)R}}{c_{(22)R}}\right|,\left|\frac{c_{(21)R}}{c_{(22)R}}\right|\leq\mathcal{O}(10^{-5});\;\left|\frac{c_{(13)R}}{c_{(22)R}}\right|,\;\left|\frac{c_{(31)R}}{c_{(22)R}}\right|,\;\left|\frac{c_{(23)R}}{c_{(22)R}}\right|,\;\left|\frac{c_{(32)R}}{c_{(22)R}}\right|\leq\mathcal{O}(10^{-2}).$
(56)

Normally, our numerical scan gives a relation that
$|c^{H_{k}}_{(22)R}|/|c^{H_{k}}_{(ab)R}|\leq\mathcal{O}(10)$ with $a\neq b$.
As a result, the huge destructive correlation between charged Higgs
contributions to guarantee simultaneously the experimental constraints of
Br$(\mu\rightarrow e\gamma)$ and $\Delta a_{\mu}$. Also, the two cLFV decays
of $\tau\rightarrow e\gamma,\mu\gamma$ also need smaller but still large
destructive charged Higgs contributions to satisfy the upper experimental
bounds because some of these particular contributions often satisfy
$|c^{H_{k}}_{(13)R}|/|c^{H_{k}}_{(22)R}|,|c^{H_{k}}_{(23)R}|/|c^{H_{k}}_{(22)R}|\geq
0.1$. While $|c^{H_{k}}_{(31)R}|,|c^{H_{k}}_{(32)R}|\ll
10^{-10}\;[\mathrm{GeV}^{-2}]$, consequently they are subdominant to the cLFV
decays where their branching ratios are close to the upper experimental
constraints. The mentioned properties are very important for us to point out
the validation of the allowed regions.

For convenience in estimating qualitatively the above properties, we define
new important quantities determining the correlations between two charged
Higgs contributions in a physical process as follows:

$\displaystyle R^{X}_{ab}$
$\displaystyle\equiv\left|\frac{\mathrm{Re}[c^{X}_{(ab)R}]}{\mathrm{Re}[c_{(ab)R}]}\right|;\;a,b=1,2,3;\;X=W,Y,H^{\pm}_{1},H^{\pm}_{2},$
(57) $\displaystyle R^{-}_{ab}$
$\displaystyle\equiv\left|\frac{\mathrm{Re}[c^{H_{1}}_{(ab)R}+c^{H_{2}}_{(ab)R}]}{\mathrm{Re}[c_{(ab)R}]}\right|.$
(58)

The first ratio $R^{X}_{ab}$ shows the relative contribution from the particle
$X$ in the loop to the total contribution. The second one shows the relative
contributions of both singly charged Higgs bosons. In the 331ISS model, we
will see that the relations $R^{W}_{ab},R^{Y}_{ab}\ll R^{H_{k}}_{ab}$ often
happens. The interesting possibility we would like to discuss is that large
contributions of $H_{k}$ for large $\Delta a_{\mu}$ with $|R^{H_{k}}_{22}|\sim
10^{-1}$, while the huge destructive correlations of these two Higgs bosons
$\mathrm{Re}[c^{H_{1}}_{(ab)R}/c^{H_{2}}_{(ab)R}]\simeq-1$ will allow small
cLFV constraints. Quantitatively, we estimate that $R^{H_{k}}_{ab}\gg 1$ and
$R^{-}_{ab}\simeq 1$, with $a\neq b$. The details of numerical investigation
will be shown below.

##  IV  Additional singly charged Higgs boson for an explanation of
$(g-2)_{\mu}$ data at 1 $\sigma$ deviation

The appearance of the gauge singlet $X_{R}$ leads to a possibility that, a new
singly charged Higgs bosons $h_{3}^{\pm}\sim(1,1,\pm 1)$ can be included in
the 331ISS model so that they can give one-loop contributions to both $\Delta
a_{\mu}$ and cLFV amplitudes through the following Yukawa interactions:

$\mathcal{L}^{Y}_{h_{3}}=-Y^{3}_{ab}\overline{(X_{aR})^{c}}e_{bR}h^{+}_{3}+\mathrm{H.c.}=-Y^{3}_{ab}U^{\nu*}_{(a+6)i}\overline{(n_{i})}P_{R}e_{b}h^{+}_{3}+\mathrm{H.c.}.$
(59)

The new contributions to the cLFV decays and $\Delta
a^{331\mathrm{ISS}}_{\mu}$ is

$\displaystyle c^{h_{3}}_{(ab)R}$
$\displaystyle=\frac{em_{e_{a}}}{16\pi^{2}m_{e_{b}}m^{2}_{h_{3}}}\sum_{i=1}^{9}\sum_{c=1}^{3}Y^{3}_{ca}Y^{3*}_{cb}U^{\nu*}_{(a+6)i}U^{\nu}_{(b+6)i}\tilde{F}_{LHH}\left(\frac{m^{2}_{n_{i}}}{m^{2}_{h_{3}}}\right),$
$\displaystyle c^{h_{3}}_{(ba)R}$
$\displaystyle=\frac{e}{16\pi^{2}m^{2}_{h_{3}}}\sum_{i=1}^{9}\sum_{c=1}^{3}Y^{3*}_{ca}Y^{3}_{cb}U^{\nu}_{(a+6)i}U^{\nu*}_{(b+6)i}\tilde{F}_{LHH}\left(\frac{m^{2}_{n_{i}}}{m^{2}_{h_{3}}}\right).$
(60)

Although the contributions of these singly charged Higgs bosons to $\Delta
a_{\mu}$ are normally small and negative, the contributions to the cLFV
amplitudes may be significantly large. Consequently, they can affect
destructively the total cLFV decay amplitudes. These properties will keep
$\Delta a^{331\mathrm{ISS}}_{\mu}$ reaching the experimental constraint given
in Eq. (1), while keeping all other cLFV branching ratios well below the
experimental constraints. In this work, we consider the simplest case that
$h^{\pm}_{3}$ does not mix with the other singly charged Higgs bosons in the
331RN, and the mass is another free parameter. All of these properties can be
derived easily from the total Higgs potential, hence it will be ignored in
this work.

##  V  Numerical discussion

###  V.1 Without contributions from additional singly charged Higgs bosons
$h^{\pm}_{3}$

The numerical experimental parameters are taken from Ref. Zyla:2020zbs :

$\displaystyle G_{F}$ $\displaystyle=1.663787\times
10^{-5}\;\mathrm{GeV}^{-2},\;g=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_{\tau}=1.776\;\mathrm{GeV},\;m_{W}=80.385\;\mathrm{GeV},$
$\displaystyle\mathrm{Br}(\mu\rightarrow e\overline{\nu_{e}}\nu_{\mu})$
$\displaystyle\simeq 1.,\;\mathrm{Br}(\tau\rightarrow
e\overline{\nu_{e}}\nu_{\tau})\simeq
0.1782,\;\mathrm{Br}(\tau\rightarrow\mu\overline{\nu_{\mu}}\nu_{\tau})\simeq
0.1739.$ (61)

Before discussing on the allowed regions that satisfy all experimental
constraints of cLFV decays $e_{b}\rightarrow e_{a}\gamma$ as well as
$(g-2)_{\mu}$ data, we give some important crude estimation on the allowed
regions of parameter space constrained by both large $\Delta
a_{\mu}^{331\mathrm{ISS}}\geq\mathcal{O}(10^{-9})$ and small
Br$(e_{b}\rightarrow e_{a}\gamma)$. The way to derive the total mass matrix to
calculate numerically the masses and total neutrino mixing matrix $U^{\nu}$
were presented in the previous section. We have checked that the input changes
of $\Delta m^{2}_{ij}$ and $s^{2}_{ij}$ in the allowed ranges given in Eq.
(II.1) do not change significantly the final results, so we will fix these
quantities at their best-fit points. An exception that the Dirac phase
$\delta=180$ [Deg.] is considered so that the imagine parts of $c_{(ab)R}$ are
zeros, leading to a simple case of destruction among the one-loop
contributions from charged Higgs bosons.

In the numerical scan, the points in the allowed regions also satisfy
simultaneously the following conditions:

  1. 1.

The condition
$\mathrm{Re}[c^{H_{1}}_{(ab)R}]/\mathrm{Re}[c^{H_{2}}_{(ab)R}]<0$ with $a\neq
b$, will give a possibility that
$\mathrm{Re}[c^{H_{1}}_{(21)R}]+\mathrm{Re}[c^{H_{2}}_{(21)R}]\sim 0$, which
will result in valid regions of the parameter space in which two charged Higgs
bosons contributions can cancel each others. Therefore, these regions will
contain points which give the very suppressed total contributions to guarantee
the small Br$(e_{b}\rightarrow e_{a}\gamma)$. We will use this condition in
our numerical investigation.

  2. 2.

A crude numerical scan shows that the condition
$\mathrm{Re}[c^{H_{1}}_{(22)R}]/\mathrm{Re}[c^{H_{2}}_{(22)R}]>0$ so that the
two charged Higgs bosons contributions to $\Delta a^{331\mathrm{ISS}}_{\mu}$
always have the same sign, i.e., they give constructive contributions.
Therefore the values of $\Delta a^{331\mathrm{ISS}}_{\mu}$ are remained in the
original orders of $\mathcal{O}(10^{-9})$. Another case giving large $\Delta
a_{\mu}$ is that
$|\mathrm{Re}[c^{H_{i}}_{(22)R}]|\ll|\mathrm{Re}[c^{H_{j}}_{(22)R}]|$ with
$i\neq j$ when they have opposite signs.

First, we consider the simplest cases of all zero values of off-diagonal
elements $k_{ij}=0$ with $i\neq j$. The numerical investigation shows that we
cannot obtain any allowed points satisfying simultaneously both experimental
data of cLFV constraints and $\Delta a_{\mu}$. The reason is that there always
exists a strict relation that
$\mathrm{Re}[c^{H_{1}}_{(22)R}]/\mathrm{Re}[c^{H_{2}}_{(22)R}]$ and
$\mathrm{Re}[c^{H_{1}}_{(21)R}]/\mathrm{Re}[c^{H_{2}}_{(21)R}]$ are always
negative leading to small Br$(\mu\rightarrow e\gamma)$. As a consequence,
charged Higgs contributions to $\Delta a^{331\mathrm{ISS}}_{\mu}$ are always
destructive. Hence, the derived values are smaller than the experimental data.
A requirement of Br$(\mu\rightarrow e\gamma)\leq\mathcal{O}(10^{-8})$ gives
largest values of $\Delta a^{331\mathrm{ISS}}_{\mu}<10^{-9}$.

From a crude numerical scan, we can find the allowed regions of the parameter
space satisfying both conditions that Br$(\mu\rightarrow e\gamma)<4.2\times
10^{-13}$ and large $\Delta
a_{\mu}^{331\mathrm{ISS}}\geq\mathcal{O}(10^{-9})$. These allowed regions will
be used to collect the allowed points satisfying the remaining cLFV
constraints. The following ranges of the parameter space will be chosen as the
necessary conditions of free parameters when scanning to collect allowed
points:

$\displaystyle t_{\beta}\in[0.3,60],\;0.6\;\mathrm{[TeV]}\leq
m_{H_{1}},\;m_{H_{2}}\leq 3\;\mathrm{[TeV]},$ $\displaystyle\;|k_{ij}|\times
zc_{\beta}<\sqrt{4\pi}w=5.3\;\mathrm{[TeV]},\;10\;\mathrm{[GeV]}\leq z\leq
1223\;\mathrm{[TeV]}.$ (62)

Numerical values of $k_{ij}$ will be chosen so that they give active neutrino
masses and $U_{\mathrm{PMNS}}$ consistent with neutrino oscillation data. The
value of $5.3$ TeV is fixed from the lower bound $w$ obtained from the
experimental data of the heavy $Z^{\prime}$ boson mass $m_{Z^{\prime}}$. But
it can be relaxed with larger $w$ without any changes of final conclusions in
this work.

Without contributions of the additional singly charged Higgs boson, our
numerical investigation shows that the largest values of $\Delta a_{\mu}$
satisfying all cLFV constraints is $\Delta a_{\mu}\leq 108.5\times 10^{-11}$,
see an illustration shown in Fig. 1.

|
---|---
Figure 1:  The left panel shows $\Delta a_{\mu}$ vs. Br$(e_{b}\rightarrow
e_{a}\gamma)\sim\left(|c_{(ab)R}|^{2}+|c_{(ba)R}|^{2}\right)$ in the free
parameter ranges given in Table 2. The right panel shows
Br$(abR)\sim|c_{(ab)R}|^{2}$ with $a<b$.

The corresponding ranges of the free parameters are shown in Table 2,

Notation | $k_{11}$ | $k_{22}$ | $k_{33}$ | $k_{12}$ | $k_{13}$ | $k_{23}$ | $k_{21}$ | $k_{31}$ | $k_{32}$ | $t_{\beta}$ |  $z$ [GeV] |  $m_{H_{1}}$[GeV] |  $m_{H_{2}}$[GeV]
---|---|---|---|---|---|---|---|---|---|---|---|---|---
Min | -3.99 | -50.2 | 509. | -29.9 | 15.4 | -80.4 | 121. | 21.2 | 29.4 | 29.0 | 885. | 705 | 769
Max | 2.47 | -35.1 | 528. | -20.6 | 24.9 | -66.4 | 135. | 36.8 | 45.9 | 40.0 | 1150 | 893 | 962
Example | -3.26 | -49.7 | 509. | -28.6 | 23.0 | -77.8 | 124. | 25.3 | 37.1 | 36.9 | 969. | 754 | 825
Table 2:  Numerical values of free parameters for large $\Delta
a^{331\mathrm{ISS}}_{\mu}\geq 10^{-9}$ satisfying all experimental constraints
of the cLFV decays $e_{b}\rightarrow e_{a}\gamma$.

where the right panel shows the only contributions from $c_{(ab)R}$ ($a<b$) to
the decay rates, namely

$\mathrm{Br}(abR)=\frac{48\pi^{2}}{G_{F}^{2}}\left|c_{(ab)R}\right|^{2}\mathrm{Br}(e_{b}\rightarrow
e_{a}\overline{\nu_{a}}\nu_{b}).$

Here the two first lines show the respective minimum and maximum values of the
free parameters. The third line shows a particular example of the set of the
parameters giving large $\Delta a_{\mu}\simeq 108.1\times 10^{-11}$. The other
quantities are shown in Table 3, which will be discussed more later. In the
left panel of Fig. 1, only Br$(\tau\rightarrow\mu\gamma)$ always enhances with
increasing $\Delta a_{\mu}$. The upper constraint
Br$(\tau\rightarrow\mu\gamma)<4.4\times 10^{-8}$ gives the largest value of
$\Delta a_{\mu}\simeq 108.5\times 10^{-11}$. From the right panel of Fig. 1,
we see that $|c_{(23)R}|<|c_{(32)R}|$ in the region predicting large $\Delta
a_{\mu}$, because the contribution from $|c_{(23)R}|$ to
Br$(\tau\rightarrow\mu\gamma)$ denoted as Br$(23R)$ is small, namely
Br$(23R)\leq 0.2\times 10^{-8}$ with $\Delta a_{\mu}\geq 108\times 10^{-11}$.
This is in contrast to other normal cases, as we will discuss based on the
Table 3.

Table 3 illustrates particular values of $c_{(ab)R}$ and large $\Delta
a_{\mu}^{331\mathrm{ISS}}\simeq 108.1\times 10^{-11}$, corresponding to a set
of free parameters given in the third line of Table 2.

Notations | $c^{W}_{(ab)R}-c^{W,\mathrm{SM}}_{(ab)R}$ | $c^{Y}_{(ab)R}$ | $c^{H_{1}}_{(ab)R}$ | $c^{H_{2}}_{(ab)R}$ | $c_{(ab)R}$ | Process
---|---|---|---|---|---|---
$\Delta a_{\mu}:\;c_{(22)R}\times 10^{10}$ | 5.22 | -0.499 | -82.07 | 3.11 | -74.24 | $\Delta a_{\mu}=10.81\times 10^{-10}$
$\mu\rightarrow e\gamma:\;c_{(12)R}\times 10^{13}$ | 422.13 | 29.645 | -29086. | 28636. | 1.6960 | $\;\mathrm{Br}(12R)=1.002\times 10^{-13}$
$\mu\rightarrow e\gamma:\;c_{(21)R}\times 10^{13}$ | 2.010 | 0.1412 | -138.5 | 138.9 | 2.568 | $\;\mathrm{Br}(21R)=2.296\times 10^{-13}$
$\tau\rightarrow e\gamma:\;c_{(13)R}\times 10^{10}$ | -0.031 | 0.01941 | 13.60 | -15.63 | -2.039 | $\;\mathrm{Br}(13R)=257.9\times 10^{-10}$
$\tau\rightarrow e\gamma:\;c_{(31)R}\times 10^{10}$ | $\simeq 0$ | $\simeq 0$ | 0.004 | 0.031 | 0.035 | $\;\mathrm{Br}(31R)=0.076\times 10^{-10}$
$\tau\rightarrow\mu\gamma:\;c_{(23)R}\times 10^{10}$ | -0.02505 | -0.03235 | -0.3305 | 0.5170 | 0.1291 | $\;\mathrm{Br}(23R)=1.009\times 10^{-10}$
$\tau\rightarrow\mu\gamma:\;c_{(32)R}\times 10^{10}$ | -0.001481 | -0.001913 | -0.01954 | -2.656 | -2.679 | $\;\mathrm{Br}(32R)=434.7\times 10^{-10}$
Table 3:  Particular contributions $c^{X}_{(ab)R}\;[\mathrm{GeV}^{-2}]$ to the
$\Delta a_{\mu}$ and Br$(e_{b}\rightarrow e_{a}\gamma)$ with the free
parameters given in the third line of Table 2. The last column shows values of
$\Delta a_{\mu}$ and Br$(e_{b}\rightarrow e_{a}\gamma)$.

The numerical results given in Table 3 show that the experimental constraint
from Br$(\tau\rightarrow\mu\gamma)<4.4\times 10^{-8}$ does not allow large
$\Delta a_{\mu}^{331\mathrm{ISS}}>108.5\times 10^{-11}$. More particular,
$c_{(32)R}$ gives the dominant contribution to
Br$(\tau\rightarrow\mu\gamma)<4.4\times 10^{-8}$, with
$|c^{H_{2}}_{(32)R}|\gg|c^{H_{1}}_{(32)R}|$. In contrast, the remaining cLFV
decays have some common properties that $|c_{(ab)R}|>|c_{(ba)R}|$ with $a<b$,
$|c^{H_{k}}_{(ab)R}|\gg|c_{(ab)R}|$, and the huge destructive correlation
between two charged Higgs boson contributions. They are very important to
guarantee small Br$(\tau\rightarrow e\gamma)$ and Br$(\mu\rightarrow
e\gamma)$. On the other hand, they allow large and/or constructive
$c^{H_{k}}_{(22)R}$, which are the dominant contributions resulting in large
$\Delta a^{331\mathrm{ISS}}_{\mu}\geq 10^{-9}$.

The above properties are also true for the allowed region of the parameter
space given in Table 2. They are summarized in Table 4 through the quantities
defined in Eqs. (57) and (58).

| $R^{W}_{22}$ | $R^{H_{1}}_{22}$ | $R^{H_{2}}_{22}$ | $R^{W}_{12}$ | $R^{Y}_{12}$ | $R^{H_{1}}_{12}$ | $R^{H_{2}}_{12}$ | $R^{-}_{12}$ | $R^{W}_{21}$ | $R^{Y}_{21}$ | $R^{H_{1}}_{21}$ | $R^{H_{2}}_{21}$ | $R^{-}_{21}$ | $R^{W}_{13}$ | $R^{Y}_{13}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
Min | 0.06 | 0.96 | 0 | 54 | 7 | $\sim 10^{3}$ | $\sim 10^{3}$ | 60 | 0.3 | 0.04 | 34 | 34 | 0.01 | 0.01 | 0
Max | 0.08 | 1.11 | 0.1 | $\sim 10^{5}$ | $\sim 10^{4}$ | $\sim 10^{7}$ | $\sim 10^{7}$ | $\sim 10^{6}$ | 299 | 37 | $\sim 10^{4}$ | $\sim 10^{4}$ | 336 | 4. | 3
| $R^{H_{1}}_{13}$ | $R^{H_{2}}_{13}$ | $R^{-}_{13}$ | $R^{H_{1}}_{31}$ | $R^{H_{2}}_{31}$ | $R^{-}_{31}$ | $R^{W}_{23}$ | $R^{Y}_{23}$ | $R^{H_{1}}_{23}$ | $R^{H_{2}}_{23}$ | $R^{-}_{23}$ | $R^{H_{1}}_{32}$ | $R^{H_{2}}_{32}$ | $R^{-}_{32}$ |
Min | 6 | 7 | 0.05 | 0.1 | 0.8 | $\simeq 1$ | 0.01 | 0.03 | 0.5 | 0 | 0.02 | 0 | 0.96 | 0.998 |
Max | $\sim 10^{3}$ | $\sim 10^{3}$ | 2.6 | 0.2 | 0.89 | $\simeq 1$ | 16.4 | 15.4 | 222 | 255 | 33 | 0.04 | 0.996 | 1.002 |
Table 4:  Correlations between different contributions to $c_{(ab)R}$ with
ranges of free parameters given in Table 2, where we denote $0\simeq
R^{Y}_{(22)R},R^{W}_{(31)R},R^{Y}_{(31)R},R^{W}_{(32)R},R^{Y}_{(32)R}\leq\mathcal{O}(10^{-3})$
.

We can see that
$R^{-}_{32}=|\mathrm{Re}[c^{H_{2}}_{(32)R}+c^{H_{1}}_{(32)R}]/\mathrm{Re}[c_{(32)R}]|\rightarrow
1$ implies that sum of the two contributions of the charged Higgs bosons to
$c_{(32)R}$ is dominant. In addition
$R^{H_{2}}_{32}\equiv|\mathrm{Re}[c^{H_{2}}_{(32)R}]/\mathrm{Re}[c_{(32)R}]|\simeq
1$ indicates that the contributions of the charged Higgs boson $H_{2}$ is
dominant, hence the destructive correlation is small. This is not enough to
keep the cLFV constraint Br$(\tau\rightarrow\mu\gamma)<4.4\times 10^{-8}$ for
larger $\Delta a_{\mu}>108.5\times 10^{-11}$. All contributions of the two
decays $\mu\rightarrow e\gamma$ and $\tau\rightarrow e\gamma$ do not have
properties mentioned here. In the next discussion, we will show that new
destructive contributions from additional singly charged Higgs bosons will
relax the sum of the contributions from two charged Higgs bosons $H_{1,2}$ to
a larger values, while allow both $\Delta a^{331\mathrm{ISS}}_{\mu}$ and
Br$(\tau\to\mu\gamma)$ satisfying the experimental constraints given in Eqs.
(1) and (I).

From the above discussion, we can see that max$[\Delta a_{\mu}]\simeq
108.5\times 10^{-11}$ predicted by the 331ISS model comes from the
experimental constraints Br$(\tau\rightarrow\mu\gamma)<4.4\times 10^{-8}$,
which gets main contributions from $c_{(32)R}$. On the other hand,
Br$(\tau\rightarrow e\gamma)$ can reach to zero for large $\Delta a_{\mu}\geq
10^{-9}$, which is different from the normal behave of these two branching
ratios Br$(\tau\rightarrow\mu\gamma)\sim$Br$(\tau\rightarrow
e\gamma)\sim|c_{(23)R}|^{2},|c_{(13)R}|^{2}\gg|c_{(32)R}|^{2},|c_{(31)R}|^{2}$.
After some numerical checks, we see that the this difference is originated
mainly from the following property: each quantity
Br$(\tau\rightarrow\mu\gamma)$ or Br$(\tau\rightarrow e\gamma)$ contains only
one type of terms with a factor $\frac{m_{\mu}}{m_{\tau}}$ or
$\frac{m_{e}}{m_{\tau}}$ appearing in $c_{(32)R}$ or $c_{(31)R}$,
respectively. These terms are normally suppressed because of many other large
terms contained in $c^{H_{k}}_{(23)R},c^{H_{k}}_{(13)R}\gg
c^{H_{k}}_{(32)R},c^{H_{k}}_{(31)R}$. But when huge destructive correlations
between two charged Higgs contributions and gauge contributions happen, there
appears a situation that $|c_{(13)R}|,|c_{(23)R}|\rightarrow 0$, and also for
other normal large terms in $|c_{(31)R}|,|c_{(32)R}|$. Now, the terms with
factors $\frac{m_{\mu}}{m_{\tau}}$ and $\frac{m_{e}}{m_{\tau}}$ become
significant, leading to the consequence that Br$(\tau\rightarrow
e\gamma)\sim\frac{m^{2}_{e}}{m^{2}_{\tau}}$ can be close to 0, while
Br$(\tau\rightarrow\mu\gamma)\sim\frac{m^{2}_{\mu}}{m^{2}_{\tau}}\gg$Br$(\tau\rightarrow
e\gamma)$. It is reasonable to think that the terms with factor
$\frac{m_{\mu}}{m_{\tau}}$ and $\Delta a^{331\mathrm{ISS}}_{\mu}$ get similar
contributions relating to $\mu$, hence both of them must be large if $\Delta
a^{331\mathrm{ISS}}_{\mu}$ is required to be large in order to reach the
experimental constraints. Our explanation is confirmed by a numerical check,
where we change $m_{\mu}\to m_{e}$ in only the formula of $c_{(32)R}$. We saw
that Br$(\tau\to\mu\gamma,e\gamma)$ can reach small values
Br$(\tau\to\mu\gamma,e\gamma)<10^{-9}$ with $\Delta
a^{331\mathrm{ISS}}_{\mu}>125\times 10^{-11}$. Other numerical checks also
show that the lower bound of Br$(\tau\rightarrow\mu\gamma)$ depends strictly
on the lepton mixing matrix $U_{\mathrm{PMNS}}$, which is the only cLFV source
in the 331ISS model. First, the case of large $\tau-e$ mixing inputs
$s^{2}_{13}=s^{2}_{23}=0.547$ can give large $\Delta
a^{331\mathrm{ISS}}_{\mu}>115\times 10^{-11}$ and both small
Br$(\tau\rightarrow\mu\gamma,e\gamma)\to 0$. Second, the small input
$s^{2}_{23}=0.0216$ and the large input $s^{2}_{13}=0.547$ will result in that
max$[\Delta a^{331}_{\mu}]\simeq 90\times 10^{-11}$. In both cases,
max$[\Delta a^{331\mathrm{ISS}}_{\mu}]$ is still constrained by
Br$(\tau\to\mu\gamma)<4.4\times 10^{-8}$. In conclusion, the regions of the
parameter space giving max$[a^{331\mathrm{ISS}}_{\mu}]$ allows all small
$c_{(ab)R}$ except the terms with factor $\frac{m_{\mu}}{m_{\tau}}$ in
$c_{(32)R}$.

###  V.2  New contributions from additional singly charged Higgs bosons
$h^{\pm}_{3}$

Adding contributions of the new singly charged Higgs boson, the allowed values
of $\Delta a_{\mu}\equiv\Delta a^{331\mathrm{ISS}}_{\mu}\geq 192\times
10^{-11}$ corresponding to the lower bound of the $1\sigma$ confidence level
are explained successfully, see an illustration shown in Fig. 2,

|
---|---
Figure 2:  Correlations between $\Delta a_{\mu}\equiv\Delta
a^{331\mathrm{ISS}}_{\mu}$ with $\Delta a_{\mu}(h_{3})$ and
Br$(\tau\rightarrow\mu\gamma)[h_{3}]$.

where $\Delta a_{\mu}(h_{3})$ and Br$(\tau\rightarrow\mu\gamma)[h_{3}]$ show
the respective one-loop contributions from only $h^{\pm}_{3}$ to $\Delta
a_{\mu}$ and Br$(\tau\rightarrow\mu\gamma)$, which are defined as follows:

$\displaystyle\Delta a_{\mu}[h_{3}]$
$\displaystyle=-\frac{4m^{2}_{\mu}}{e}\mathrm{Re}[c^{h_{3}}_{(22)R}],$
$\displaystyle\mathrm{Br}(e_{b}\rightarrow e_{a}\gamma)[h_{3}]$
$\displaystyle=\frac{48\pi^{2}}{G_{F}^{2}}\left(\left|c^{h_{3}}_{(ab)R}\right|^{2}+\left|c^{h_{3}}_{(ba)R}\right|^{2}\right)\mathrm{Br}(e_{b}\rightarrow
e_{a}\overline{\nu_{a}}\nu_{b}).$ (63)

The corresponding benchmark is calculated numerically with 30 digits of
precision number. The numerical values of the free parameters are

$\displaystyle k_{11}$
$\displaystyle\simeq-19.19,\;k_{22}\simeq-94.53,\;k_{33}\simeq
428.75,\;k_{12}\simeq-89.46,$ $\displaystyle k_{13}$ $\displaystyle\simeq
29.47,\;k_{23}\simeq-211.84,\;k_{21}\simeq
60.09,\;k_{31}\simeq-262.44,\;k_{32}\simeq 30.53,$ $\displaystyle t_{\beta}$
$\displaystyle=49.86,\;z=1169\;\mathrm{GeV},\;m_{H_{1}}=657.1\;\mathrm{GeV},\;m_{H_{2}}=734\;\mathrm{GeV}.$
(64)

In this case, the heavy neutrino masses are $m_{n_{4}}=m_{n_{5}}=137.2$ GeV,
$m_{n_{6}}=m_{n_{7}}=4709.4$ GeV, $m_{n_{8}}=m_{n_{9}}=11958$ GeV. For
simplicity we assume that
$Y^{3}_{11}=Y^{3}_{12}=Y^{3}_{21}=Y^{3}_{13}=Y^{3}_{31}=0$, therefore the
contribution from $h_{3}$ does not change the two cLFV decays
Br$(\mu\rightarrow e\gamma)\simeq 3.93\times 10^{-13}$ and Br$(\tau\rightarrow
e\gamma)\simeq 1.11\times 10^{-8}$. They always satisfy the experimental data.
The non-zero Yukawa couplings are scanned in the ranges
$Y^{3}_{ab}\in[-3.5,3.5]$ that satisfy the perturbative limit. This results in
the following allowed range of the charged Higgs boson mass
$500\;\mathrm{GeV}\leq m_{h_{3}}\leq\;1158\mathrm{GeV}$. Numerical values of
$c_{(ab)R}$ is shown in Table 5.

Notations | $c^{W}_{(ab)R}-c^{W,\mathrm{SM}}_{(ab)R}$ | $c^{Y}_{(ab)R}$ | $c^{H_{1}}_{(ab)R}$ | $c^{H_{2}}_{(ab)R}$ | $c^{h_{3}}_{(ab)R}$ | $c_{(ab)R}$ | Process
---|---|---|---|---|---|---|---
$\Delta a_{\mu}:\;c_{(22)R}\times 10^{10}$ | 5.3 | -0.386 | -211. | 61.1 | 3.7 | -141.1 | $\Delta a_{\mu}=20.5\times 10^{-10}$
$\mu\rightarrow e\gamma:\;c_{(12)R}\times 10^{13}$ | 449.16 | 61.536 | -75957. | 75443. | 0 | -2.5234 | $\;\mathrm{Br}(12R)=2.2174\times 10^{-13}$
$\mu\rightarrow e\gamma:\;c_{(21)R}\times 10^{13}$ | 2.1388 | 0.29303 | -361.70 | 357.43 | 0 | -1.8329 | $\;\mathrm{Br}(21R)=1.1699\times 10^{-13}$
$\tau\rightarrow e\gamma:\;c_{(13)R}\times 10^{10}$ | -0.00510 | 0.0540 | 4.25 | -2.96 | 0 | 1.34 | $\;\mathrm{Br}(13R)=111.\times 10^{-10}$
$\tau\rightarrow e\gamma:\;c_{(31)R}\times 10^{10}$ | $\sim 0$ | $\sim 0$ | 0.00120 | 0.0664 | 0 | 0.0676 | $\;\mathrm{Br}(31R)=0.284\times 10^{-10}$
$\tau\rightarrow\mu\gamma:\;c_{(23)R}\times 10^{10}$ | -0.00721 | -0.0445 | 1.20 | -2.51 | 0.164 | -1.20 | $\;\mathrm{Br}(23R)=86.7\times 10^{-10}$
$\tau\rightarrow\mu\gamma:\;c_{(32)R}\times 10^{10}$ | -0.000426 | -0.00263 | 0.0708 | -5.18 | 2.77 | -2.33 | $\;\mathrm{Br}(32R)=330.\times 10^{-10}$
Table 5:  Particular contributions $c^{X}_{(ab)R}[\mathrm{GeV}^{-2}]$ to
$\Delta a_{\mu}$ and Br$(e_{b}\rightarrow e_{a}\gamma)$ with the free
parameters shown in Eq. (V.2). The last column shows values of $\Delta
a_{\mu}$ and Br$(e_{b}\rightarrow e_{a}\gamma)$.

The numerical results shown in Fig. 2 have some interesting properties. In the
left panel, the contributions from $h^{\pm}_{3}$ to $\Delta a_{\mu}$ are
always negative, but much smaller than the total one: $0<-\Delta
a_{\mu}(h^{\pm}_{3})\leq 1.5\times 10^{-10}\ll 200\times 10^{-11}\sim\Delta
a_{\mu}$. On the other hand, the one-loop contributions $c^{h_{3}}_{(32)R}$
and $c^{H_{2}}_{(32)R}$ have the same order, but opposite signs. Therefore,
the total $|c_{(32)R}|$ is small enough to guarantee that
Br$(\tau\rightarrow\mu\gamma)<4.4\times 10^{-8}$. This is reason why in the
right panel, we see that $|c_{(32)R}|<|c^{h_{3}}_{(32)R}|$, i.e.
Br$(\tau\rightarrow\mu\gamma)<$ Br$(\tau\rightarrow\mu\gamma)[h_{3}]$ may
happen. More specifically, this property can be seen from a particular
numerical illustration presented in Table 5. We can see a property that
$|c_{(22)R}|\gg|c^{h_{3}}_{(22)R}|\sim|c^{h_{3}}_{(32)R}|\sim|c^{H_{2}}_{(32)R}|\sim|c_{(32)R}|$,
which explains why the contributions from $h_{3}$ affect strongly
Br$(\tau\rightarrow\mu\gamma)$ but weakly $\Delta a_{\mu}$.

The allowed regions of parameters allowing $\Delta a_{\mu}^{331\mathrm{ISS}}$
around the value $200\times 10^{-11}$ can be found easily in the ranges given
in Eq. (V.1). The allowed regions with larger $\Delta
a_{\mu}^{331\mathrm{ISS}}$ are shown in Fig. 3, where charged Higgs masses
have to be smaller than 600 GeV.

|
---|---
Figure 3:  Correlations between $\Delta a_{\mu}\equiv\Delta
a^{331\mathrm{ISS}}_{\mu}\geq 240\times 10^{-11}$ with $\Delta a_{\mu}(h_{3})$
and Br$(\tau\rightarrow\mu\gamma)[h_{3}]$.

It is noted that large $\Delta a_{\mu}^{331\mathrm{ISS}}>300\times 10^{-11}$
require light charged Higgs boson masses $m_{H_{1}}\to 500$ GeV, $z\rightarrow
1223$ GeV, and large $t_{\beta}\to 60$. The region of parameter space
corresponding to the Fig 3 is:

$\displaystyle k_{11}$
$\displaystyle\in\left[-21.77,\;-17.84\right],\;k_{22}\in\left[-101.9,\;-93.76\right],\;k_{33}\in\left[420.1,\;429.4\right],$
$\displaystyle k_{12}$
$\displaystyle\in\left[-96.22,\;-88.92\right],\;k_{13}\in\left[26.95,\;31.12\right],\;k_{23}\in\left[-220.2,\;-210.4\right],$
$\displaystyle k_{21}$
$\displaystyle\in\left[59.19,\;66.55\right],\;k_{31}\in\left[-268.6,\;-262.9\right],\;k_{32}\in\left[25.35,\;33.64\right],$
$\displaystyle t_{\beta}$
$\displaystyle\in\left[41.68,\;59.97\right],\;z\in\left[1051,\;1223\right]\;\mathrm{GeV},\;m_{H_{1}}\in\left[500.6,\;631.3\right]\;\mathrm{GeV},$
$\displaystyle m_{H_{2}}$
$\displaystyle\in\left[571.3,\;703.8\right]\;\mathrm{GeV},\;m_{h_{3}}\in\left[500.5,\;778.6\right]\;\mathrm{GeV},\;|Y_{22}|\in\left[0.11,\;3.49\right],$
$\displaystyle|Y_{23}|$
$\displaystyle\in\left[0.51,\;3.5\right],\;|Y_{32}|\in\left[0.06,\;3.49\right],\;|Y_{33}|\in\left[0.009,\;3.5\right].$
(65)

The heavy neutrino masses are in the following ranges:
$m_{n_{4}}=m_{n_{5}}\in\left[109.2,\;172.3\right]$ GeV,
$m_{n_{6}}=m_{n_{7}}\in\left[3.66,\;5.87\right]$ TeV,
$m_{n_{8}}=m_{n_{9}}\in\left[8.99,\;14.92\right]$ TeV. The cLFV branching
ratios are in the following ranges: Br$(\mu\rightarrow e\gamma)\times
10^{13}\in[5.8\times 10^{-16},\;4.2\times 10^{-13}]$, Br$(\tau\rightarrow
e\gamma)\in[4\times 10^{-11},\;3.3\times 10^{-8}]$, and
Br$(\tau\rightarrow\mu\gamma)\in[1.6\times 10^{-12},\;4.4\times 10^{-8}]$.

##  VI Conclusion

In this work, we have pointed out that the one of the versions of the 3-3-1RN
model, namely the 331ISS model, can predict large values of $\Delta
a_{\mu}\simeq 108\times 10^{-11}$ under the recent constraint of all cLFV
decays $e_{b}\rightarrow e_{a}\gamma$. This large value corresponds to the
upper bound Br$(\tau\rightarrow\mu\gamma)\simeq 4.4\times 10^{-8}$, while the
two remaining decay branching ratios are still well below the recent
experimental constraints. This model predicts the existence of the two charged
Higgs bosons which can give large contributions of the order
$\mathcal{O}(10^{-9})-\mathcal{O}(10^{-8})$ to the $\Delta a_{\mu}$, so that
it can reach the maximal values around $10^{-9}$, which is still much smaller
than the allowed values given by the recent experimental data. On the other
hand, the two other charged Higgs bosons contributions to Br($e_{b}\rightarrow
e_{a}\gamma$) will be at the orders of
$\mathcal{O}(10^{-10})-\mathcal{O}(10^{-9})[\mathrm{GeV}^{-2}]$. But the huge
destructive correlations can happen between these contributions, leading to a
small values of Br$(e_{b}\rightarrow e_{a}\gamma)$. Although the model
contains many free parameters, maybe the antisymmetry of the Dirac mass matrix
$m_{D}$ does not allow large destruction enough to keep the
Br$(\tau\rightarrow\mu\gamma)$ below the experimental constraint, while allow
large $\Delta a_{\mu}^{331\mathrm{ISS}}\geq 192\times 10^{-11}$. The model
needs to include an additional singly charged Higgs boson so that all
experimental data of $\Delta a_{\mu}$ and the cLFV decays can be explained
simultaneously. As a consequence, all of the cLFV decays $e_{b}\rightarrow
e_{a}\gamma$ are predicted that their branching ratios can be large closely
the recent experimental bounds. Therefore, our model can also explain
simultaneously all cLFV decays $e_{b}\rightarrow e_{a}\gamma$ once they are
observed by upcoming experiments.

## Acknowledgments

We are grateful to Prof. Martin Hoferichter for introducing us the latest
result of the SM prediction of $\Delta a_{\mu}$. We thank Prof. Hidezumi
Terazawa, Dr. Wen Yin, and Dr. Pengxuan Zhu for their useful information. We
would like to express our sincere gratitude to the referee for correcting the
electron mass in the original draft, leading to the new numerical illustration
in the new version. This research is funded by Vietnam National Foundation for
Science and Technology Development (NAFOSTED) under grant number
103.01-2018.331.

## References

  * (1) K. Hagiwara, R. Liao, A. D. Martin, D. Nomura and T. Teubner, J. Phys. G 38, 085003 (2011) [arXiv:1105.3149 [hep-ph]].
  * (2) M. Davier, A. Hoecker, B. Malaescu, and Z. Zhang, Eur. Phys. J. C 77, 827 (2017), arXiv:1706.09436.
  * (3) R. H. Parker, C. Yu, W. Zhong, B. Estey, and H. Müller, Science 360 191 (2018), arXiv:1812.04130,
  * (4) P. A. Zyla et al. [Particle Data Group], Prog.Theor. Exp. Phys. 2020, 083C01 (2020)
  * (5) A. Keshavarzi, D. Nomura and T. Teubner, Phys. Rev. D 97, 114025 (2018) [arXiv:1802.02995 [hep-ph]].
  * (6) G. Colangelo, M. Hoferichter and P. Stoffer, J. High. Energy. Phys. 02 (2019) 006 [arXiv:1810.00007 [hep-ph]].
  * (7) M. Hoferichter, B. L. Hoid and B. Kubis, J. High. Energy. Phys. 08 (2019) 137 [arXiv:1907.01556 [hep-ph]].
  * (8) M. Davier, A. Hoecker, B. Malaescu and Z. Zhang, Eur. Phys. J. C 80, 241 (2020); 80, 410(E) (2020) [arXiv:1908.00921 [hep-ph]].
  * (9) A. Keshavarzi, D. Nomura and T. Teubner, Phys. Rev. D 101, 014029 (2020) [arXiv:1911.00367 [hep-ph]].
  * (10) A. Kurz, T. Liu, P. Marquard and M. Steinhauser, Phys. Lett. B 734, 144 (2014) [arXiv:1403.6400 [hep-ph]].
  * (11) K. Melnikov and A. Vainshtein, Phys. Rev. D 70, 113006 (2004) [arXiv:hep-ph/0312226 [hep-ph]].
  * (12) P. Masjuan and P. Sanchez-Puertas, Phys. Rev. D 95, 054026 (2017) [arXiv:1701.05829 [hep-ph]].
  * (13) G. Colangelo, M. Hoferichter, M. Procura and P. Stoffer, J. High. Energy. Phys. 04 (2017) 161 [arXiv:1702.07347 [hep-ph]].
  * (14) M. Hoferichter, B. L. Hoid, B. Kubis, S. Leupold and S. P. Schneider, J. High. Energy. Phys. 10 (2018) 141 [arXiv:1808.04823 [hep-ph]].
  * (15) A. Gérardin, H. B. Meyer and A. Nyffeler, Phys. Rev. D 100, 034520 (2019) [arXiv:1903.09471 [hep-lat]].
  * (16) J. Bijnens, N. Hermansson-Truedsson and A. Rodríguez-Sánchez, Phys. Lett. B 798, 134994 (2019) [arXiv:1908.03331 [hep-ph]].
  * (17) G. Colangelo, F. Hagelstein, M. Hoferichter, L. Laub and P. Stoffer, J. High. Energy. Phys. 03 (2020) 101 [arXiv:1910.13432 [hep-ph]].
  * (18) G. Colangelo, M. Hoferichter, A. Nyffeler, M. Passera and P. Stoffer, Phys. Lett. B 735, 90 (2014) [arXiv:1403.7512 [hep-ph]].
  * (19) T. Blum, N. Christ, M. Hayakawa, T. Izubuchi, L. Jin, C. Jung and C. Lehner, Phys. Rev. Lett. 124, 132002 (2020) [arXiv:1911.08123 [hep-lat]].
  * (20) T. Aoyama, M. Hayakawa, T. Kinoshita and M. Nio, Phys. Rev. Lett. 109, 111808 (2012) [arXiv:1205.5370 [hep-ph]].
  * (21) T. Aoyama, T. Kinoshita and M. Nio, Atoms 7, 28 (2019)
  * (22) A. Czarnecki, W. J. Marciano and A. Vainshtein, Phys. Rev. D 67, 073006 (2003); 73, 119901(E) (2006) [arXiv:hep-ph/0212229 [hep-ph]].
  * (23) C. Gnendiger, D. Stöckinger and H. Stöckinger-Kim, Phys. Rev. D 88, 053005 (2013) [arXiv:1306.5546 [hep-ph]].
  * (24) T. Aoyama, N. Asmussen, M. Benayoun, J. Bijnens, T. Blum, M. Bruno, I. Caprini, C. M. Carloni Calame, M. Cè, G. Colangelo, et al. Phys. Rep. 887, 1 (2020) [arXiv:2006.04822 [hep-ph]].
  * (25) B. Abi et al. (Muon g-2 Collaboration), Phys. Rev. Lett. 126, 141801 (2021) [arXiv:2104.03281 [hep-ex]].
  * (26) B. Aubert et al. (BABAR Collaboration), Phys. Rev. Lett. 104, 021802 (2010) arXiv:0908.2381 [hep-ex].
  * (27) A. M. Baldini et al. (MEG CollaboratioN), Eur. Phys. J. C 76, 434 (2016) [arXiv:1605.05081 [hep-ex]].
  * (28) N. A. Ky, H. N. Long and D. V. Soa, Phys. Lett. B 486, 140 (2000), arXiv:hep-ph/0007010 [hep-ph].
  * (29) C. Kelso, H. N. Long, R. Martinez and F. S. Queiroz, Phys. Rev. D 90, 113011 (2014) arXiv:1408.6203 [hep-ph].
  * (30) Dinh Thanh Binh, Do Thi Huong, Le Tho Hue, and Hoang Ngoc Long, Commun. Phys. 25, 29 (2015)
  * (31) A. S. De Jesus, S. Kovalenko, F. S. Queiroz, C. Siqueira and K. Sinha, Phys. Rev. D 102, 035004 (2020) arXiv:2004.01200 [hep-ph].
  * (32) Á. S. de Jesus, S. Kovalenko, C. A. de S. Pires, F. S. Queiroz and Y. S. Villamizar, Phys. Lett. B 809, 135689 (2020) arXiv:2003.06440 [hep-ph].
  * (33) M. Lindner, M. Platscher and F. S. Queiroz, Phys. Rep. 731, 1 (2018) [arXiv:1610.06587 [hep-ph]].
  * (34) A. E. Cárcamo Hernández, D. T. Huong and H. N. Long, Phys. Rev. D 102, 055002 (2020) [arXiv:1910.12877 [hep-ph]].
  * (35) A. E. Cárcamo Hernández, Y. Hidalgo Velásquez, S. Kovalenko, H. N. Long, N. A. Pérez-Julve and V. V. Vien, Eur. Phys. J. C 81, 191 (2021) [arXiv:2002.07347 [hep-ph]].
  * (36) Le Tho Hue, Nguyen Thanh Phong, and Tran Dinh Tham, Commun. Phys. 30, 221 (2020).
  * (37) B. Sánchez-Vega, E. Schmitz and J. Montero, Eur. Phys. J. C 78, 166 (2018) arXiv:1612.03827 [hep-ph].
  * (38) D. Chang and H. N. Long, Phys. Rev. D 73, 053006 (2006) arXiv:hep-ph/0603098 [hep-ph].
  * (39) M. E. Catano, R. Martinez and F. Ochoa, Phys. Rev. D 86, 073015 (2012) [arXiv:1206.1966 [hep-ph]].
  * (40) A. G. Dias, C. A. de S.Pires, P. S. Rodrigues da Silva and A. Sampieri, Phys. Rev. D 86, 035007 (2012) [arXiv:1206.2590 [hep-ph]].
  * (41) S. M. Boucenna, J. W. F. Valle and A. Vicente, Phys. Rev. D 92, 053001 (2015) arXiv:1502.07546 [hep-ph].
  * (42) C. A. de Sousa Pires, F. Ferreira De Freitas, J. Shu, L. Huang and P. Wagner Vasconcelos Olegário, Phys. Lett. B 797, 134827 (2019) arXiv:1812.10570 [hep-ph].
  * (43) A. M. Baldini, F. Cei, C. Cerri, S. Dussoni, L. Galli, M. Grassi, D. Nicolo, F. Raffaelli, F. Sergiampietri, G. Signorelli, et al. [arXiv:1301.7225 [physics.ins-det]].
  * (44) T. Aushev, W. Bartel, A. Bondar, J. Brodzicka, T. E. Browder, P. Chang, Y. Chao, K. F. Chen, J. Dalseno, A. Drutskoy, et al. [arXiv:1002.5012 [hep-ex]].
  * (45) J. Cao, J. Lian, L. Meng, Y. Yue and P. Zhu, Phys. Rev. D 101, 095009 (2020) [arXiv:1912.10225 [hep-ph]].
  * (46) J. Cao, Y. He, J. Lian, D. Zhang and P. Zhu, [arXiv:2102.11355 [hep-ph]].
  * (47) T. Nomura, H. Okada and P. Sanyal, [arXiv:2103.09494 [hep-ph]].
  * (48) T. Mondal and H. Okada, [arXiv:2103.13149 [hep-ph]].
  * (49) T.Phong Nguyen, T. Thuy Le, T. T. Hong, and L. T. Hue, Phys. Rev. D 97, 073003 (2018) 073003, arXiv:1802.00429 [hep-ph].
  * (50) R. Foot, H. N. Long, and T. A. Tran, Phys. Rev. D 50, R34 (1994) arXiv:hep-ph/9402243 [hep-ph].
  * (51) A. J. Buras, F. De Fazio, J. Girrbach and M. V. Carlucci, J. High. Energy. Phys. 02 (2013) 023 arXiv:1211.1237 [hep-ph].
  * (52) A. J. Buras, F. De Fazio and J. Girrbach-Noe, J. High Energy Phys. 08 (2014) 039 [arXiv:1405.3850 [hep-ph]].
  * (53) L. Hue and L. Ninh, Eur. Phys. J. C 79, 221 (2019) [arXiv:1812.07225 [hep-ph]].
  * (54) L. T. Hue, H. N. Long, T. T. Thuc and T.Phong Nguyen, Nucl. Phys. B907, 37 (2016) arXiv:1512.03266 [hep-ph].
  * (55) A. E. Cárcamo Hernández, S. Kovalenko, H. N. Long and I. Schmidt, J. High. Energy. Phys. 07, 144 (2018) [arXiv:1705.09169 [hep-ph]].
  * (56) M. Tanabashi et al. [Particle Data Group], Phys. Rev. D 98, 030001 (2018).
  * (57) K. Abe et al. (T2K Collaboration), Nature (London) 580, 339 (2020); 583, E16 (2020) arXiv:1910.03887 [hep-ex].
  * (58) W. Yin, [arXiv:2103.14234 [hep-ph]].
  * (59) M. J. Baker, P. Cox and R. R. Volkas, J. High. Energy. Phys. 05 (2021) 174 [arXiv:2103.13401 [hep-ph]].
  * (60) L. D. Ninh and H. N. Long, Phys. Rev. D 72, 075004 (2005) arXiv:hep-ph/0507069 [hep-ph].
  * (61) L. T. Hue, L. D. Ninh, T. T. Thuc and N. Dat, Eur. Phys. J. C 78, 128 (2018) [arXiv:1708.09723 [hep-ph]].
  * (62) A. Crivellin, M. Hoferichter and P. Schmidt-Wellenburg, Phys. Rev. D 98, 113002 (2018) [arXiv:1807.11484 [hep-ph]].
  * (63) G. W. Bennett et al. (Muon (g-2) Collboration), Phys. Rev. D 80, 052008 (2009) [arXiv:0811.1207 [hep-ex]].
  * (64) F. Jegerlehner and A. Nyffeler, Phys. Rep. 477, 1 (2009) [arXiv:0902.3360 [hep-ph]].
  * (65) T. Aoyama, M. Hayakawa, T. Kinoshita and M. Nio, Phys. Rev. Lett. 109, 111807 (2012) [arXiv:1205.5368 [hep-ph]].
  * (66) S. Laporta, Phys. Lett. B 772, 232 (2017) [arXiv:1704.06996 [hep-ph]].
  * (67) H. Terazawa, Nonlin. Phenom. Complex Syst. 21, 268 (2018)
  * (68) S. Volkov, Phys. Rev. D 100, 096004 (2019) [arXiv:1909.08015 [hep-ph]].