#  $(g-2)_{e,\mu}$ and decays $e_{b}\to e_{a}\gamma$ in a
$\mbox{SU}(4)_{L}\otimes\mbox{U}(1)_{X}$ model with inverse seesaw neutrinos

N. H. Thao<EMAIL_ADDRESS>Department of Physics, Hanoi Pedagogical
University 2, no 32 Nguyen Van Linh, Phuc Yen, Vinh Phuc, Vietnam  D. T. Binh
<EMAIL_ADDRESS>Institute of Theoretical and Applied Research,
Duy Tan University, Hanoi, Vietnam  Faculty of Natural Science, Duy Tan
University, Da Nang, Vietnam  T. T. Hong<EMAIL_ADDRESS>An Giang
University, VNU - HCM, Ung Van Khiem Street, Long Xuyen, An Giang, Vietnam  L.
T. Hue<EMAIL_ADDRESS>Subatomic Physics Research Group, Science and
Technology Advanced Institute, Van Lang University, Ho Chi Minh City, Vietnam
Faculty of Applied Technology, School of Technology, Van Lang University, Ho
Chi Minh City, Vietnam  D. P. Khoi 111corresponding author
<EMAIL_ADDRESS>Department of Physics, Vinh University, 182 Le Duan,
Vinh City, Nghe An, Vietnam

###### Abstract

In this paper, we will show that the 3-4-1 model with heavy right-handed
neutrinos can explain the recent experimental data of $(g-2)_{\mu,e}$ of
charged leptons and neutrino oscillations through the inverse seesaw
mechanism. In addition, the model can predict large lepton flavor violating
decay rates $\mu\to e\gamma$ and $\tau\to\mu\gamma,e\gamma$ up to the recent
experimental sensitivities.

##  I Introduction

The 3-4-1 model with right-handed neutrinos (341RHN) was discussed in Ref.
Foot:1994ym ; Pisano:1994tf  as a natural extension that new right-handed
neutrinos are assigned into the same left-handed lepton quadruplets. For a
complete study of the highest possible gauge symmetry in the electroweak
sector Voloshin:1987qy , various 3-4-1 extensions were introduced with
different electric charges of new exotic leptons Ponce:2003uu ; Sanchez:2004uf
; Ponce:2006vw ; Sanchez:2008qv ; Riazuddin:2008yx ; Nam:2009tr ; Long:2016lmj
; Palacio:2016mam ; Djouala:2019hhn . It was proved that these original 3-4-1
models can not explain the recent data of the $(g-2)$ anomaly of muon unless
they must be extended such as adding new inert scalars Cogollo:2020nrc . A
solution applied to the 3-3-1 models Hue:2021zyw , that adding new singly
charged Higgs bosons and inverse seesaw (ISS) neutrinos, is another viable
approach. In this work, we will investigate the possibility of whether this
approach can work in the 3-4-1 model framework, which can accommodate data of
charged lepton anomalies $a_{e_{a}}\equiv(g-2)_{e_{a}}/2$, neutrino
oscillation data, and the recent bounds of the lepton flavor violating decays
of charged leptons (cLFV) $e_{b}\to e_{a}\gamma$. The explanations of neutrino
oscillation data were mentioned previously in various 3-4-1 models, including
the ISS mechanism Palcu:2009uk ; Palcu:2015ica ; Palacio:2016mam , but they
did not relate to the $(g-2)$ data and cLFV decays. In the ISS models, new
gauge contributions to $(g-2)$ are suppressed Pinheiro:2021mps ; Palcu:2015ica
, especially the 3-3-1 and 3-4-1 models, because the new gauge bosons must be
heavy to guarantee the recent lower bounds from experimental searches
Lee:2014kna . Therefore, we will study the appearance of new singly charged
Higgs bosons and their mixing with the $SU(4)_{L}$ ones that results in large
chirally-enhanced one-loop Higgs contributions to $a_{e_{a}}$ enough to be
consistent with experiments Crivellin:2018qmi . On the other hand, the ISS
mechanism may lead to large one-loop contributions to cLFV rates. In this
work, the numerical investigations to determine the allowed regions of
parameter space satisfying all experimental constraints of $(g-2)$ anomaly and
cLFV decays $e_{b}\to e_{a}\gamma$ will be discussed precisely in the 3-4-1
framework. Our work will confirm the reality of the 3-4-1 models, besides many
available models beyond the SM Han:2018znu ; Endo:2019bcj ; Li:2020dbg ;
DelleRose:2020oaa ; Bigaran:2020jil ; Botella:2020xzf ; Bharadwaj:2021tgp ;
Han:2021gfu ; Arbelaez:2020rbq ; Chun:2020uzw ; Chen:2020tfr ; Dutta:2020scq ;
Wang:2022yhm ; Hernandez:2021tii ; Hernandez:2021xet ; Li:2022zap ;
Botella:2022rte ; Kriewald:2022erk ; Barman:2021xeq ; Dermisek:2022hgh ;
Chowdhury:2022jde ; Chen:2023eof  under various experimental data including
the $(g-2)_{e,\mu}$ anomalies.

The latest experimental measurement for muon anomaly $a_{\mu}$ has been
reported from the combination of the two experimental results at Fermilab
Muong-2:2021ojo  and Brookhaven National Laboratory (BNL) E82 Muong-2:2006rrc
: $a^{\mathrm{exp}}_{\mu}=116592061(41)\times 10^{-11}$. It leads to a
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 ; Danilkin:2016hnh ; 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 ; Pauk:2014rta ; Jegerlehner:2017gek ;
Knecht:2018sci ; Eichmann:2019bqf ; Roig:2019reh . 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 ;
Laporta:2017okg ; Aoyama:2017uqe ; Terazawa:2018pdc ; Volkov:2019phy ;
Gerardin:2020gpp . In this work, we accept the following value:

$\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)

The ISS mechanism may result in large values of not only $(g-2)_{e,\mu}$ but
also Br$(e_{b}\to e_{a}\gamma)$, which are constrained by recent experiments
BaBar:2009hkt ; MEG:2016leq ; Belle:2021ysv

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

The future sensitivities for these decays are Br$(\mu\to e\gamma)<6\times
10^{-14}$ MEGII:2018kmf  and Br$(\tau\to\mu\gamma)<2.7\times 10^{-8}$ Belle-
II:2018jsg . Hence the correlations between $a_{e_{a}}$ and cLFV decays
$e_{b}\to e_{a}\gamma$ may give new predictions on both of these kinds of
processes, namely whether large $a_{e_{a}}$ will exclude the more strict
experimental constraint of cLFV decays in the near future.

Our paper is arranged as follows. Section II presents all ingredients of a
3-4-1 model to calculate the $(g-2)_{e_{a}}$ data and cLFV decays. Section III
introduces the 341ISS model to accommodate the recent $(g-2)_{e_{a}}$ data. In
this model, the Yukawa Lagrangian and Higgs potential must respect a new
global $U(1)_{\mathcal{L}}$ symmetry in order to guarantee the appearance of
the ISS mechanism, the mixing between singly charged Higgs bosons, and the
Yukawa couplings resulting in large chirally-enhanced one-loop contributions
to $(g-2)_{e_{a}}$ anomalies. Section IV will present detailed numerical
results to determine the allowed regions of the parameter space that explain
both experimental results of two $(g-2)_{e,\mu}$ anomalies and cLFV decays.
Section V summarizes important results.

##  II The model with Dirac active neutrinos

###  II.1 Yukawa couplings and masses for fermions

In this work, we will study the 3-4-1 model with heavy right-handed neutrinos
and new singly charged leptons assigned in the three left-handed quadruplets
Sanchez:2004uf ; Palacio:2016mam . The electric operator is defined as
follows:

$Q=T_{3}+\frac{1}{\sqrt{3}}T_{8}-\frac{2}{\sqrt{6}}T_{15}+X\mathbb{I}.$ (4)

The lepton sector consists of three left-handed quadruplets and respective
right-handed singlets, namely

$\displaystyle
L_{a}=(\nu^{\prime}_{a}\,,e^{\prime}_{a}\,,E^{\prime}_{a}\,,N^{\prime}_{a})_{L}^{T}\sim\left(1,4,-\frac{1}{2}\right)\,,\;a=1,2,3,$
$\displaystyle e^{\prime}_{aR},\;E^{\prime}_{aR}\sim(1,1,-1),\hskip
14.22636pt\nu^{\prime}_{aR},\;N^{\prime}_{aR}\sim(1,1,0).$ (5)

The Higgs multiplets and non-zero vacuum expectation values (VEV) of neutral
components needed for generating all fermion masses are:

$\displaystyle\chi$
$\displaystyle=\left(\chi_{1}^{0}\,,\chi_{2}^{-}\,,\chi_{3}^{-}\,,\chi_{4}^{0}\right)^{T}\sim\left(1,4,-\frac{1}{2}\right),\;\langle\chi\rangle=\left(0\,,0\,,0\,,\frac{V}{\sqrt{2}}\right)^{T},$
$\displaystyle\phi$
$\displaystyle=\left(\phi_{1}^{+}\,,\phi_{2}^{0}\,,\phi_{3}^{0}\,,\phi_{4}^{+}\right)^{T}\,\sim\left(1,4,\frac{1}{2}\right),\;\langle\phi\rangle=\left(0\,,0\,,\frac{\omega}{\sqrt{2}}\,,0\right)^{T},$
$\displaystyle\rho$
$\displaystyle=\left(\rho_{1}^{+}\,,\rho_{2}^{0}\,,\rho_{3}^{0}\,,\rho_{4}^{+}\right)^{T}\sim\left(1,4,\frac{1}{2}\right),\;\langle\rho\rangle=\left(0\,,\frac{v_{1}}{\sqrt{2}}\,,0\,,0\right)^{T},$
$\displaystyle\eta$
$\displaystyle=\left(\eta_{1}^{0}\,,\eta_{2}^{-}\,,\eta_{3}^{-}\,,\eta_{4}^{0}\right)^{T}\,\sim\left(1,4,-\frac{1}{2}\right),\;\langle\eta\rangle=\left(\frac{v_{2}}{\sqrt{2}}\,,0\,,0\,,0\right)^{T}.$
(6)

The lepton masses are generated from the following Yukawa interactions

$\displaystyle-\mathcal{L}_{y}=$
$\displaystyle\;Y^{N}_{ab}\overline{L_{a}}\chi
N^{\prime}_{bR}+Y^{E}_{ab}\overline{L_{a}}\phi
E^{\prime}_{bR}+Y^{e}_{ab}\overline{L_{a}}\rho
e^{\prime}_{bR}+Y^{\nu}_{ab}\overline{L_{a}}\eta\nu^{\prime}_{bR}+\mathrm{H.c.}.$
(7)

The model consists of quark multiplets that must be arranged to cancel the
gauge anomalies, see for example a discussion in Ref. Long:2016lmj . It can be
seen that the quark masses can be constructed to satisfy the recent
experimental data. This sector is irrelevant to our work. The zero VEV values
of some neutral Higgs components can be explained by considering a global
symmetry called the general lepton number ${\cal L}$ defined as follows:

$L=\frac{4}{\sqrt{3}}\left(T_{8}+\frac{1}{\sqrt{2}}T_{15}\right)+{\cal L}\,,$
(8)

where $L$ is the normal lepton number. This formula is an extension of that
introduced for a 3-3-1 model Chang:2006aa . As a consequence, all singlets
have $L(\mathrm{singlet})={\cal L}(\mathrm{singlet)}$, namely
$\mathcal{L}(u_{aR})=\mathcal{L}(d_{aR})=0$,
$\mathcal{L}(e_{aR})=\mathcal{L}(\nu^{\prime}_{aR})=1=-\mathcal{L}(E_{aR})=-\mathcal{L}(N^{\prime}_{aR})$.
The normal lepton $L$ for a quadruplet is computed as follows

$L=\textrm{diag}\left(1+{\cal L}\,,1+{\cal L}\,,-1+{\cal L}\,,-1+{\cal
L}\right).$ (9)

In this work, we adopt only the Yukawa couplings respecting the generalized
lepton number $\mathcal{L}$ including the Lagrangian (7) for leptons. The
particular values of $\mathcal{L}$ of all multiplets are listed in Tables 1.

Multiplet |  $\chi$ |  $\phi$ |  $\eta$ |  $\rho$ |  $L_{a}$ |  $\nu^{\prime}_{aR}$ |  $e^{\prime}_{aR}$ |  $E^{\prime}_{aR}$ | $N^{\prime}_{aR}$ |  $Q_{\alpha L}$ |  $Q_{3L}$ |  $u_{aR}$ |  $d_{aR}$ |  $D_{\alpha R}$ |  $U_{\alpha R}$ |  $U_{3R}$ |  $D_{3R}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
$\cal L$ charge |  $1$ |  $1$ |  $-1$ |  $-1$ |  $0$ |  1 |  1 |  -1 | -1 |  $1$ |  $-1$ |  $0$ |  $0$ |  $2$ |  $2$ |  $-2$ |  $-2$
Table 1:  ${\cal L}$ charges for multiplets in the 341RH.

They result in consistent values of the normal lepton numbers for all SM
leptons $L(\nu^{\prime}_{aL,R})=L(e^{\prime}_{aL,R})=1$, and all SM quarks
have $L=0$. The remaining non-zero lepton numbers $L$ are listed in Table 2,

Fields |  $\chi_{1}^{0}$ |  $\chi_{2}^{-}$ |  $\phi^{+}_{1}$ |  $\phi_{2}^{0}$ |  $\rho^{0}_{3}$ |  $\rho_{4}^{+}$ |  $\eta_{3}^{-}$ |  $\eta_{4}^{0}$ |  $E^{\prime}_{aL,R}$ |  $N^{\prime}_{aL,R}$ |  $U_{\alpha L,R}$ |  $D_{\alpha L,R}$ |  $U_{3L,R}$ |  $D_{3L,R}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
$L$ |  $2$ |  $2$ |  $2$ |  $2$ |  $-2$ |  $-2$ |  $-2$ |  $-2$ |  -1 |  -1 |  2 |  2 |  -2 |  -2
Table 2:  Nonzero lepton number $L$ of all fields in the 341RH.

because of the requirement that the total lepton number $L$ is always
conservative.

The mass terms of all leptons are:

$\displaystyle(M_{\nu})_{ab}=Y^{\nu}_{ab}\frac{v_{2}}{\sqrt{2}}\,,\;(M_{e})_{ab}=Y^{e}_{ab}\frac{v_{1}}{\sqrt{2}}\,,\;(M_{E})_{ab}=Y^{E}_{ab}\frac{\omega}{\sqrt{2}},\;(M_{N})_{ab}=Y^{N}_{ab}\frac{V}{\sqrt{2}}.$
(10)

The active Dirac neutrino masses and mixing are constructed from the mass
matrix $M_{\nu}$. But, these tiny masses do not affect significantly the one-
loop contributions to $a_{e_{a}}$.

Now we focus on the lepton sector in the Yukawa part of Eq. (7). The normal
lepton mass matrix $M_{e}$ given in Eq. (10) is assumed to be diagonal for
simplicity. As a result, the flavor basis of the charged leptons
$e^{\prime}_{a}$ is the mass basis $e_{aL,R}\equiv e^{\prime}_{aL,R}$, namely

$m_{e_{a}}=Y^{e}_{ab}\delta_{ab}\frac{v_{1}}{\sqrt{2}}\Rightarrow
Y^{e}_{ab}=\delta_{ab}\frac{\sqrt{2}m_{e_{a}}}{v_{1}}.$ (11)

Three other base
$f^{\prime}_{L,R}\equiv(f^{\prime}_{1},f^{\prime}_{2},f^{\prime}_{3})^{T}_{L,R}$
with $f=\nu,E,N$ are transformed into the corresponding mass base
$f^{\prime}_{L,R}$ through the following relations:

$\displaystyle U^{f\dagger}_{L}M_{\nu}U^{f}_{R}$
$\displaystyle=\hat{M}_{f}=\mathrm{diag}(m_{f_{1}},m_{f_{2}},m_{f_{3}}),\;f^{\prime}_{L,R}=U^{f}_{L,R}f_{L,R}.$
(12)

###  II.2  Gauge boson masses and mixing

Gauge boson masses arise from the covariant kinetic term of Higgs multiplets,
namely

$L_{\mathrm{Higgs}}=\sum_{H}^{4}\left(D^{\mu}\langle
H\rangle\right)^{\dagger}D_{\mu}\langle H\rangle,$ (13)

where $H=\chi,\phi,\eta,\rho$. The covariant derivative is defined as

$\displaystyle D_{\mu}$
$\displaystyle=\partial_{\mu}-ig\sum_{a=1}^{15}W_{a\mu}T_{a}-ig_{X}XB^{\prime\prime}_{\mu}T_{16}\equiv\partial_{\mu}-igP_{\mu}^{NC}-igP_{\mu}^{CC}\,,$
(14)

where $g,g_{X}$ and $W_{a\mu},B^{\prime\prime}_{\mu}$ are gauge couplings and
fields of the gauge groups $SU(4)_{L}$ and $U(1)_{X}$, respectively. and the
two parts $P_{\mu}^{NC}$ and $P_{\mu}^{CC}$ relating with the neutral and non-
hermitian currents Long:2016lmj . For quadruplet,
$T_{16}=\frac{1}{2\sqrt{2}}\textrm{diag}(1,1,1,1)$, and

$\displaystyle P_{\mu}^{CC}=$
$\displaystyle\frac{1}{\sqrt{2}}\left(\begin{array}[]{cccc}0&W^{+}&W_{13}^{+}&W_{14}^{0}\\\
W^{-}&0&W_{23}^{0}&W_{24}^{-}\\\ W_{13}^{-}&W_{23}^{0*}&0&W_{34}^{-}\\\
W_{14}^{0*}&W_{24}^{+}&W_{34}^{+}&0\end{array}\,\right)_{\mu}\,,$ (19)

where $t\equiv g_{X}/g$ and $\sqrt{2}\,W^{\mu}_{ij}\equiv
W^{\mu}_{i}-iW^{\mu}_{j}$ with $i<j$. The upper subscripts label the electric
charges of gauge bosons. The relation between the original basis
$(W_{3},\;W_{8},W_{15},B^{\prime\prime})$ and the mass basis
$(A,Z,Z_{3},Z_{4})$ of all real neutral boson was determined previously
Long:2016lmj . The masses of non-Hermitian (charged) gauge bosons are given by

$\displaystyle m^{2}_{W}$
$\displaystyle=\frac{g^{2}(v_{1}^{2}+v_{2}^{2})}{4}\,,m^{2}_{W_{13}}=\frac{g^{2}(v_{2}^{2}+\omega^{2})}{4}\,,m^{2}_{W_{23}}=\frac{g^{2}(v_{1}^{2}+\omega^{2})}{4}\,,$
$\displaystyle m^{2}_{W_{14}}$
$\displaystyle=\frac{g^{2}(v_{2}^{2}+V^{2})}{4}\,,m^{2}_{W_{24}}=\frac{g^{2}(v_{1}^{2}+V^{2})}{4}\,,m^{2}_{W_{34}}=\frac{g^{2}(\omega^{2}+V^{2})}{4}\,.$
(20)

By spontaneous symmetry breaking (SSB), the following relation should be in
order: $V\gg\omega\gg v_{1},v_{2}$. A consequence from Eq. (20) is that
$W^{\pm}$ must be identified with the singly charged SM gauge boson, namely

$v_{1}^{2}+v_{2}^{2}=v^{2}=246^{2}\;\textrm{GeV}^{2}.$ (21)

Then neutral gauge boson masses are

$\displaystyle m^{2}_{Z}$
$\displaystyle\simeq\frac{m^{2}_{W}}{c^{2}_{W}},\;m^{2}_{Z_{3}}\simeq,\;m^{2}_{Z_{4}}\simeq$
(22)

The above formula implies that $\eta$ and $\rho$ play roles as two Higgs
doublets in the well-known two Higgs doublet models after the breaking steps
to the SM gauge group $SU(2)_{L}\times U(1)_{Y}$. Then we define the mixing
angle $\beta$ as follows

$t_{\beta}\equiv\tan\beta=\frac{v_{2}}{v_{1}},\;v_{1}=vc_{\beta},\;v_{2}=vs_{\beta},$
(23)

where $t_{\beta}\geq 0.4$ as the perturbative limit of the top quark Yukawa
coupling $|Y^{u}_{33}|\simeq\sqrt{2}m_{t}/(vs_{\beta})<\sqrt{4\pi}$. The upper
bound of $t_{\beta}$ may comes from the tau mass
$Y^{\tau}\simeq\sqrt{2}m_{\tau}/(vc_{\beta})<\sqrt{4\pi}$, equivalently
$t_{\beta}<390$.

###  II.3 Higgs potential and Higgs spectrum

The most general Higgs potential including the appearance of the singly
charged Higgs boson $\sigma^{\pm}\sim(1,1,\pm 1)$ is:

$\displaystyle V_{h}=$
$\displaystyle\;V(\eta,\rho,\phi,\chi)+V(\sigma^{\pm}),$ $\displaystyle
V(\eta,\rho,\phi,\chi)=$
$\displaystyle\;\mu^{2}_{1}\eta^{\dagger}\eta+\mu^{2}_{2}\rho^{\dagger}\rho+\mu^{2}_{3}\phi^{\dagger}\phi+\mu^{2}_{4}\chi^{\dagger}\chi$
$\displaystyle+\lambda_{1}(\eta^{\dagger}\eta)^{2}+\lambda_{2}(\rho^{\dagger}\rho)^{2}+\lambda_{3}(\phi^{\dagger}\phi)^{2}+\lambda_{4}(\chi^{\dagger}\chi)^{2}$
$\displaystyle+(\eta^{\dagger}\eta)[\lambda_{5}(\rho^{\dagger}\rho)+\lambda_{6}(\phi^{\dagger}\phi)+\lambda_{7}(\chi^{\dagger}\chi)]$
$\displaystyle+\
(\rho^{\dagger}\rho)[\lambda_{8}(\phi^{\dagger}\phi)+\lambda_{9}(\chi^{\dagger}\chi)]+\lambda^{\prime}_{9}(\phi^{\dagger}\phi)(\chi^{\dagger}\chi)$
$\displaystyle+\lambda_{10}(\rho^{\dagger}\eta)(\eta^{\dagger}\rho)+\lambda_{11}(\rho^{\dagger}\phi)(\phi^{\dagger}\rho)+\lambda_{12}(\rho^{\dagger}\chi)(\chi^{\dagger}\rho)$
$\displaystyle+\lambda_{13}(\phi^{\dagger}\eta)(\eta^{\dagger}\phi)+\lambda_{14}\,(\chi^{\dagger}\eta)(\eta^{\dagger}\chi)+\lambda_{15}\,(\chi^{\dagger}\phi)(\phi^{\dagger}\chi)$
$\displaystyle+(f\epsilon^{ijkl}\eta_{i}\rho_{j}\phi_{k}\chi_{l}+\mathrm{H.c.}),$
(24) $\displaystyle V(\sigma^{\pm})=$
$\displaystyle\;\mu_{5}^{2}\sigma^{+}\sigma^{-}+\lambda_{\sigma}\left(\sigma^{+}\sigma^{-}\right)^{2}+\left(f_{\sigma}\sigma^{+}\rho^{\dagger}\eta+\mathrm{H.c.}\right)$
$\displaystyle+(\sigma^{+}\sigma^{-})\left(\lambda_{\sigma\eta}\eta^{\dagger}\eta+\lambda_{\sigma\rho}\rho^{\dagger}\rho+\lambda_{\sigma\phi}\phi^{\dagger}\phi+\lambda_{\sigma\chi}\chi^{\dagger}\chi\right).$
(25)

Here we consider the Higgs potential respecting the generalized lepton number
$\mathcal{L}$ and the new singly charged Higgs boson has
$\mathcal{L}(\sigma^{\pm})=0$. We also ignore the triple Higgs self couplings
$(f^{\prime}_{\sigma}\chi^{\dagger}\phi+\mathrm{h.c.})$ because it results in
unnecessary mixing in our calculation. The detailed discussion to derive the
masses and mixing parameters of the Higgs bosons were presented previously in
Ref. Long:2016lmj  without $\sigma^{\pm}$. The Higgs potential consisting of
$\sigma^{\pm}$ was discussed in 3-3-1 models Hue:2021zyw ; Hue:2021xzl , in
which the detailed calculation to derive Higgs masses and mixing were
performed. We collect only the most important results relating to our work.

The mixing of $\chi_{2}^{\pm}$ and $\rho^{\pm}_{4}$ results in two massless
Goldstone boson $G^{\pm}_{24}$ and two singly charged Higgs bosons
$h^{\pm}_{3}$:

$\displaystyle\begin{pmatrix}\chi_{2}^{\pm}\\\
\rho_{4}^{\pm}\end{pmatrix}=\left(\begin{array}[]{cc}c_{\theta_{1}}&s_{\theta_{1}}\\\
-s_{\theta_{1}}&c_{\theta_{1}}\\\
\end{array}\right)\begin{pmatrix}G_{24}^{\pm}\\\
h_{3}^{\pm}\end{pmatrix},\;M^{2}_{h^{+}_{3}}=\left(c_{\beta}{}^{2}v^{2}+V^{2}\right)\left(\frac{\lambda_{12}}{2}-\frac{ft_{\beta}w}{2V}\right),$
(28)

and

$s_{\theta 1}\equiv\sin\theta_{1},\;c_{\theta
1}\equiv\cos\theta_{1},\;\tan\theta_{1}=\frac{vc_{\beta}}{V}.$ (29)

We consider here the mixing of $\phi^{0*}_{2}$ and $\rho^{0}_{3}$, which leads
to a non-hermitian Goldstone boson $G^{0}_{23}\neq G^{0*}_{23}$ and a physical
state $H^{0}_{1}\neq H^{0*}_{1}$

$\displaystyle\begin{pmatrix}\phi_{2}^{0*}\\\
\rho_{3}^{0}\end{pmatrix}=\left(\begin{array}[]{cc}c_{\theta_{2}}&s_{\theta_{2}}\\\
-s_{\theta_{2}}&c_{\theta_{2}}\\\
\end{array}\right)\begin{pmatrix}G_{23}^{0}\\\
H_{1}^{0}\end{pmatrix},\;M^{2}_{H^{0}_{1}}=\left(c_{\beta}{}^{2}v^{2}+w^{2}\right)\left(\frac{\lambda_{11}}{2}-\frac{ft_{\beta}V}{2w}\right),$
(32)

and

$s_{\theta_{2}}\equiv\sin\theta_{2},\;c_{\theta_{2}}\equiv\cos\theta_{2},\;\tan\theta_{2}=\frac{vc_{\beta}}{w}.$
(33)

Three singly charged Higgs bosons
$(\rho^{\pm}_{1},\;\eta^{\pm}_{2},\sigma^{\pm})$ are changed into the two
physical states $h^{\pm}_{1,2}$ and the Goldstone bosons $G^{\pm}_{W}$ of
$W^{\pm}$ as follows:

$\displaystyle\rho^{\pm}_{1}$
$\displaystyle=c_{\beta}\phi^{\pm}_{W}+s_{\beta}\left(c_{\alpha}h^{\pm}_{1}+s_{\alpha}h^{\pm}_{2}\right),$
$\displaystyle\eta^{\pm}_{2}$
$\displaystyle=-s_{\beta}\phi^{\pm}_{W}+c_{\beta}\left(c_{\alpha}h^{\pm}_{1}+s_{\alpha}h^{\pm}_{2}\right),$
$\displaystyle\sigma^{\pm}$
$\displaystyle=-s_{\alpha}h^{\pm}_{1}+c_{\alpha}h^{\pm}_{2}.$ (34)

The relations in Eq. (34) were given in Refs. Hue:2021zyw ; Hue:2021xzl ,
which consist of the same part $V(\sigma)$ given in Eq. (25) in the Higgs
potential. The mixing parameter $\alpha$ and Higgs boson masses
$m_{h^{+}_{1}}$, $m_{h^{+}_{2}}$ will be investigated as free parameters,
while three dependent parameters are:

$\displaystyle\mu_{5}^{2}$
$\displaystyle=c_{\alpha}{}^{2}m_{h_{2}^{+}}^{2}+s_{\alpha}{}^{2}m_{h_{1}^{+}}^{2}-\frac{1}{2}\left(c_{\beta}{}^{2}\lambda_{\sigma\rho}v^{2}+\lambda_{\sigma\chi}V^{2}+\lambda_{\sigma\phi}w^{2}+\lambda_{\sigma\eta}s_{\beta}{}^{2}v^{2}\right),$
$\displaystyle f$
$\displaystyle=-\frac{c_{\beta}s_{\beta}\left(2c_{\alpha}{}^{2}m_{h_{1}^{+}}^{2}-\lambda_{10}v^{2}+2s_{\alpha}{}^{2}m_{h_{2}^{+}}^{2}\right)}{Vw},$
$\displaystyle f_{\sigma}$
$\displaystyle=-\frac{\sqrt{2}c_{\alpha}s_{\alpha}(m_{h_{1}^{+}}^{2}-m_{h_{2}^{+}}^{2})}{v}.$
(35)

###  II.4 Analytic formulas for one-loop contributions to $a_{e_{a}}$ with
active Dirac neutrinos

In this section, we consider the simple case that $\nu_{1,2,3}$ are Dirac
ones, and no mixing between $\sigma^{\pm}$ and other singly charged Higgs
bosons, i.e. $s_{\alpha}=0$, $c_{\alpha}=1$, then only $h_{2}^{\pm}$ couples
with active neutrinos. The conclusion with $s_{\alpha}\neq 0$ is unchanged
because the tiny neutrino masses result in suppressed one-loop contributions
to $(g-2)_{e_{a}}$ anomalies.

The relevant Lagrangian giving one-loop contributions to $a_{e_{a}}$ is:

$\displaystyle-\mathcal{L}^{l}_{y}$
$\displaystyle=\frac{g}{\sqrt{2}m_{W_{23}}}\sum_{a,c=1}^{3}U_{L,ac}^{E*}\left[m_{E_{c}}t_{\theta_{2}}P_{L}+m_{e_{a}}t^{-1}_{\theta_{2}}P_{R}\right]e_{a}H_{1}^{0}$
$\displaystyle+\frac{g}{\sqrt{2}m_{W_{24}}}\sum_{a,c=1}^{3}U_{L,ac}^{N*}\left[m_{N_{c}}t_{\theta_{1}}P_{L}+m_{e_{a}}t^{-1}_{\theta_{1}}P_{R}\right]e_{a}h_{3}^{+}$
$\displaystyle+\frac{g}{\sqrt{2}m_{W}}\sum_{a,c=1}^{3}\overline{\nu_{c}}U^{\nu*}_{L,ac}\left(m_{\nu_{c}}t^{-1}_{\beta}P_{L}+m_{e_{a}}t_{\beta}P_{R}\right)e_{a}h^{+}_{1}+\mathrm{H.c.}+\dots.$
(36)

To avoid large cLFV decays $e_{b}\to e_{a}\gamma$, which may be ruled out by
experiments, we will pay attention to the limit that
$U^{N}_{L}=U^{E}_{L}=I_{3}$, and $m_{E_{a}}=m_{E}$, $m_{N_{a}}=m_{N}$ for all
$a=1,2,3$. For active neutrinos having tiny masses, the respective one-loop
contributions to $\Delta a_{e_{a}}$ are suppressed. The non-zero form factors
relevant with one-loop contributions to $a_{e_{a}}$ are Crivellin:2018qmi

$\displaystyle a_{e_{a}}(h_{3}^{+})$
$\displaystyle=-\frac{g^{2}m^{2}_{e_{a}}}{8\pi^{2}m^{2}_{W}}\left\\{\frac{v^{2}}{V^{2}}\times
x_{N}f_{\Phi}(x_{N})+\left[\frac{v^{4}c^{2}_{\beta}}{V^{4}}x_{N}+\frac{m^{2}_{e_{a}}c^{2}_{\theta_{1}}}{m^{2}_{h^{\pm}_{3}}c^{2}_{\beta}}\right]\tilde{f}_{\Phi}(x_{N})\right\\},$
$\displaystyle a_{e_{a}}(h_{1}^{+})$
$\displaystyle=-\frac{g^{2}m^{2}_{e_{a}}}{8\pi^{2}m^{2}_{W}}\left(\frac{m_{e_{a}}t_{\beta}}{m_{h^{+}_{1}}}\right)^{2}\tilde{f}_{\Phi}(0),$
(37) $\displaystyle a_{e_{a}}(H_{1}^{0})$
$\displaystyle=-\frac{g^{2}m^{2}_{e_{a}}}{8\pi^{2}m^{2}_{W}}\left\\{\frac{v^{2}}{w^{2}}x_{E}\left[f_{\Phi}(x_{E})-g_{\Phi}(x_{E})\right]+\left(\frac{v^{4}s^{2}_{\beta}}{w^{4}}x_{E}+\frac{m^{2}_{e_{a}}c^{2}_{\theta_{2}}}{m^{2}_{H^{0}_{1}}c^{2}_{\beta}}\right)\left[\tilde{f}_{\Phi}(x_{E})-\tilde{g}_{\Phi}(x_{E})\right]\right\\},$

where $x_{N}=m^{2}_{N}/m^{2}_{h^{+}_{3}}$ and
$x_{E}=m^{2}_{E}/m^{2}_{H^{0}_{1}}$, and

$\displaystyle f_{\Phi}(x)$
$\displaystyle=2\tilde{g}_{\Phi}(x)=\frac{x^{2}-1-2x\ln x}{4(x-1)^{3}},$
$\displaystyle g_{\Phi}$ $\displaystyle=\frac{x-1-\ln x}{2(x-1)^{2}},$
$\displaystyle\tilde{f}_{\Phi}(x)$
$\displaystyle=\frac{2x^{3}+3x^{2}-6x+1-6x^{2}\ln x}{24(x-1)^{4}}.$ (38)

The total contribution to $\Delta a_{e_{a}}$ from all Higgs bosons are:

$\displaystyle a_{e_{a}}(H)$ $\displaystyle=\sum_{X}a_{e_{a}}(X),$ (39)

where $X=h^{+}_{1},h^{+}_{3},H^{0}_{1}$.

The functions relating to one-loop contributions of Higgs bosons are shown in
Fig. 1.

Figure 1:  Properties of the scalar functions relating to one-loop
contributions of heavy Higgs bosons to $(g-2)_{e_{a}}$ anomalies in the 3-4-1
model with active Dirac neutrinos.

Because $xf_{\Phi}(x)$, $\tilde{f}_{\Phi}(x)$, $x\tilde{f}_{\Phi}(x)>0$ for
all $x>0$, the one-loop contributions from singly charged Higgs bosons
$H^{\pm}_{1,2}$ are in opposite signs with $\Delta a^{\mathrm{NP}}_{\mu}$,
hence they should be small. On the other hand, the one-loop contribution from
the neutral Higgs boson $H^{0}_{1}$ consists of negative functions
$\left(f_{\Phi}(x)-g_{\Phi}(x)\right)$, and
$\left(\tilde{f}_{\Phi}(x)-\tilde{g}_{\Phi}(x)\right)$, hence the final
contributions support large $\Delta a_{\mu}$.

The couplings of leptons and non-hermitian gauge bosons are

$\displaystyle\mathcal{L}^{\ell\ell V}=$ $\displaystyle
i\overline{L_{aL}}\gamma^{\mu}D_{\mu}L_{aL}$ $\displaystyle=$
$\displaystyle\frac{g}{\sqrt{2}}\left[U^{\nu*}_{ai}\overline{\nu_{i}}\gamma^{\mu}P_{L}e_{a}W^{+}_{\mu}+U^{N*}_{ai}\overline{N_{i}}\gamma^{\mu}P_{L}e_{a}W^{+}_{24\mu}+\overline{e_{aL}}\gamma^{\mu}E^{\prime}_{aL}W^{0}_{23\mu}+\dots,\right]+\mathrm{H.c.},$
(40)

where the last line collects only terms giving one-loop contributions to cLFV
amplitudes and $(g-2)$ anomalies, resulting in the following formulas
Crivellin:2018qmi ,

$\displaystyle a_{e_{a}}(W)$
$\displaystyle=-\frac{g^{2}m^{2}_{e_{a}}}{8\pi^{2}m_{W}^{2}}\sum_{i=1}^{3}U_{ai}^{\nu}U_{ai}^{\nu*}{\tilde{f}_{V}\left({x_{\nu_{i}}}\right)}\simeq-\frac{g^{2}m^{2}_{e_{a}}}{8\pi^{2}m_{W}^{2}}{\tilde{f}_{V}\left(0\right)},$
(41) $\displaystyle a_{e_{a}}(W_{24})$
$\displaystyle=-\frac{g^{2}m^{2}_{e_{a}}}{8\pi^{2}m_{W}^{2}}\sum_{b=1}^{3}U_{ab}^{N}U_{ab}^{N*}\frac{m_{W}^{2}}{m^{2}_{W_{24}}}{\tilde{f}_{V}\left({x_{N_{b}}}\right)}=-\frac{g^{2}m^{2}_{e_{a}}}{8\pi^{2}m_{W}^{2}}\times\frac{m_{W}^{2}}{m^{2}_{W_{24}}}{\tilde{f}_{V}\left({x_{N}}\right)},$
(42) $\displaystyle a_{e_{a}}(W_{23})$
$\displaystyle=-\frac{g^{2}m^{2}_{e_{a}}}{8\pi^{2}m_{W}^{2}}\sum_{b=1}^{3}U_{ab}^{E}U_{ab}^{E*}\frac{m_{W}^{2}}{m^{2}_{W_{23}}}{\tilde{f}_{V}\left({x_{E_{b}}}\right)}=-\frac{g^{2}m^{2}_{e_{a}}}{8\pi^{2}m_{W}^{2}}\times\frac{m_{W}^{2}}{m^{2}_{W_{23}}}{\tilde{f}_{V}\left({x_{E}}\right)},$
(43)

where

$\tilde{f}_{V}(x)=\frac{-4x^{4}+49x^{3}-78x^{2}+43x-10-18x^{3}\ln
x}{24(x-1)^{4}}.$ (44)

The total one contribution from new heavy gauge bosons to the $a_{e_{a}}$ is

$\displaystyle a_{e_{a}}(V)=a_{e_{a}}(W_{23})+a_{e_{a}}(W_{24}).$ (45)

The deviation of $a_{\mu}$ between predictions by the two models 3-4-1 and SM
are

$\displaystyle\Delta a^{\mathrm{341}}_{e_{a}}$ $\displaystyle\equiv\Delta
a_{e_{a}}=\Delta a_{e_{a}}({W})+a_{e_{a}}(V)+a_{e_{a}}(H),\quad\Delta
a^{W}_{e_{a}}=a^{W}_{e_{a}}-a^{\mathrm{SM}}_{e_{a}}(W),$ (46)

where $a^{\mathrm{SM}}_{\mu}(W)=3.887\times 10^{-9}$ Jegerlehner:2009ry , and
$a^{\mathrm{SM}}_{e_{a}}(W)=a^{\mathrm{SM}}_{\mu}(W)\times(m_{e_{a}}^{2}/m^{2}_{\mu})$
are the SM’s prediction for the one-loop contribution from $W$ boson. In this
work, $\Delta a^{\mathrm{341}}_{e_{a}}=\Delta a_{e_{a}}$ will be considered as
new physics (NP) predicted by the 3-4-1 models, which must satisfy the
experimental data given in Eqs. (1) and (2) in the numerical investigation.

Numerical illustrations are shown in Fig. 1 with fixed $w=5$ TeV, $V=10$ TeV,
$m_{H^{\pm}_{1}}=m_{H^{\pm}_{2}}=m_{H^{0}_{1}}=1$ TeV, and $t_{\beta}=50$.
Different contributions are considered as functions of $m_{N}=m_{E}=m_{f}$
with numerical illustrations given in Fig. 2.

|
---|---
Figure 2:  One-loop contributions of new $SU(4)_{L}$ particles to $a_{\mu}$
with $w=2.$ TeV (left panel) and $w=4$ TeV (right panel).

Both $a_{e_{a}}(W)$ and $a_{e_{a}}(H^{\pm}_{1})$ are independent with $m_{f}$,
$w$, and $V$. In addition, the one-loop contribution from $W^{\pm}$ gauge
boson is $a_{e_{a}}(W)=a^{\mathrm{SM}}_{e_{a}}(W)$, hence it does not affect
$\Delta a_{\mu}$. The two conditions $t_{\beta}\leq 100$ and
$m_{H^{\pm}_{1}}\geq 200$ GeV give $0<a_{\mu}(H_{1}^{+})\leq 10^{-12}$ hence
$a_{\mu}(H_{1}^{+})$ gives suppressed contributions to $\Delta a_{\mu}$. The
$\Delta a_{\mu}$ depends strongly on $w$, which lower bound is constrained
strictly from the heavy neutral gauge bosons $Z_{3,4}$, namely $m_{Z_{3,4}}>3$
TeV from LHC searches Lee:2014kna . As a result, large values of $w\geq 2$ TeV
gives $\Delta a_{\mu}\leq 7\times 10^{-13}\ll 192\times 10^{-13}\leq
a^{\mathrm{NP}}_{\mu}$. In conclusion, all the one-loop contributions
mentioned above are much smaller than $a^{\mathrm{NP}}_{\mu}$.

##  III  The 341ISS with ISS neutrinos

Now we consider an extension of the above 341RHN model that may explain
successfully the $(g-2)_{e_{a}}$ data. This version consists of six right-
handed neutrinos $\nu_{aR},X_{aR}\sim(1,0)$, $a=1,2,3$ generating active
neutrino masses through the ISS mechanism, and a singly charged Higgs boson
$\sigma^{+}\sim(1,1,1)$ needed to give large one-loop contributions to AMM. We
call this model the 3-4-1 model with ISS neutrinos (341ISS). The generalized
lepton numbers are $\mathcal{L}(\nu_{aR})=1$, $\mathcal{L}(X_{aR})=-1$, and
$\mathcal{L}(\sigma^{+})=0$. Now tree-level neutrino masses and mixing angles
arise from the ISS mechanism. Requiring that $\mathcal{L}$ is only softly
broken, the additional Yukawa part is

$\displaystyle-\mathcal{L}_{Y,\nu_{R}}=$ $\displaystyle
Y^{\nu}_{ab}\overline{\nu_{aR}}\eta^{\dagger}L_{b}+(M_{R})_{ab}\overline{\nu_{aR}}(X_{bR})^{c}+\frac{1}{2}(\mu_{X})_{ab}\overline{X_{aR}}(X_{bR})^{c}$
$\displaystyle+Y^{\sigma}_{ab}\overline{(X_{aR})^{c}}e_{bR}\sigma^{+}+\mathrm{h.c.},$
(47)

where $Y^{\nu}$, $M_{R}$, $\mu_{X}$, and $Y^{\sigma}$ are $3\times 3$
matrices. The first term in Eq. (47) is similar to that given in Eq. (7), but
it will generate the Dirac neutrino mass matrix $M_{D}$ instead of the active
neutrino masses.

Notations for flavor states of active left-handed neutrinos are
$\nu_{L}=(\nu^{\prime}_{1},\nu^{\prime}_{2},\nu^{\prime}_{3})^{T}_{L}$ and
$\nu_{R}=(\nu_{1},\nu_{2},\nu_{3})^{T}_{R}$,
$X_{R}=(X_{1},X_{2},X_{3})^{T}_{R}$, the neutrino mass terms derived from (47)
are written in the following ISS form ParticleDataGroup:2020ssz :

$\displaystyle-\mathcal{L}^{\nu}_{\mathrm{mass}}$
$\displaystyle\equiv\frac{1}{2}\left(\overline{(\nu_{L})^{c}},\;\overline{\nu_{R}},\;\overline{X_{R}}\right)\mathcal{M}^{\nu}\begin{pmatrix}\nu_{L}\\\
\left(\nu_{R}\right)^{c}\\\
\left(X_{R}\right)^{c}\end{pmatrix}+\mathrm{h.c.},$
$\displaystyle\mathcal{M}^{\nu}$
$\displaystyle=\left(\begin{array}[]{cc}\mathcal{O}_{3\times 3}&M_{D}^{T}\\\
M_{D}&M_{N}\\\ \end{array}\right),\;M_{D}=\begin{pmatrix}m_{D}\\\
\mathcal{O}_{3\times
3}\end{pmatrix},\;M_{N}=\left(\begin{array}[]{cc}\mathcal{O}_{3\times
3}&M_{R}\\\ M^{T}_{R}&\mu_{X}\\\ \end{array}\right),$ (52)

where $\mathcal{O}_{3\times 3}$ is a zero matrix and $m_{D}=Y^{\nu}\times
v_{2}/\sqrt{2}$. The analytic form of the Dirac mass matrix was chosen
generally following Ref. Casas:2001sr .

The total unitary mixing matrix $U^{\nu}$ is defined as follows

$\displaystyle U^{\nu
T}\mathcal{M}^{\nu}U^{\nu}=\hat{\mathcal{M}}^{\nu}=\mathrm{diag}(m_{n_{1}},m_{n_{2}},m_{n_{3}},m_{n_{4}},...,m_{n_{9}})\equiv\mathrm{diag}\left(\hat{m}_{\nu},\;\hat{M}_{N}\right),$
(53)

where $m_{n_{i}}$ ($i=1,2,...,9$) are eigenvalues of the $9$ mass eigenstates
$n_{iL}$, including three light active neutrinos $n_{aL}$ ($a=1,2,3$) with
mass matrix $\hat{m}_{\nu}$ and six other heavy neutrinos with mass matrix
$\hat{M}_{N}$. The relation between the flavor and mass eigenstates are

$\begin{pmatrix}\nu_{L}\\\ \left(\nu_{R}\right)^{c}\\\
\left(X_{R}\right)^{c}\end{pmatrix}=U^{\nu}n_{L},\hskip
14.22636pt\mathrm{and}\;\begin{pmatrix}\left(\nu_{L}\right)^{c}\\\ \nu_{R}\\\
X_{R}\end{pmatrix}=U^{\nu*}n_{R},$ (54)

where $n_{L}\equiv(n_{1L},n_{2L},...,n_{9L})^{T}$ and
$n_{R}=\left(n_{L}\right)^{c}$. The neutrino mixing matrix is parameterised in
the following form:

$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}),$ (55)

where $U_{\mathrm{PMNS}}$ is the $3\times 3$ Pontecorvo-Maki-Nakagawa-Sakata
(PMNS) matrix Pontecorvo:1957cp ; Maki:1962mu , $V$ is a $6\times 6$ unitary
matrix, and $R$ is a $3\times 6$ matrix satisfying $|R_{aI}|\ll 1$ for all
$a=1,2,3$ and $I=1,2,...,6$. In the ISS framework we considered here, $m_{D}$
is parameterized in terms of many free parameters, hence it is enough to
choose a simple form of $\mu_{X}=\mu_{0}I_{3}$ and
$M_{R}=\hat{M}_{R}=M_{0}I_{3}$ Arganda:2014dta ; Thao:2017qtn . The formulas
of $m_{D}$ and mixing parameters are Casas:2001sr

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

where max$[\left|\left(\hat{x}_{\nu}\right)_{aa}\right|]\ll 1$ for all
$a=1,2,3$.

The ISS conditions $|\hat{m}_{\nu}|\ll|\mu_{0}|\ll|m_{D}|\ll M_{0}$ so that
$\frac{\sqrt{\mu_{0}\hat{m}_{\nu}}}{M_{0}}\simeq 0$, the mixing matrix and
Majorana mass term are

$\displaystyle\hat{M}_{N}=\left(\begin{matrix}\hat{M}_{R}&0\\\
0&\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),$ (57)

which give $V^{*}\hat{M}_{N}V^{\dagger}\simeq M_{N}$, i.e., $m_{n_{i}}\simeq
M_{0}$ for all $i=4,5,...,9$.

###  III.1 Analytic formulas for one-loop contributions to $a_{e_{a}}$ and
Br$(e_{b}\to e_{a}\gamma)$

We note that except the contributions from ISS neutrino couplings with singly
charged Higgs bosons given in Eq. (34), all other contributions are the same
results as those discussed in the case of Dirac active neutrinos mentioned
above. Hence, we just focus here on the Higgs contributions relating to the
couplings with ISS neutrinos. The relevant couplings are listed in the
following Lagrangian

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

where

$\displaystyle\lambda^{L,1}_{ia}$
$\displaystyle=\sum_{I=1}^{6}t_{\beta}^{-1}M_{D,Ia}c_{\alpha}U^{\nu}_{(I+3)i}\simeq\frac{t_{\beta}^{-1}c_{\alpha}}{\sqrt{2}}\times\left[\begin{array}[]{cc}0,&\quad
i\leq 3\\\
-iM_{0}\left(U^{*}_{\mathrm{PMNS}}\hat{x}^{1/2}_{\nu}\right)_{a(i-3)}&\quad
3<i\leq 6\\\
M_{0}\left(U^{*}_{\mathrm{PMNS}}\hat{x}^{1/2}_{\nu}\right)_{a(i-6)}&\quad
i\geq 7\end{array}\right.,$ (62) $\displaystyle\lambda^{L,2}_{ia}$
$\displaystyle\simeq\lambda^{L,1}_{ia}t_{\alpha},$
$\displaystyle\lambda^{R,1}_{ia}$
$\displaystyle=m_{e_{a}}t_{\beta}c_{\alpha}U^{\nu*}_{ai}-\sum_{I=1}^{6}\frac{v}{\sqrt{2}}s_{\alpha}Y^{\sigma}_{Ia}U^{\nu*}_{(I+3)i}$
$\displaystyle\simeq\left[\begin{array}[]{cc}m_{e_{a}}t_{\beta}c_{\alpha}\left[U^{*}_{\mathrm{PMNS}}\left(I_{3}-\frac{1}{2}\hat{x}_{\nu}\right)\right]_{ai}+\frac{vs_{\alpha}}{\sqrt{2}}\left(Y^{\sigma
T}\hat{x}_{\nu}^{1/2}\right)_{ai}&\quad i\leq 3\\\
\frac{-i}{\sqrt{2}}m_{e_{a}}t_{\beta}c_{\alpha}\left(U^{*}_{\mathrm{PMNS}}\hat{x}_{\nu}^{1/2}\right)_{a(i-3)}+\frac{ivs_{\alpha}}{2}\left[Y^{\sigma
T}\left(I_{3}-\frac{1}{2}\hat{x}_{\nu}\right)\right]_{a(i-3)}&\quad 4\leq
i<7\\\
\frac{1}{\sqrt{2}}m_{e_{a}}t_{\beta}c_{\alpha}\left(U^{*}_{\mathrm{PMNS}}\hat{x}_{\nu}^{1/2}\right)_{a(i-6)}-\frac{vs_{\alpha}}{2}\left[Y^{\sigma
T}\left(I_{3}-\frac{1}{2}\hat{x}_{\nu}\right)\right]_{a(i-6)}&\quad i\geq
7\end{array}\right.,$ (66) $\displaystyle\lambda^{R,2}_{ia}$
$\displaystyle=\lambda^{R,1}_{ia}\left[s_{\alpha}\to-
c_{\alpha},\;c_{\alpha}\to s_{\alpha}\right].$ (67)

The branching ratios of the cLFV decays are formulated as follows
Lavoura:2003xp ; Hue:2017lak ; Crivellin:2018qmi :

$\displaystyle\mathrm{Br}(e_{b}\to
e_{a}\gamma)=\frac{48\pi^{2}}{G_{F}^{2}m_{b}^{2}}\left(\left|c_{(ab)R}\right|^{2}+\left|c_{(ba)R}\right|^{2}\right)\mathrm{Br}(e_{b}\to
e_{a}\overline{\nu_{a}}\nu_{b}),$ (68)

where $G_{F}=g^{2}/(4\sqrt{2}m_{W}^{2})$, Br$(\mu\to e\gamma)\simeq 1$,
Br$(\tau\to e\gamma)\simeq 0.1782$, Br$(\tau\to\mu\gamma)\simeq 0.1739$
ParticleDataGroup:2020ssz , and

$\displaystyle c_{(ab)R}$
$\displaystyle=\sum_{k=1}^{2}c_{(ab)R}\left(h^{\pm}\right)+c_{(ab)R}(W),\quad
c_{(ba)R}=c_{(ab)R}[a\to b,\;b\to a],$ $\displaystyle
c_{(ab)R}\left(h^{\pm}\right)$
$\displaystyle=\sum_{k=1}^{2}c_{(ab)R}\left(h^{\pm}_{k}\right),$
$\displaystyle c_{(ab)R}\left(h^{\pm}_{k}\right)$
$\displaystyle=\frac{g^{2}e\;}{32\pi^{2}m^{2}_{W}m^{2}_{h^{\pm}_{k}}}$
$\displaystyle\times\sum_{i=1}^{9}\left[\lambda^{L,k*}_{ia}\lambda^{R,k}_{ib}m_{n_{i}}f_{\Phi}(x_{i,k})+\left(m_{e_{b}}\lambda^{L,k*}_{ia}\lambda^{L,k}_{ib}+m_{e_{a}}\lambda^{R,k*}_{ia}\lambda^{R,k}_{ib}\right)\tilde{f}_{\Phi}(x_{i,k})\right],$
(69)

with $x_{i,k}\equiv m^{2}_{n_{i}}/m^{2}_{h^{\pm}_{k}}$.

Up to the order $\mathcal{O}(R^{2})$ of the neutrino mixing matrix given in
Eq. (55), the non-zero one-loop contributions relating to $h^{\pm}_{1,2}$ are

$\displaystyle c_{(ab)R}\left(h^{\pm}_{1}\right)$
$\displaystyle=\frac{eG_{F}m_{e_{b}}}{4\sqrt{2}\pi^{2}}$
$\displaystyle\times\left\\{\left[c^{2}_{\alpha}\left(U_{\mathrm{PMNS}}\hat{x}_{\nu}U_{\mathrm{PMNS}}^{\dagger}\right)_{ab}-\frac{vt_{\beta}^{-1}c_{\alpha}s_{\alpha}}{m_{e_{b}}\sqrt{2}}\left(U_{\mathrm{PMNS}}\hat{x}_{\nu}^{1/2}Y^{\sigma}\right)_{ab}\right]x_{1}f_{\Phi}(x_{1})\right.$
$\displaystyle\quad+t_{\beta}^{-2}c^{2}_{\alpha}\left(U_{\mathrm{PMNS}}\hat{x}_{\nu}U_{\mathrm{PMNS}}^{\dagger}\right)_{ab}x_{1}\tilde{f}_{\Phi}(x_{1})$
$\displaystyle\quad+\frac{m^{2}_{e_{a}}t_{\beta}^{2}c_{\alpha}^{2}}{m^{2}_{h^{\pm}_{1}}}\left[\frac{\delta_{ab}}{24}-\left(U_{\mathrm{PMNS}}\hat{x}_{\nu}U_{\mathrm{PMNS}}^{\dagger}\right)_{ab}\left(\frac{1}{24}-\tilde{f}_{\Phi}(x_{1})\right)\right]$
$\displaystyle\quad+\frac{m_{e_{a}}v^{2}s_{\alpha}^{2}}{2m_{e_{b}}m^{2}_{h^{\pm}_{1}}}\left[\left(Y^{\sigma\dagger}\hat{x}_{\nu}Y^{\sigma}\right)_{ab}\left(\frac{1}{24}-\tilde{f}_{\Phi}(x_{1})\right)+\left(Y^{\sigma\dagger}Y^{\sigma}\right)_{ab}\tilde{f}_{\Phi}(x_{1})\right]$
$\displaystyle\quad\left.+\frac{vm_{e_{a}}t_{\beta}s_{2\alpha}}{2\sqrt{2}m^{2}_{h^{\pm}_{1}}}\left(\frac{1}{24}-\tilde{f}_{\Phi}(x_{1})\right)\left[\frac{m_{e_{a}}}{m_{e_{b}}}\left(U_{\mathrm{PMNS}}\hat{x}_{\nu}^{1/2}Y^{\sigma}\right)_{ab}+\left(U_{\mathrm{PMNS}}\hat{x}_{\nu}^{1/2}Y^{\sigma}\right)^{*}_{ba}\right]\right\\},$
$\displaystyle c_{(ab)R}(h^{\pm}_{2})$
$\displaystyle=c_{(ab)R}(h^{\pm}_{1})\left[x_{1}\to\;x_{2},\;s_{\alpha}\to-
c_{\alpha},\;c_{\alpha}\to s_{\alpha}\right],$ $\displaystyle c_{(ab)R}(W)$
$\displaystyle\simeq\frac{eG_{F}m_{e_{b}}}{4\sqrt{2}\pi^{2}}\left[-\frac{5\delta_{ab}}{12}+\left(U_{\mathrm{PMNS}}\hat{x}_{\nu}U_{\mathrm{PMNS}}^{\dagger}\right)_{ab}\times\left(\tilde{f}_{V}\left(x_{W}\right)+\frac{5}{12}\right)\right].$
(70)

The one-loop contributions from the $h^{\pm}_{1,2}$ exchanges to $a_{e_{a}}$
are:

$\displaystyle a_{e_{a}}(h^{\pm}_{1})$
$\displaystyle=-\frac{4m_{e_{a}}}{e}\mathrm{Re}[c_{(aa)R}(h^{\pm}_{1})],$
$\displaystyle a_{e_{a}}(h^{\pm}_{2})$
$\displaystyle=a_{\mu}(h^{\pm}_{1})\left[x_{1}\to\;x_{2},\;s_{\alpha}\to-
c_{\alpha},\;c_{\alpha}\to s_{\alpha}\right],$ $\displaystyle
a_{e_{a}}(h^{\pm})$
$\displaystyle=a_{e_{a}}(h^{\pm}_{1})+a_{e_{a}}(h^{\pm}_{2}),$ (71)

where $x_{k}=M_{0}^{2}/m^{2}_{h^{\pm}_{k}}$ and $a_{e_{a}}(h^{\pm})$ is the
total contribution from these two Higgs bosons.

For qualitative estimation, the main contribution to $a_{e_{a}}(h^{\pm})$ is
the chirally-enhanced part Crivellin:2018qmi . From Eq. (III.1), this part is
determined as follows

$a_{e_{a},0}=\frac{G_{F}m^{2}_{e_{a}}}{\sqrt{2}\pi^{2}}\times\mathrm{Re}\left\\{\left[\frac{vt_{\beta}^{-1}c_{\alpha}s_{\alpha}}{\sqrt{2}m_{e_{a}}}U_{\mathrm{PMNS}}\hat{x}_{\nu}^{1/2}Y^{\sigma}\right]_{aa}\left[x_{1}f_{\Phi}(x_{1})-x_{2}f_{\Phi}(x_{2})\right]\right\\}.$
(72)

But this term may also give large contributions to $c_{(ab)R}$ and
$c_{(ba)R}$, which may be excluded by the cLFV constraints with $b\neq a$. To
avoid this problem, we start from the diagonal form of $c_{(ab)R,0}$ as
follows

$c_{(ab)R,0}\sim\left[U_{\mathrm{PMNS}}\hat{x}_{\nu}^{1/2}Y^{\sigma}\right]_{ab}\sim\delta_{ab}.$
(73)

Correspondingly, the formula of $a_{e_{a},0}$ is proportional to a diagonal
matrix $Y^{d}$ satisfying:

$\displaystyle
U_{\mathrm{PMNS}}\times\mathrm{diag}\left(\frac{m_{n_{1}}}{m_{n_{3}}},\;\frac{m_{n_{2}}}{m_{n_{3}}},\;1\right)^{1/2}Y^{\sigma}=Y^{d}\equiv\mathrm{diag}\left(Y^{d}_{11},\;Y^{d}_{22},Y^{d}_{33}\right),$
(74)

Hence, the diagonal entries will give main contributions to $a_{e_{a}}$,
namely

$a_{e_{a},0}=\frac{G_{F}m^{2}_{e_{a}}\sqrt{x_{0}}}{\sqrt{2}\pi^{2}}\times\mathrm{Re}\left[\frac{vt_{\beta}^{-1}c_{\alpha}s_{\alpha}}{\sqrt{2}m_{e_{a}}}Y^{d}\right]_{aa}\left[x_{1}f_{\Phi}(x_{1})-x_{2}f_{\Phi}(x_{2})\right].$
(75)

Then $a_{e_{a},0}$ may be large, provided that $t_{\beta}$ should not too
large, $t_{1}\neq t_{2}$, and the following quantities are large enough:
$x_{0}$, $|s_{2\alpha}=2s_{\alpha}c_{\alpha}|$, and $|Y^{d}_{aa}|$ with
$a=1,2,3$. In contrast, cLFV amplitudes do not get any contributions from
$c_{(ab)R,0}$. Numerical investigations will be done to check this conclusion.
For simplicity, we use the approximation that all tiny contributions are
ignored in the numerical results. Namely, $\Delta a_{e_{a}}\equiv
a_{e_{a}}(h^{\pm})$ given in Eq. (III.1), instead of the total formula given
in Eq. (46). In contrast, the one-loop contribution from $W$ must be included
in the formulas of Br$(e_{b}\to e_{a}\gamma)$. This conclusion was confirmed
based on the qualitative estimation discussed above. The numerical checks have
been performed, which are consistent with previous discussions on 3-3-1 models
Hue:2021xap ; Hong:2022xjg .

##  IV  Numerical discussion

We will use the best-fit values of the neutrino osculation 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}].$ (76)

The active mixing matrix and neutrino masses are determined as follows

$\displaystyle\hat{m}_{\nu}$
$\displaystyle=\left(\hat{m}^{2}_{\nu}\right)^{1/2}=\mathrm{diag}\left(m_{n_{1}},\;\sqrt{m^{2}_{n_{1}}+\Delta
m^{2}_{21}},\;\sqrt{m^{2}_{n_{1}}+\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}e^{i\delta}s_{13}-c_{12}s_{23}&c_{13}c_{23}\\\
\end{array}\right).$ (80)

This choice of active neutrino masses also satisfies the constraint from
Plank2018 Planck:2018vyg , $\sum_{i=a}^{3}m_{n_{a}}\leq 0.12\;\mathrm{eV}$.

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 phenomenology such as electroweak precision Fernandez-
Martinez:2016lgt ; Agostinho:2017wfs ; Coutinho:2019aiy , leading to a very
strict constraint of
$\eta\equiv\frac{1}{2}\left|RR^{\dagger}\right|\sim\hat{x}_{\nu}$ in the ISS
framework Mondal:2021vou ; Biggio:2019eeo ; Escribano:2021css . The reasonable
constraint of the non-unitary part of the neutrino mixing matrix is as follows

$x_{0}\equiv\frac{m_{n_{3}}}{\mu_{0}}\leq 10^{-3}.$ (81)

The other well-known numerical parameters are ParticleDataGroup:2020ssz

$\displaystyle g$ $\displaystyle=0.652,\;G_{F}=1.1664\times
10^{-5}\;\mathrm{GeV},\;s^{2}_{W}=0.231,\;m_{W}=80.385\;\mathrm{GeV},$
$\displaystyle m_{e}$ $\displaystyle=5\times
10^{-4}\;\mathrm{GeV},\;m_{\mu}=0.105\;\mathrm{GeV},\;m_{\tau}=1.776\;\mathrm{GeV}.$
(82)

For the free parameters of the 341ISS model, the numerical scanning ranges are

$\displaystyle M_{0}$
$\displaystyle\in\left[0.1,\;10\right]\;\mathrm{TeV},\;m_{h^{\pm}_{1,2}}\in\left[0.8,\;10\right]\;\mathrm{TeV},$
$\displaystyle t_{\beta}$
$\displaystyle\in\left[0.3,50\right],\;x_{0}\in\left[10^{-6},10^{-3}\right],\;s_{\alpha}\in\left[-1,1\right],\;|Y^{d}_{ab}|\leq
3.5\;\forall a,b=1,2,3.$ (83)

In addition, we will fix $m_{n_{1}}=0.01$ eV, and check the perturbative limit
of all Yukawa couplings $Y^{\sigma}$ and $Y^{\nu}$, namely
$|Y^{\sigma}_{ab}|,|Y^{\nu}_{ab}|\leq 3.5$ must be satisfied. As we discussed
above, the diagonal form of $Y^{d}$ will allow large $(g-2)_{e_{a}}$ anomalies
and small Br$(e_{b}\to e_{a}\gamma)$ satisfying the recent experimental
constraints. In the numerical investigation, we will consider the general case
that $Y^{d}_{ab}\neq 0$. The allowed regions of parameters we imply below must
guarantee both the experimental data of $1\sigma$ ranges of $(g-2)_{e,\mu}$
and the cLFV constraints of Br$(e_{b}\to e_{a}\gamma)$.

Now, we consider some particular different choices of zero entries of $Y^{d}$.
The allowed regions of the parameter space predict some interesting
properties. First, even with only $Y^{d}_{11},Y^{d}_{22}\neq 0$, the allowed
regions satisfying two $(g-2)_{e,\mu}$ data are still strictly constrained by
the cLFV decay rate Br$(\mu\to e\gamma)<4.2\times 10^{-13}$, namely the
constraints $|Y^{d}_{12}|,|Y^{d}_{12}|<10^{-4}$ must guarantee. Therefore, we
will fix the $|Y^{d}_{12}|,|Y^{d}_{12}|=0$ as the default values in our
numerical investigations. Second, we give comments on the three following
cases:

  1. 1.

$Y^{d}_{23}=Y^{d}_{32}=Y^{d}_{13}=Y^{d}_{31}=0$, which is the case of $Y^{d}$
being diagonal, as given in Eq. (74). The allowed regions of non-zero entries
of $Y^{d}$ are $0.02\leq|Y^{d}_{11}|\leq 0.25$, $0.72\leq|Y^{d}_{22}|\leq 3$
and there is no more strict constraint on $Y^{d}_{33}$, see the left panel of
Fig. 3.

|
---|---
Figure 3:  The correlations between $\Delta a_{\mu}(h^{\pm})$ vs $Y^{d}_{aa}$
(left panel) and $\Delta a_{e}(h^{\pm})$ vs Br$(e_{b}\to e_{a}\gamma)$ (right
panel) in the case $Y^{d}_{23}=Y^{d}_{32}=Y^{d}_{13}=Y^{d}_{31}=0$.

In this allowed region, Br$(\tau\to\mu\gamma)$ can reach the order of
$\mathcal{O}(10^{-9})$, but predicts suppressed branching ratios Br$(\tau\to
e\gamma)<10^{-11}$, see the right panel of Fig. 3. We can see that
Br$(e\to\mu\gamma)$ can reach the recent constraint of $\mathcal{O}(10^{-13})$
even when $Y^{d}_{12}=Y^{d}_{21}=0$. This property is different from two other
cLFV decays of $\tau$, which are suppressed if $Y^{d}_{33}=0$ is fixed, namely
Br$(\tau\to\mu\gamma)<10^{-11}$. In addition, there exist regions that allow
both small values of Br$(e\to\mu\gamma)\sim\mathcal{O}(10^{-14})$
corresponding to the future experimental sensitivity and large $\Delta
a_{e,\mu}$. In conclusion, we confirm that the formula of $a_{e,0}$ given in
Eq. (72) is the main one-loop contribution to $(g-2)_{e_{a}}$ data and cLFV
decay amplitudes. Hence, the diagonal form of $Y^{d}$ results in small cLFV
decay rates, even though they get contributions from other terms in Eq.
(III.1).

  2. 2.

$Y^{d}_{32}=Y^{d}_{23}=Y^{d}_{33}=0$ while $Y^{d}_{31},Y^{d}_{31}\neq 0$. The
two cLFV rates Br$(\tau\to e\gamma)$ and Br$(\mu\to e\gamma)$ can reach recent
experimental upper bound, while Br$(\tau\to\mu\gamma)<3\times 10^{-12}$. The
allowed ranges of $Y^{d}_{11}$ and $Y^{d}_{22}$ are almost unchanged. The
constraints of $Y^{d}_{13,31}$ are $|Y^{d}_{13}|<0.5$ and $|Y^{d}_{31}|<0.3$.
Illustrations for relations between entries of $Y^{d}$ and Br$(\tau\to
e\gamma)$ are presented in Fig. 4.

|
---|---
Figure 4:  The correlations between Br$(\tau\to e\gamma)$ vs $Y^{d}_{13,31}$
(left panel) and $Y^{d}_{aa}$ (right panel) in the second case
$Y^{d}_{23}=Y^{d}_{32}=Y^{d}_{33}=0$.

We conclude that Br$(\tau\to e\gamma)$ depends strongly on $Y^{d}_{13}$ and
$Y^{d}_{31}$. In contrast, the allowed ranges of $Y^{d}_{11}$ and $Y^{d}_{22}$
are mainly arising from the $(g-2)$ data.

  3. 3.

$Y^{d}_{13}=Y^{d}_{31}=Y^{d}_{33}=0$ while $Y^{d}_{32},Y^{d}_{2,3}\neq 0$. Two
cLFV decays Br$(\mu\to e\gamma)$ and Br$(\tau\to\mu\gamma)$ can reach recent
experimental bounds, but very suppressed Br$(\tau\to e\gamma)<10^{-12}$. The
constraints of the non-zero entries of $Y^{d}$ are $|Y^{d}_{23,32}|\leq 1$.
Illustrations between Br$(\tau\to\mu\gamma)$ vs. $Y^{d}_{23,32}$ and
$a_{\mu,e}(h^{\pm})$ are similar to the case 2, hence we do not show
explicitly here.

We comment here on some properties of entries of $Y^{d}$ derived from studying
three particular cases mentioned above. Firstly, the diagonal entries of
$Y^{d}$ give main contributions to $(g-2)_{e_{a}}$ anomalies, while the non-
zero entries $Y^{d}_{13,31}$ and $Y^{d}_{23,32}$ affect strongly in the decay
rates of Br$(\tau\to e\gamma)$ and Br$(\tau\to\mu\gamma)$, respectively. The
Br$(\mu\to e\gamma)$ depends strongly on $Y^{d}_{12,21}$ and the combination
of the remaining contributions. The recent experimental bound of Br$(\mu\to
e\gamma)<4.2\times 10^{-13}$ results in tiny values $|Y^{d}_{12,21}|<10^{-4}$,
hence we always fix these two entries being zeros.

In the final illustration, we consider the more general case that the only two
zero entries of $Y^{d}$ are $Y^{d}_{12}=Y^{d}_{21}=0$. The correlations of
important parameters vs. $\Delta a_{\mu}(h^{\pm})$ are shown in Fig. 5.

|
---|---
|
|
Figure 5:  The correlations between different free parameters vs $\Delta
a_{\mu}(h^{\pm})$ with $Y_{21}=Y_{12}=0$.

The corresponding allowed ranges of the free parameters are

$\displaystyle t_{\beta}$
$\displaystyle\in\left[0.35,22.59\right],\;s_{\alpha}\in\left[-0.991,-0.06\right]\cup\left[0.075,0.994\right],\;x_{0}\in\left[1.82\times
10^{-6},2.5\times 10^{-4}\right],$ $\displaystyle\;M_{0}$
$\displaystyle\in\left[0.691,10\right]\;\text{[TeV]},\;M_{1,2}\in\left[0.8,10\right]\;\text{[TeV]},$
$\displaystyle Y_{11}^{d}$
$\displaystyle\in\left[-0.243,-0.03\right]\cup\left[0.025,0.241\right],\;Y_{22}^{d}\in\left[-3.5,-1.08\right]\cup\left[0.92,3.498\right],\;|Y_{33}^{d}|\leq
3.5,$ $\displaystyle Y_{13}^{d}$
$\displaystyle\in\left[-0.84,\;0.82\right],\;Y_{31}^{d}\in\left[-0.23,\;0.21\right],\;Y_{23}^{d}\in\left[-0.94,\;0.91\right],\;Y_{32}^{d}\in\left[-1.08,\;1.13\right].$
(84)

We can see that many allowed ranges are stricter than the scanning ones given
in Eq. (IV). For example, large $(g-2)_{e,\mu}$ requires small $t^{\beta}$ and
large $x_{0}$ so that the new upper bound of $t_{\beta}$ and lower bound of
$x_{0}$ are determined. As a result, in Fig. 5, the allowed regions of small
$t_{\beta}$ and large $x_{0}$ are favored. Similarly, $a_{e_{a},0}\sim
s_{\alpha}c_{\alpha}=s_{2\alpha}/2$, therefore, the allowed values of large
$s_{2\alpha}$ close to 1 are supported. In the bottom left panel, the allowed
regions of $Y^{d}_{11,22}$ are the same as those predicted in Fig. 3, in which
all non-diagonal entries are zeros. This implies that these regions are
independent of non-diagonal entries of $Y^{d}$. Instead, large values of these
entries favor the small $a_{\mu}(h^{\pm})$, see the bottom right panel. As a
result, small $(g-2)_{\mu}$ may predict large Br$(\tau\to\mu\gamma,\mu\gamma)$
and vice versa. In addition to explain both $(g-2)_{e_{a}}$ data of $\Delta
a_{\mu}(h^{\pm})$, the condition of $|m_{h^{\pm}_{1}}-m_{h^{\pm}_{2}}|\geq
397$ GeV is required, which is condition derived from $a_{e_{a},0}$ that
$x_{1}\neq x_{2}$. We emphasize that this model allows the existence of heavy
charged Higgs bosons, which did not appear in some recent discussions on
$(g-2)$ anomalies DelleRose:2020oaa ; Botella:2020xzf .

The relation of $Y^{d}_{ij}$ vs. $\Delta a_{e}(h^{\pm})$ are shown in Fig. 6.

|
---|---
Figure 6:  The correlations between $Y^{d}_{ij}$ vs $\Delta a_{e}(h^{\pm})$
with $Y_{21}=Y_{12}=0$.

In the left panel, $a_{e}(h^{\pm})\sim|Y^{d}_{11}|$ confirms the relation
given in Eq. (75). The $a_{\mu}(h^{\pm})\sim|Y^{d}_{22}|$ is not very clear,
because many points are excluded by the perturbative limit. The right panel
shows that $a_{e}(h^{\pm})$ does not give any prediction on $Y^{d}_{ij}$ like
the case of $a_{\mu}(h^{\pm})$.

The correlations of $Y^{d}_{ij}$ vs Br$(\mu\to e\gamma)$ are shown in Fig. 7.

|
---|---
Figure 7:  The correlations between $Y^{d}_{ij}$ vs Br$(\tau\to
e\gamma,\mu\gamma)$ with $Y_{21}=Y_{12}=0$.

We see again that the constraints of $Y^{d}_{ii}$ with $i=1,2$ are similar to
those shown in the left panel of Fig. 3. This confirms the conclusion that the
allowed ranges of $Y^{d}_{11,22}$ are mainly controlled by the $(g-2)_{e_{a}}$
data. In the right panel of Fig. 7, large Br$(\mu\to e\gamma)$ prefers large
non-diagonal entries $Y^{d}_{ij}$ but the small values are still allowed
because of the destructive correlations between contributions from different
parameters.

The correlations between $\Delta a_{e,\mu}$ and Br$(e_{b}\to e_{a}\gamma)$ in
the allowed regions of parameters listed in Eq. (IV) are shown in Fig. 8.

|
---|---
Figure 8:  The correlations between Br$(e_{b}\to e_{a}\gamma)$ vs $\Delta
a_{e,\mu}(h^{\pm})$ with $Y_{21}=Y_{12}=0$.

We can see that all allowed values of $a_{\mu,e}(h^{\pm})$ in the $1\sigma$
experimental ranges still predict large Br$(e_{b}\to e_{a}\gamma)$ near the
current upper experimental bounds. Therefore, the future results of cLFV
experiments, $(g-2)$ data, and neutrino oscillation will give more strict
constraints on the allowed regions of the parameter space.

##  V  Conclusion

We have discussed a solution to explain the recent experimental data of the
$(g-2)_{e_{a}}$ anomalies in the 341ISS framework. We have constructed the
Yukawa Lagrangian of leptons and Higgs potential obeying the generalized
lepton number $\mathcal{L}$ that keeps necessary terms generating the ISS
mechanism and large chirally-enhanced one-loop contributions to
$(g-2)_{e_{a}}$ anomalies. Although the ISS mechanism may result in large cLFV
decays $e_{b}\to e_{a}\gamma$, we have shown numerically that there always
exist allowed regions of the parameter space guaranteeing these experimental
bounds. In addition, these allowed regions will not be excluded totally if the
future sensitivities of the cLFV experiments are updated, and no cLFV
significations are found. The model can also explain successfully the
existence of at least one of the cLFV decays $\tau\to\mu\gamma,e\gamma$ or
$\mu\to e\gamma$ once they are detected by incoming experiments.

## Acknowledgments

This research is funded by Vietnam Ministry of Education and Training and
Hanoi Pedagogical University 2 under grant number B.2021-SP2-05.

## References

  * (1)
  * (2) R. Foot, H. N. Long and T. A. Tran, Phys. Rev. D 50 (1994) no.1, R34-R38 [arXiv:hep-ph/9402243 [hep-ph]].
  * (3) F. Pisano and V. Pleitez, Phys. Rev. D 51 (1995), 3865-3869 [arXiv:hep-ph/9401272 [hep-ph]].
  * (4) M. B. Voloshin, Sov. J. Nucl. Phys. 48 (1988), 512 ITEP-87-215.
  * (5) W. A. Ponce, D. A. Gutierrez and L. A. Sanchez, Phys. Rev. D 69 (2004), 055007 [arXiv:hep-ph/0312143 [hep-ph]].
  * (6) L. A. Sanchez, F. A. Perez and W. A. Ponce, Eur. Phys. J. C 35 (2004), 259-265 [arXiv:hep-ph/0404005 [hep-ph]].
  * (7) W. A. Ponce and L. A. Sanchez, Mod. Phys. Lett. A 22 (2007), 435-448 [arXiv:hep-ph/0607175 [hep-ph]].
  * (8) L. A. Sanchez, L. A. Wills-Toro and J. I. Zuluaga, Phys. Rev. D 77 (2008), 035008 [arXiv:0801.4044 [hep-ph]].
  * (9) Riazuddin and Fayyazuddin, Eur. Phys. J. C 56 (2008), 389-394 [arXiv:0803.4267 [hep-ph]].
  * (10) S. h. Nam, K. Y. Lee and Y. Y. Keum, Phys. Rev. D 82 (2010), 105027 [arXiv:0909.3770 [hep-ph]].
  * (11) G. Palacio, Int. J. Mod. Phys. A 31 (2016) no.25, 1650142 [arXiv:1608.08676 [hep-ph]].
  * (12) H. N. Long, L. T. Hue and D. V. Loi, Phys. Rev. D 94 (2016) no.1, 015007 [arXiv:1605.07835 [hep-ph]].
  * (13) M. Djouala, N. Mebarki and H. Aissaoui, Int. J. Mod. Phys. A 36 (2021) no.17, 17 [arXiv:1911.04887 [hep-ph]].
  * (14) D. Cogollo, Y. M. Oviedo-Torres and Y. S. Villamizar, Int. J. Mod. Phys. A 35 (2020) no.23, 2050126 [arXiv:2004.14792 [hep-ph]].
  * (15) L. T. Hue, K. H. Phan, T. P. Nguyen, H. N. Long and H. T. Hung, Eur. Phys. J. C 82 (2022) no.8, 722 [arXiv:2109.06089 [hep-ph]].
  * (16) A. Palcu, Mod. Phys. Lett. A 24 (2009), 2589-2600 [arXiv:0908.1636 [hep-ph]].
  * (17) A. Palcu, Int. J. Theor. Phys. 56 (2017) no.2, 403-414 [arXiv:1510.06717 [hep-ph]].
  * (18) 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]].
  * (19) K. Y. Lee and S. h. Nam, J. Phys. G 42 (2015) no.12, 125003 [arXiv:1412.1541 [hep-ph]].
  * (20) A. Crivellin, M. Hoferichter and P. Schmidt-Wellenburg, Phys. Rev. D 98 (2018) no.11, 113002 [arXiv:1807.11484 [hep-ph]].
  * (21) X. F. Han, T. Li, L. Wang and Y. Zhang, Phys. Rev. D 99 (2019) no.9, 095034 [arXiv:1812.02449 [hep-ph]].
  * (22) M. Endo and W. Yin, JHEP 08 (2019), 122 [arXiv:1906.08768 [hep-ph]].
  * (23) E. J. Chun and T. Mondal, JHEP 11 (2020), 077 [arXiv:2009.08314 [hep-ph]].
  * (24) L. Delle Rose, S. Khalil and S. Moretti, Phys. Lett. B 816 (2021), 136216 [arXiv:2012.06911 [hep-ph]].
  * (25) F. J. Botella, F. Cornet-Gomez and M. Nebot, Phys. Rev. D 102 (2020) no.3, 035023 [arXiv:2006.01934 [hep-ph]].
  * (26) S. P. Li, X. Q. Li, Y. Y. Li, Y. D. Yang and X. Zhang, JHEP 01 (2021), 034 [arXiv:2010.02799 [hep-ph]].
  * (27) I. Bigaran and R. R. Volkas, Phys. Rev. D 102 (2020) no.7, 075037 doi:10.1103/PhysRevD.102.075037 [arXiv:2002.12544 [hep-ph]].
  * (28) X. F. Han, T. Li, H. X. Wang, L. Wang and Y. Zhang, Phys. Rev. D 104 (2021) no.11, 115001 [arXiv:2104.03227 [hep-ph]].
  * (29) H. Bharadwaj, S. Dutta and A. Goyal, JHEP 11 (2021), 056 [arXiv:2109.02586 [hep-ph]].
  * (30) C. Arbeláez, R. Cepedello, R. M. Fonseca and M. Hirsch, Phys. Rev. D 102 (2020) no.7, 075005 [arXiv:2007.11007 [hep-ph]].
  * (31) K. F. Chen, C. W. Chiang and K. Yagyu, JHEP 09 (2020), 119 [arXiv:2006.07929 [hep-ph]].
  * (32) B. Dutta, S. Ghosh and T. Li, Phys. Rev. D 102 (2020) no.5, 055017 [arXiv:2006.01319 [hep-ph]].
  * (33) A. E. C. Hernández, S. F. King and H. Lee, Phys. Rev. D 103 (2021) no.11, 115024 [arXiv:2101.05819 [hep-ph]].
  * (34) A. E. C. Hernández, D. T. Huong and I. Schmidt, Eur. Phys. J. C 82 (2022) no.1, 63 [arXiv:2109.12118 [hep-ph]].
  * (35) S. Li, Z. Li, F. Wang and J. M. Yang, Nucl. Phys. B 983 (2022), 115927 [arXiv:2205.15153 [hep-ph]].
  * (36) F. J. Botella, F. Cornet-Gomez, C. Miró and M. Nebot, Eur. Phys. J. C 82 (2022), 915 [arXiv:2205.01115 [hep-ph]].
  * (37) L. Wang, J. M. Yang and Y. Zhang, Commun. Theor. Phys. 74 (2022) no.9, 097202 [arXiv:2203.07244 [hep-ph]].
  * (38) J. Kriewald, J. Orloff, E. Pinsard and A. M. Teixeira, Eur. Phys. J. C 82 (2022) no.9, 844 [arXiv:2204.13134 [hep-ph]].
  * (39) R. K. Barman, R. Dcruz and A. Thapa, JHEP 03 (2022), 183 [arXiv:2112.04523 [hep-ph]].
  * (40) R. Dermisek, Moscow Univ. Phys. Bull. 77 (2022) no.2, 102-107 [arXiv:2201.06179 [hep-ph]].
  * (41) T. A. Chowdhury, M. Ehsanuzzaman and S. Saad, JCAP 08 (2022), 076 [arXiv:2203.14983 [hep-ph]].
  * (42) C. H. Chen, C. W. Chiang and C. W. Su, “Top-quark FCNC decays, LFVs, lepton $g-2$, and $W$ mass anomaly with inert charged Higgses,” [arXiv:2301.07070 [hep-ph]].
  * (43) B. Abi et al. [Muon g-2], Phys. Rev. Lett. 126 (2021) no.14, 141801 [arXiv:2104.03281 [hep-ex]].
  * (44) G. W. Bennett et al. [Muon g-2], Phys. Rev. D 73 (2006), 072003 [arXiv:hep-ex/0602035 [hep-ex]].
  * (45) 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-166 (2020) [arXiv:2006.04822 [hep-ph]].
  * (46) 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]].
  * (47) T. Aoyama, M. Hayakawa, T. Kinoshita and M. Nio, Phys. Rev. Lett. 109, 111808 (2012) [arXiv:1205.5370 [hep-ph]].
  * (48) T. Aoyama, T. Kinoshita and M. Nio, Atoms 7, no.1, 28 (2019)
  * (49) 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]].
  * (50) C. Gnendiger, D. Stöckinger and H. Stöckinger-Kim, Phys. Rev. D 88, 053005 (2013) [arXiv:1306.5546 [hep-ph]].
  * (51) I. Danilkin and M. Vanderhaeghen, Phys. Rev. D 95, no.1, 014019 (2017) [arXiv:1611.04646 [hep-ph]].
  * (52) M. Davier, A. Hoecker, B. Malaescu and Z. Zhang, Eur. Phys. J. C 77, no.12, 827 (2017) [arXiv:1706.09436 [hep-ph]].
  * (53) A. Keshavarzi, D. Nomura and T. Teubner, Phys. Rev. D 97, no.11, 114025 (2018) [arXiv:1802.02995 [hep-ph]].
  * (54) G. Colangelo, M. Hoferichter and P. Stoffer, JHEP 02, 006 (2019) [arXiv:1810.00007 [hep-ph]].
  * (55) M. Hoferichter, B. L. Hoid and B. Kubis, JHEP 08, 137 (2019) [arXiv:1907.01556 [hep-ph]].
  * (56) 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]].
  * (57) A. Keshavarzi, D. Nomura and T. Teubner, Phys. Rev. D 101, no.1, 014029 (2020) [arXiv:1911.00367 [hep-ph]].
  * (58) A. Kurz, T. Liu, P. Marquard and M. Steinhauser, Phys. Lett. B 734, 144-147 (2014) [arXiv:1403.6400 [hep-ph]].
  * (59) K. Melnikov and A. Vainshtein, Phys. Rev. D 70, 113006 (2004) [arXiv:hep-ph/0312226 [hep-ph]].
  * (60) P. Masjuan and P. Sanchez-Puertas, Phys. Rev. D 95, no.5, 054026 (2017) [arXiv:1701.05829 [hep-ph]].
  * (61) G. Colangelo, M. Hoferichter, M. Procura and P. Stoffer, JHEP 04, 161 (2017) [arXiv:1702.07347 [hep-ph]].
  * (62) M. Hoferichter, B. L. Hoid, B. Kubis, S. Leupold and S. P. Schneider, JHEP 10, 141 (2018) [arXiv:1808.04823 [hep-ph]].
  * (63) A. Gérardin, H. B. Meyer and A. Nyffeler, Phys. Rev. D 100, no.3, 034520 (2019) [arXiv:1903.09471 [hep-lat]].
  * (64) J. Bijnens, N. Hermansson-Truedsson and A. Rodríguez-Sánchez, Phys. Lett. B 798, 134994 (2019) [arXiv:1908.03331 [hep-ph]].
  * (65) G. Colangelo, F. Hagelstein, M. Hoferichter, L. Laub and P. Stoffer, JHEP 03, 101 (2020) [arXiv:1910.13432 [hep-ph]].
  * (66) 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]].
  * (67) G. Colangelo, M. Hoferichter, A. Nyffeler, M. Passera and P. Stoffer, Phys. Lett. B 735, 90-91 (2014) [arXiv:1403.7512 [hep-ph]].
  * (68) V. Pauk and M. Vanderhaeghen, Eur. Phys. J. C 74, no.8, 3008 (2014) [arXiv:1401.0832 [hep-ph]].
  * (69) F. Jegerlehner, Springer Tracts Mod. Phys. 274, pp.1-693 (2017)
  * (70) M. Knecht, S. Narison, A. Rabemananjara and D. Rabetiarivony, Phys. Lett. B 787, 111-123 (2018) [arXiv:1808.03848 [hep-ph]].
  * (71) G. Eichmann, C. S. Fischer and R. Williams, Phys. Rev. D 101, no.5, 054015 (2020) [arXiv:1910.06795 [hep-ph]].
  * (72) P. Roig and P. Sanchez-Puertas, Phys. Rev. D 101, no.7, 074019 (2020) [arXiv:1910.02881 [hep-ph]].
  * (73) D. Hanneke, S. Fogwell and G. Gabrielse, Phys. Rev. Lett. 100 (2008), 120801 [arXiv:0801.1134 [physics.atom-ph]].
  * (74) R. H. Parker, C. Yu, W. Zhong, B. Estey and H. Müller, Science 360 (2018), 191 [arXiv:1812.04130 [physics.atom-ph]].
  * (75) L. Morel, Z. Yao, P. Cladé and S. Guellati-Khélifa, Nature 588 (2020) no.7836, 61-65
  * (76) T. Aoyama, M. Hayakawa, T. Kinoshita and M. Nio, Phys. Rev. Lett. 109 (2012), 111807 [arXiv:1205.5368 [hep-ph]].
  * (77) S. Laporta, Phys. Lett. B 772 (2017), 232-238 [arXiv:1704.06996 [hep-ph]].
  * (78) T. Aoyama, T. Kinoshita and M. Nio, Phys. Rev. D 97 (2018) no.3, 036001 [arXiv:1712.06060 [hep-ph]].
  * (79) H. Terazawa, Nonlin. Phenom. Complex Syst. 21 (2018) no.3, 268-272
  * (80) S. Volkov, Phys. Rev. D 100 (2019) no.9, 096004 [arXiv:1909.08015 [hep-ph]].
  * (81) A. Gérardin, Eur. Phys. J. A 57 (2021) no.4, 116 [arXiv:2012.03931 [hep-lat]].
  * (82) B. Aubert et al. [BaBar], Phys. Rev. Lett. 104 (2010), 021802 [arXiv:0908.2381 [hep-ex]].
  * (83) A. M. Baldini et al. [MEG], Eur. Phys. J. C 76 (2016) no.8, 434 [arXiv:1605.05081 [hep-ex]].
  * (84) A. Abdesselam et al. [Belle], JHEP 10 (2021), 19 [arXiv:2103.12994 [hep-ex]].
  * (85) A. M. Baldini et al. [MEG II], Eur. Phys. J. C 78 (2018) no.5, 380 [arXiv:1801.04688 [physics.ins-det]].
  * (86) E. Kou et al. [Belle-II], PTEP 2019 (2019) no.12, 123C01 [erratum: PTEP 2020 (2020) no.2, 029201] [arXiv:1808.10567 [hep-ex]].
  * (87) D. Chang and H. N. Long, Phys. Rev. D 73 (2006), 053006 [arXiv:hep-ph/0603098 [hep-ph]].
  * (88) L. T. Hue, A. E. Cárcamo Hernández, H. N. Long and T. T. Hong, Nucl. Phys. B 984 (2022), 115962 [arXiv:2110.01356 [hep-ph]].
  * (89) F. Jegerlehner and A. Nyffeler, Phys. Rept. 477, 1-110 (2009) [arXiv:0902.3360 [hep-ph]].
  * (90) P. A. Zyla et al. [Particle Data Group], PTEP 2020, no.8, 083C01 (2020)
  * (91) J. A. Casas and A. Ibarra, Nucl. Phys. B 618, 171-204 (2001) [arXiv:hep-ph/0103065 [hep-ph]].
  * (92) B. Pontecorvo, Sov. Phys. JETP 6, 429 (1957)
  * (93) Z. Maki, M. Nakagawa and S. Sakata, Prog. Theor. Phys. 28, 870-880 (1962)
  * (94) E. Arganda, M. J. Herrero, X. Marcano and C. Weiland, Phys. Rev. D 91, no.1, 015001 (2015) [arXiv:1405.4300 [hep-ph]].
  * (95) N. H. Thao, L. T. Hue, H. T. Hung and N. T. Xuan, Nucl. Phys. B 921, 159-180 (2017) [arXiv:1703.00896 [hep-ph]].
  * (96) L. Lavoura, Eur. Phys. J. C 29, 191-195 (2003) [arXiv:hep-ph/0302221 [hep-ph]].
  * (97) L. T. Hue, L. D. Ninh, T. T. Thuc and N. T. T. Dat, Eur. Phys. J. C 78, no.2, 128 (2018) [arXiv:1708.09723 [hep-ph]].
  * (98) L. T. Hue, H. T. Hung, N. T. Tham, H. N. Long and T. P. Nguyen, Phys. Rev. D 104 (2021) no.3, 033007 [arXiv:2104.01840 [hep-ph]].
  * (99) T. T. Hong, N. H. T. Nha, T. P. Nguyen, L. T. T. Phuong and L. T. Hue, PTEP 2022 (2022) no.9, 093B05 [arXiv:2206.08028 [hep-ph]].
  * (100) N. Aghanim et al. [Planck], Astron. Astrophys. 641, A6 (2020) [erratum: Astron. Astrophys. 652, C4 (2021)] [arXiv:1807.06209 [astro-ph.CO]].
  * (101) E. Fernandez-Martinez, J. Hernandez-Garcia and J. Lopez-Pavon, JHEP 08, 033 (2016) [arXiv:1605.08774 [hep-ph]].
  * (102) 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]].
  * (103) A. M. Coutinho, A. Crivellin and C. A. Manzari, Phys. Rev. Lett. 125, no.7, 071802 (2020) [arXiv:1912.08823 [hep-ph]].
  * (104) C. Biggio, E. Fernandez-Martinez, M. Filaci, J. Hernandez-Garcia and J. Lopez-Pavon, JHEP 05, 022 (2020) [arXiv:1911.11790 [hep-ph]].
  * (105) T. Mondal and H. Okada, Nucl. Phys. B 976 (2022), 115716 [arXiv:2103.13149 [hep-ph]].
  * (106) P. Escribano, J. Terol-Calvo and A. Vicente, Phys. Rev. D 103, no.11, 115018 (2021) [arXiv:2104.03705 [hep-ph]].