$(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

We will show that the 3-4-1 model with heavy right-handed neutrinos can explain the recent experimental data of $(g-2)_{e, \mu}$ anomalies 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.


Introduction
The 3-4-1 model with right-handed neutrinos (341RHN) was discussed in Ref. [1,2] as a natural extension that new right-handed neutrinos are assigned into the same lefthanded lepton quadruplets. For a complete study of the highest possible gauge symmetry in the electroweak sector [3], various 3-4-1 extensions were introduced with different electric charges of new exotic leptons [4][5][6][7][8][9][10][11][12]. 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 [13]. A solution applied to the 3-3-1 models [14], 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 ea ≡ (g − 2) ea /2, neutrino oscillation data, and the recent bounds of the lepton flavor violating decays of charged leptons (cLFV) e b → e a γ. The explanations of neutrino oscillation data were mentioned previously in various 3-4-1 models, including the ISS mechanism [10,15,16], 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 [16,17], 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 [18]. 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 ea enough to be consistent with experiments [19]. 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 → e a γ will be discussed precisely in the 3-4-1 framework. Besides many available models beyond the standard model (BSM)  under various experimental data including the (g − 2) e,µ anomalies, our work will confirm the reality of the 3-4-1 models.
The latest experimental measurement for muon anomaly a µ has been reported from the combination of the two experimental results at Fermilab [43] and Brookhaven National Laboratory (BNL) E82 [44]: a exp µ = 116592061(41) × 10 −11 . It leads to a standard deviation of 4.2 σ (standard deviation) from the Standard Model (SM) prediction, namely ∆a NP µ ≡ a exp µ − a SM µ = (2.51 ± 0.59) × 10 −9 , where a SM µ = 116591810(43) × 10 −11 is the SM prediction [45] combined from various different contributions based on the dispersion approach . Another larger SM value calculated using the lattice-QCD technique implies a smaller value of ∆a NP µ than that given in Eq. (1) [107]. In this work we will accept the experimental constraint from Eq. (1) consisting of both results. On the other hand, the recent experimental a e data was reported from different groups [73][74][75]108], leading to the two inconsistent deviations between experiments and the SM prediction [76][77][78][79][80][81]. In this work, we accept the following value 1 : ∆a NP e ≡ a exp e − a SM e = (3.4 ± 1.6) × 10 −13 , where we have used the latest experimental result of a exp e reported in Ref. [108], consistent with the 2008 result of the same group [73], and the a SM e value reported in Ref. [75], derived indirectly from the measurement of the fine-structure constant α using the Rb atoms, corresponding to the 2.1σ deviation. There is another a SM e value corresponding to the measurement of the fine-structure constant of Cs-133 atoms [74], which is inconsistent with the earlier, namely ∆a NP e = (−10.2 ± 2.6) × 10 −13 , with the 3.9σ deviation. Both results can be explained in this work, see details later in the numerical investigation The ISS mechanism may result in large values of not only (g − 2) e,µ but also Br(e b → e a γ), which are constrained by recent experiments [82][83][84] Br(τ → µγ) < 4.4 × 10 −8 , Br(µ → eγ) < 4.2 × 10 −13 .
Our paper is arranged as follows. Section 2 presents all ingredients of a 3-4-1 model to calculate the (g − 2) ea data and cLFV decays. Section 3 introduces the 341ISS model to accommodate the recent (g − 2) ea data. In this model, the Yukawa Lagrangian and Higgs potential must respect a new global U (1) 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) ea anomalies. Section 4 will present detailed numerical results to determine the allowed regions of the parameter space that explain both experimental results of two (g − 2) e,µ anomalies and cLFV decays. Section 5 summarizes important results. 2 The model with Dirac active neutrinos

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 [5,10]. This model is constructed based on the gauge symmetry SU (3) c × SU (4) L × U (1) X , implying 16 electroweak (EW) gauge bosons. In addition, there are 4 neutral gauge bosons corresponding to the 4 diagonal generators of the EW group. Normally, the EW group is assumed to break to the final electric group through the following pattern: i.e, the 3-4-1 model can be considered as the extended version of 3-3-1 models [112,113]. The electric operator is defined as follows: where T 3,8,15 are the diagonal generators of the SU (4) group and X is the U (1) charge, see precise forms corresponding to the quadruplet in Ref. [114], for example.
The Higgs multiplets and non-zero vacuum expectation values (VEV) of neutral components needed for generating all fermion masses are: The lepton masses are generated from the following Yukawa interactions −L y = Y N ab L a χN ′ bR + Y E ab L a ϕE ′ bR + Y e ab L a ρe ′ bR + Y ν ab L a ην ′ bR + H.c..
The model consists of quark multiplets that must be arranged to cancel the gauge anomalies, see for example a discussion in Ref. [11]. It can be seen that the quark masses can be constructed to satisfy the recent experimental data.
The zero VEV values of some neutral Higgs components can be explained by considering a global symmetry called the general lepton number L defined as follows: where L is the normal lepton number. This formula is an extension of that introduced for a 3-3-1 model [87].
leptons L(ν ′ aL,R ) = L(e ′ aL,R ) = 1, and all SM quarks have L = 0. The remaining non-zero lepton numbers L are listed in Table 2, because of the requirement that the total lepton Table 2 Nonzero lepton number L of all fields in the 341RH.
number L is always conservative.
The mass terms of all leptons are: The active Dirac neutrino masses and mixing are constructed from the mass matrix M ν .
But, these tiny masses do not affect significantly the one-loop contributions to a ea . 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 ′ a is the mass basis e aL,R ≡ e ′ aL,R , namely Three other base f ′ L,R ≡ (f ′ 1 , f ′ 2 , f ′ 3 ) T L,R with f = ν, E, N are transformed into the corresponding mass base f ′ L,R through the following relations: Although the quark sector is irrelevant to our work, we review here the main property relating to the recent experimental constraints of K 0 −K 0 oscillation, similar to the 3-3-1 models [112,113]. The quark sector consists of two families of left-handed anti-quadruplets, one family of left-handed quadruplet, and singlets of right-handed quarks, namely The relevant Yukawa parts of quarks respecting all gauge and L symmetries are We can see that the Yukawa parts generating SM quark masses are the same as those shown in Refs. [112,113] for 3-3-1 models with right-handed neutrinos, where both ρ and χ inherit the same VEV properties. In addition, all exotic quarks do not mix which the SM quarks, consistent with the property shown in Table 2 that all SM quarks have zero L, in contrast with L = ±2 for all exotic quarks. Therefore, the couplings of the SM quarks with Higgs bosons in the models given in Refs. [112,113] are the same as those presented in the model under consideration. Accordingly, the mass split value ∆m K originated from the K 0 -K 0 oscillation predicted by the 3-4-1 model under considerations will be the sum of the new heavy Higgs and neutral gauge contributions at the tree level, in which every contribution is proportional to the inverse of the squared mass of the respective Higgs or gauge boson.
As we will see below, the 3-4-1 model predicts two new heavy neutral gauge bosons Z 3 and Z 4 with masses m 2 Z 3 ∝ w 2 and m 2 Z 4 ∝ V 2 in the limit V ≫ w. Because w plays the role of the SU (3) L scale, Ref. [112] confirms that if m Z 4 ≫ m Z 3 > 4 TeV, contributions from these two gauge bosons to ∆m K will satisfy the experimental constraint. The case of the heavy Higgs contributions are the same. Because this choice of parameters does not affect the main results in our work, this experimental constraint will be ignored in the remaining part of this work.

Gauge boson masses and mixing
Gauge boson masses arise from the covariant kinetic term of Higgs multiplets, namely where H = χ, ϕ, η, ρ. The covariant derivative is defined as where g, g X and W aµ , B ′′ µ are gauge couplings and fields of the gauge groups SU (4) L and U (1) X , respectively. The two parts P N C µ and P CC µ relate to the neutral and non-hermitian currents [11]. For quadruplet, T 16 = 1 2 √ 2 diag(1, 1, 1, 1) and where t ≡ g X /g and The upper subscripts label the electric charges of gauge bosons. The relation between the original basis (W 3 , W 8 , W 15 , B ′′ ) and the mass basis (A, Z, Z 3 , Z 4 ) of all real neutral boson was determined previously [11]. The masses of non-Hermitian (charged) gauge bosons are given by By spontaneous symmetry breaking (SSB), the following relation should be in order: V ≫ ω ≫ v 1 , v 2 . A consequence from Eq. (18) is that W ± must be identified with the singly charged SM gauge boson, namely Then neutral gauge boson masses are The above formula implies that η and ρ 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 × U (1) Y . Then we define the mixing angle β as follows where t β ≥ 0.4 as the perturbative limit of the top quark Yukawa coupling |Y u The mixing parameters of the neutral gauge bosons were summarized in appendix B, see the details in Ref. [11]. Because the current study [112] shows that the breaking scale of SU (3) L and SU (L) L is the order of TeV, the mixing parameters between Z, Z 3 , and Z 4 are much small, therefore we ignore in this work. But the loop corrections to the Zµ + µ − may be significant, especially in the model inheriting the "chiral enhancement" enough to accommodate the (g − 2) ea anomaly data. These one-loop corrections were computed in appendix B used to constrain the experimental data in the numerical investigation.

Higgs potential and Higgs spectrum
The most general Higgs potential including the appearance of the singly charged Higgs boson σ ± ∼ (1, 1, ±1) is: Here we consider the Higgs potential respecting the generalized lepton number L and the new singly charged Higgs boson has L(σ ± ) = 0. We also ignore the triple Higgs self couplings (f ′ σ χ † ϕ + 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. [11] without σ ± . The Higgs potential consisting of σ ± was discussed in 3-3-1 models [14,88], 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 χ ± 2 and ρ ± 4 results in two massless Goldstone boson G ± 24 and two singly charged Higgs bosons h ± 3 : and We consider here the mixing of ϕ 0 * 2 and ρ 0 3 , which leads to a non-hermitian Goldstone boson G 0 23 ̸ = G 0 * 23 and a physical state H 0 and Three singly charged Higgs bosons (ρ ± 1 , η ± 2 , σ ± ) are changed into the two physical states h ± 1,2 and the Goldstone bosons G ± W of W ± as follows: The relations in Eq. (28) were given in Refs. [14,88], which consist of the same part V (σ) given in Eq. (23) in the Higgs potential. The mixing parameter α and Higgs boson masses will be investigated as free parameters, while three dependent parameters are: As usual, this model must contain a Standard model-like (SM-like) Higgs boson confirmed by LHC. A detailed calculation to identify this Higgs boson in the model under consideration can be found in Ref. [11]. A simple estimation is presented in appendix A. In general, the SMlike Higgs boson gets main contributions from η and ρ, similarly to the case of the property of the well-known two Higgs doublet model (2HDM). In the model under consideration, the flavor and the physical base of charged leptons are the same, as assumed in Eq. (11).
As a result, there is no mixing between SM and heavy charged leptons, implying that the leading corrections to the SM coupling hµ + µ − is the combination of the tree-level of the Higgs mixing and the one-loop level [106,110,111,116]. As commented in Ref. [106], the recent experiments give no hint for h → e + e − , while only the evidence of h → µ + µ − was reported [117,118]. The loop corrections to the coupling hµ + µ − is smaller than the recent experimental sensitivities, but may be detected in the future [106,111,116]. Therefore, the constraint from the modification of the coupling hµ + µ − will not affect the allowed regions of the parameter space we will discuss in this work.

Analytic formulas for one-loop contributions to a ea with active Dirac neutrinos
In this section, we consider the simple case that ν 1,2,3 are Dirac ones, and no mixing between σ ± and other singly charged Higgs bosons, i.e. s α = 0, c α = 1, then only h ± 2 couples with active neutrinos. The conclusion with s α ̸ = 0 is unchanged because the tiny neutrino masses result in suppressed one-loop contributions to (g − 2) ea anomalies.
The relevant Lagrangian giving one-loop contributions to a ea is: To avoid large cLFV decays e b → e a γ, 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 Ea = m E , m Na = m N for all a = 1, 2, 3.
For active neutrinos having tiny masses, the respective one-loop contributions to ∆a ea are suppressed. The non-zero form factors relevant with one-loop contributions to a ea are [19] The total contribution to ∆a ea from all Higgs bosons are: where X = h + 1 , h + 3 , H 0 1 . The functions relating to one-loop contributions of Higgs bosons are shown in Fig. 1.
for all x > 0, the one-loop contributions from singly charged Higgs bosons h ± 1,3 are in opposite signs with ∆a NP µ , hence they should be small. On the other hand, the one-loop contribution from the neutral Higgs boson H 0 The couplings of leptons and non-hermitian gauge bosons are (34) where the last line collects only terms giving one-loop contributions to cLFV amplitudes and (g − 2) anomalies, resulting in the following formulas [19], 0.01 0.10 1 10 100 1000 Properties of the scalar functions relating to one-loop contributions of heavy Higgs bosons to (g − 2) ea anomalies in the 3-4-1 model with active Dirac neutrinos. wheref The total one contribution from new heavy gauge bosons to the a ea is The deviation of a µ between predictions by the two models 3-4-1 and SM are where a SM µ (W ) = 3.887 × 10 −9 [89], and a SM ea (W ) = a SM µ (W ) × (m 2 ea /m 2 µ ) are the SM's prediction for the one-loop contribution from W boson. In this work, ∆a 341 ea = ∆a ea will be considered as new physics (NP) predicted by the 3-4-1 models, which must satisfy the experimental data given in Eqs. ∆a µ . The ∆a µ depends strongly on w, which lower bound is constrained strictly from the masses, which are given in Eq. (20), of heavy neutral gauge bosons Z 3,4 , namely m Z 3,4 > 3 TeV from LHC searches [18]. As a result, large values of w ≥ 2 TeV give ∆a µ ≤ 7 × 10 −13 ≪ 192 × 10 −13 ≤ a NP µ . In conclusion, all the one-loop contributions mentioned above are much smaller than a NP µ .

The 341ISS with ISS neutrinos
Now we consider an extension of the above 341RHN model that may explain successfully the (g − 2) ea data. This version consists of six right-handed neutrinos ν aR , X aR ∼ (1, 0), a = 1, 2, 3 generating active neutrino masses through the ISS mechanism, and a singly charged Higgs boson σ + ∼ (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 L(ν aR ) = 1, L(X aR ) = −1, and L(σ + ) = 0. Now tree-level neutrino masses and mixing angles arise from the ISS mechanism. Requiring that L is only softly broken, the additional Yukawa part is where Y ν , M R , µ X , and Y σ are 3 × 3 matrices. The first term in Eq. (41) 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 ν L = (ν ′ 1 , ν ′ 2 , ν ′ 3 ) T L and ν R = (ν 1 , ν 2 , ν 3 ) T R , X R = (X 1 , X 2 , X 3 ) T R , the neutrino mass terms derived from (41) are written in the following ISS form [90]: where O 3×3 is a zero matrix and m D = Y ν × v 2 / √ 2. The analytic form of the Dirac mass matrix was chosen generally following Ref. [91]. The total unitary mixing matrix U ν is defined as follows 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 matrixm ν and six other heavy neutrinos with mass matrixM N . The relation between the flavor and mass eigenstates are where n L ≡ (n 1L , n 2L , ..., n 9L ) T and n R = (n L ) c . The neutrino mixing matrix is parameterized in the following form: where U PMNS is the 3 × 3 Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix [92,93], V is a 6 × 6 unitary matrix, and R is a 3 × 6 matrix satisfying |R aI | ≪ 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 convenient to choose a simple form of µ X = µ 0 I 3 and [94,95]. The formulas of m D and mixing parameters are [91] where max[|(x ν ) aa |] ≪ 1 for all a = 1, 2, 3.
≃ 0, the mixing matrix and Majorana mass term arê  We note that except the contributions from ISS neutrino couplings with singly charged Higgs bosons given in Eq. (28), 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 Up to the order O(R 2 ) of the neutrino mixing matrix given in Eq. (45), the non-zero one-loop contributions relating to h ± 1,2 are The one-loop contributions from the h ± 1,2 exchanges to a ea are: where and a ea (h ± ) is the total contribution from these two Higgs bosons. We note that all Yukawa couplings considered in this work are assumed to be real. As a result, Eq. (52) shows that only the Dirac phase results in the tiny values of Im[c aa ], implying very suppressed values of electric dipole moments d ea ≡ −2 Im[c aa ] [19], consistent with experimental constraints [123,124]. This conclusion is also confirmed by our numerical investigation, namely |d e,µ | < 10 −44 [e cm], therefore we will not discuss it further.
For qualitative estimation, the main contribution to a ea (h ± ) is the chirally-enhanced part [19]. From Eq. (52), this part is determined as follows 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 ̸ = a. To avoid this problem, we start from the diagonal form of c (ab)R,0 as follows Correspondingly, the formula of a ea,0 is proportional to a diagonal matrix Y d satisfying: Hence, the diagonal entries will give main contributions to a ea , namely Then a ea,0 may be large, provided that t β should not too large, t 1 ̸ = t 2 , and the following quantities are large enough: x 0 , |s 2α = 2s α c α |, 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, ∆a ea ≡ a ea (h ± ) given in Eq. (52), instead of the total formula given in Eq. (40). In contrast, the one-loop contribution from W must be included in the formulas of Br(e b → e a γ). 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 [98,99].

Numerical discussion
We will use the best-fit values of the neutrino osculation data [90] corresponding to the normal order (NO) scheme with m n 1 < m n 2 < m n 3 , namely This choice of active neutrino masses also satisfies the constraint from Plank2018 [100], The non-unitary of the active neutrino mixing matrix I 3 − 1 2 RR † U PMNS is constrained by other phenomenology such as electroweak precision [101][102][103], leading to a very strict constraint of η ≡ 1 2 RR † ∝x ν in the ISS framework [104][105][106]. The reasonable constraint of the non-unitary part of the neutrino mixing matrix is as follows The other well-known numerical parameters are [90] g = 0.652, G F = 1.1664 × 10 −5 GeV, s 2 W = 0.231, m W = 80.385 GeV, m e = 5 × 10 −4 GeV, m µ = 0.105 GeV, m τ = 1.776 GeV.
For the free parameters of the 341ISS model, the numerical scanning ranges are In addition, we will fix m n 1 = 0.01 eV, and check the perturbative limit of all Yukawa couplings Y σ and Y ν , namely |Y σ ab |, |Y ν ab | ≤ 3.5 must be satisfied. As we discussed above, the diagonal form of Y d will allow large (g − 2) ea anomalies and small Br(e b → e a γ) satisfying the recent experimental constraints. In the numerical investigation, we will consider the general case that Y d ab ̸ = 0. The allowed regions of parameters we imply below must guarantee both the experimental data of 1σ ranges of (g − 2) e,µ , the cLFV constraints of Br(e b → e a γ) , and the constraint of the decay Z → µ + µ − using the analytic formulas and experimental data given in appendix B. 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 ̸ = 0, the allowed regions satisfying two (g − 2) e,µ data are still strictly constrained by the cLFV decay rate Br(µ → eγ) < 4.2 × 10 −13 , namely the constraints |Y d 12 |, |Y d 21 | < 10 −4 must be guaranteed. Therefore, we will fix |Y d 12 | = |Y d 21 | = 0 as the default values in our numerical investigations. Second, we give comments on the three following cases:  Fig. 3. In this allowed region, Δa e (h ± )×10 13 Δa μ (h ± )×10 9  54) is the main one-loop contribution to (g − 2) ea 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. (52).  Fig. 4. We conclude that Br(τ → eγ) depends Br(τ→eγ)

-12
10 -11 10 -10 10 -9 10 -8 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) ea anomalies, while the non-zero entries Y d 13,31 and Y d 23,32 affect strongly in the decay rates of Br(τ → eγ) and Br(τ → µγ), respectively. The Br(µ → eγ) depends strongly on Δa μ (h ± )×10 9 x 0 Δa Δa e (h ± )×10 13 We can see that many allowed ranges are stricter than the scanning ones given in Eq. (62). 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 µ (h ± ), see the bottom right panel. As a result, small (g − 2) µ may predict large Br(τ → µγ, µγ) and vice versa. In addition to explain both (g − 2) ea data 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 [24,25].  Δa e (h ± )×10 13 Fig. 6 The correlations between Y d ij vs ∆a e (h ± ) with Y 21 = Y 12 = 0.
confirms the relation given in Eq. (57). The a µ (h ± ) ∝ |Y d 22 | is not very clear, because many points are excluded by the perturbative limit. The right panel shows that a e (h ± ) does not give any prediction on Y d ij like the case of a µ (h ± ). The correlations of Y d ij vs Br(µ → eγ) are shown in Fig. 7. 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) ea data. In the right panel of Fig. 7, large Br(µ → eγ) prefers large non-diagonal entries Y d ij but the Br(μ→eγ)  1σ experimental ranges still predict large Br(e b → e a γ) 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. Note also that future experimental constraints on other cLFV decays such as µ → 3e < 10 −16 [121] and µ − e conversion in nuclei [122] will be more strict than those considered in this work, but in general theoretical calculations on some other BSM in the presence of "chiral enhancement" show that they do not affect strongly on the allowed regions we discussed above, see for example the left-right model [115]. From the theoretical side, this can be explained by the reason that the relevant one-loop corrections originated from one-loop four-point diagrams are proportional to the products of four vertex factors, therefore their numerical values are more flexible than the theoretical constraints from the cLFV decays e b → e a γ. A more detailed investigation to predict the allowed regions corresponding to the future sensitivities of all interesting cLFV decays will be done in the future.

Conclusion
We have discussed a solution to explain the recent experimental data of the (g − 2) ea anomalies in the 341ISS framework. We have constructed the Yukawa Lagrangian of leptons and Higgs potential obeying the generalized lepton number L that keeps necessary terms generating the ISS mechanism and large chirally-enhanced one-loop contributions to (g − 2) ea anomalies. Although the ISS mechanism may result in large cLFV decays e b → e a γ, 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 τ → µγ, eγ, or µ → eγ once they are detected by incoming experiments.
As a result, the SM-like Higgs boson gets dominant contributions from the neutral Higgs basis (Re[ρ 0 2 ], Re[η 0 1 ]), corresponding to the following squared mass matrix: Similarly the 2HDM framework, this CP-even Higgs boson can be identified with the SM-like Higgs boson found from LHC. Denoting the mixing parameter α of these two Higgs bosons, the difference between the tree level couplings he a e a predicted by the model under consideration and the SM is a factor s α /c β , exactly the same as the 2HDM, see Ref. [125] for example. Therefore, the combination of this factor and the loop corrections can accommodate the experimental data.
B One-loop corrections to Zµ + µ − and hµ + µ − vertices The relations between the flavor and physical base of the neutral gauge bosons is [11]  The matching relations from the breaking step SU (4) L × U (1) X → SU (2) L × U (1) Y is W 3 = W 3 , B µ = c 32 W 8µ + c 43 s 32 W 15µ + s 43 s 32 B ′′ µ , and where B µ is the gauge boson of the gauge U (1) Y in the SM. The relation (B3) is consistent with the definition of the charge operator (4), the same as that given in the SM. This means that we can consider η and σ, ν aR and X aR as new scalars and neutral fermions giving oneloop corrections to the Zµ + µ − couplings as discussed in Ref. [109]. Because all new neutral fermions are singlets, vertices Zν aRνaR and ZX aRXaR do not appear. Consequently, only diagrams 2 of Fig. 1 in Ref. [109] is irrelevant to our model.