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

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

I.
Many recent versions of the 3-3-1 models were indicated that they are difficult to explain simultaneously all of these experimental constraints [28][29][30][31][32][33] with the very large TeV values of the SU (3) L symmetry scale. Namely, the discussion on Ref. [29] needs the cLFV constraints from experimental data to rule out large ∆a µ . The remaining models rule out large ∆a µ for large SU ( ∆a µ is a popular way to explain successful the experimental data of a µ [31,32], but there seems irrelevant with neutrino oscillation data. Some recent extensions of 3-3-1 models with discrete symmetries [34,35] need a large number of new leptons and Higgs bosons for the explanation of large ∆a µ,e consistent with experiments. On the other hand, a recent note indicated that a version of the 3-3-1 model with right-handed neutrino (331RN) with heavy neutral fermions assigned as SU (3) L gauge singlets (called the 331ISS model for short) can predict large one-loop contributions from singly charged Higgs bosons and inverse seesaw (ISS) neutrinos enough to explain the recent (g − 2) µ data [36]. More interesting, the model contains two singly charged Higgs bosons, which may result in a special possibility that two one-loop contributions to ∆a µ are large and constructive, while those relate with cLFV decay amplitudes are strongly destructive. In this work, we will pay attention to this possibility, namely we will try to answer a question whether there exist any allowed regions of the parameter space that the destructive properties of the Higgs contributions are enough to satisfy the cLFV experimental constraints given in Eq. (2), and explain successfully the recent data given in Eq. (1). We will use the 3-3-1 model with the general Higgs potential given in Ref. [37,38]. The 3-3-1 models explaining active neutrino data based on the ISS mechanism has been discussed widely previously [39][40][41][42], but the interesting regions of the parameter space allowing large ∆a µ data and consistent with recent cLFV experimental constraints were not shown. In addition, the Br(τ → µγ, eγ) were predicted to be smaller than Br(µ → eγ), which is very suppressed with the recent and upcoming experimental sensitivities of the order O(10 −9 ) [43,44]. Many other models beyond the SM with the ISS mechanism can explain consistently the experimental data of ∆a µ and cLFV constraints [45][46][47][48]. Here we analyze predictions of the 3-3-1 model with right-handed neutrinos for the above observables.
Our work is arranged as follows. We will review the 331ISS model in Sec. II, summarize the gauge, Higgs bosons and the lepton sectors. In Sec. III, we introduce the analytic formulas to calculate the muon magnetic dipole moment and the cLFV branching ratios. In Sec. IV, we discuss on the effect of a new singly charged Higgs boson that can give one-loop contributions to ∆a µ and cLFV amplitudes enough to explain successful all the experimental data under consideration. In Sec. V, illustrations for numerical results are given to indicate the existence of the allowed regions satisfying the experimental data mentioned in this work.
The conclusion is presented in the last Sec. VI, where important results will be summarized.

A. Gauge bosons and fermions
The particle content of the 331ISS model was introduced in Refs. [41,49] where active neutrino masses and oscillations are originated from the ISS mechanism. The quark sector and SU (3) C representations are irrelevant in this work, and hence they are omitted here.
We refer Ref. [41] for a quark discussion. The electric charge operator corresponding to 3 ) and a right-handed charged lepton e aR ∼ (1, −1) with a = 1, 2, 3. Each left-handed neutrino N aL = (N aR ) c is equivalent with a new right-handed neutrinos defined in previous 331RN models [50]. The only difference between the two models 331RN and 331ISS is that, the 331ISS model contains three more right-handed neutrinos transforming as gauge singlets, X aR ∼ (1, 0), a = 1, 2, 3. They couple with the SU (3) L Higgs triplets to generate the neutrino mass term relating with the ISS mechanism. The three Higgs triplets ρ = ) have the following necessary vacuum expectation values for generating all tree-level quark masses The gauge bosons get masses through the covariant kinetic term of the Higgs triplets, where the covariant derivative for the electroweak symmetry is D µ = ∂ µ −igW a µ T a −ig X T 9 XX µ , a = 1, 2, .., 8. Note that T 9 ≡ I 3 √ 6 and 1 √ 6 for (anti)triplets and singlets [51]. Matching with the SM gives e = g s W and g X g = 3 , where e and s W are respective the electric charge and sine of the Weinberg angle, s 2 W 0.231. The relation g X g is the same for both choices of triplet or antitriplets representations of the left-handed leptons [52,53]. The derivation of this relation is summarized as follows. The 3-3-1 models have two spontaneous breaking steps: The first breaking step with w = 0 generates masses for heavy particles predicted by the SU (3) L symmetry. The neutral gauge bosons will change into the basis containing the SM ones W 3 µ and B µ : . Diagonalizing the squared mass matrix of these neutral gauge bosons will get a massive eigenstate Z with m 2 Z ∼ w 2 and two SM massless states W 3 µ and B µ . The relations between the two bases before and after the first breaking step are Inserting these relations to the covariant derivation of the 3 − 3 − 1 gauge group and keeping the part used to identify with the SM one, we have which results in the consequences that g and gt √ 6+β 2 t 2 = gt W are the gauge couplings of the SM, and the U (1) Y charge of the SM is Y /2 = βT 8 + IX.
Like the 331RN model, the 331ISS model includes two pairs of singly charged gauge bosons with the following physical states W ± and Y ± and masses The bosons W ± are identified with the SM ones, leading to the consequence that The general Higgs potential relating with the 331RN model will be applied in our work with v 1 = v 2 . We will use the following parameters for this general case.
The parameter t β plays a similar role known in the well-known models with two Higgs doublet and the minimal supersymmetric Standard Model. This is different from Ref. [49], where v 1 = v 2 was assumed so that the Higgs potential given in Ref. [54] was used to find the exact physical state of the SM-like Higgs boson. This simple condition was also used in previous discussions in 3-3-1 models addressed with anomalous magnetic dipole moments [31,32]. As we will show below, large t β = 1 is one of the key condition for predicting large (g − 2) µ consistent with experiments. The reason is that the physical states of the charged Higgs bosons are determined analytically from this Higgs potential, and only these Higgs bosons contribute significantly to one-loop corrections to the (g − 2) µ .
The Yukawa Lagrangian for generating lepton masses is: Here we assumed that the model under consideration respects a new lepton number symmetry L discussed in Ref. [38] so that the term ψ aL ηX bR is not allowed in the above Yukawa Lagrangian, while the soft-breaking term (µ X ) * ba (X aR ) c X bR is allowed with small (µ X ) ba . The new lepton number L called by generalized lepton number [55] is defined as where L is the normal lepton number. The specific assignment of L is L(ρ) = −1/3, L(η) = −2/3, L(χ) = 4/3, L(ψ aL ) = 1/3, which guarantees the consistence for the well-known definition of L, namely L( ) = 1 for = e aL,R , ν aL , L( ) = −1 for = N aL , X aR , and L(q) = 0 for all SM quarks [38].
The first term in Lagrangian (6) generates charged lepton masses m ea ≡ h e ab v 1 √ 2 δ ab , i.e, the mass matrix of the charged leptons is assumed to be diagonal, hence the flavor states of the charged leptons are also the physical ones. In the basis ν L = (ν L , N L , (X R ) c ) T and (ν L ) c = ((ν L ) c , (N L ) c , X R ) T of the neutral leptons, Lagrangian (6) gives a neutrino mass term corresponding to a block form of the mass matrix [49], namely where The mass matrix M R does not appear in the 331RN. The Dirac neutrino mass matrix m D must be antisymmetric. The matrix µ X defined in Eq. (6) is symmetric and it can be diagonalized by a transformation U X : The matrix U X will be absorbed by redefinition the states X a , therefore µ X will be set as the diagonal matrix given in the right hand side of Eq. (8).
The mass matrix M ν is diagonalized by a 9 × 9 unitary matrix U ν , where m n i (i = 1, 2, ..., 9) are masses of the nine physical neutrino states n iL . They consist of three active neutrinos n aL (a = 1, 2, 3) corresponding to the mass subma- , and the six extra neutrinos n IL (I = 4, 5, .., 9) with M N = diag(m n 4 , m n 5 , ..., m n 9 ). The ISS mechanism leads to the following approximation solution of U ν , where The relations between the flavor and mass eigenstates are where n L ≡ (n 1L , n 2L , ... where In this paper, we will work on the normal ordered scheme (NO) of the active neutrino masses, which allows δ = π using in this work. The respective best fit and the confidence level of 3σ of the neutrino oscillation experimental data is given as [4] Additionally, it is easily to derive that The detailed calculation shown in Ref. [49], using the ISS relations, yields where z = √ 2v h ν 23 is assumed to be positive and real, We note that the lightest active neutrino mass is zero at the tree level, but can be nonzero when loop-corrections are included [38]. Also, the quantum effects can be considered for the charged lepton masses, so that the regions predicting large ∆a µ may be larger [58,59] than the ones discussed in this work. The perturbative limit requires that h ν 23 < √ 4π, leading to the following upper bound of z, The two formulas in Eq. (20) were found in the general symmetric from of M −1 , namely they are be found by using Eq. (12) for off-diagonal entries of m ν to determine (M −1 ) ij , then insert them into the diagonal ones. The off-diagonal elements of M −1 are determined as follows: Hence all elements of the matrix M −1 depend on only three complex parameters (M −1 ) ii (11), six parameters µ X,i and (M −1 ) ii are determined as functions of elements of M R . In this work, we will consider all elements of M R are free parameters, namely where all k ij are assumed to be real for simplicity. The ISS relations are valid with at least some |k ij | 1 and detM = 0. In the numerical investigation, m ν is determined from the 3σ neutrino oscillation data through Eq. (12). The Dirac matrix m D is then determined by Eq. (19). The free parameters k ij and z are assumed to be real, and z is positive. The three elements of the matrix µ X are determined as functions of these free parameters. The respective formulas are lengthy hence they are not written down explicitly here. In our work, we only consider the case max|µ X,i | z hence all (µ X ) i gives suppressed mixing elements in the total lepton mixing matrix U ν . This condition will always be checked numerically to derive the final results.
In the numerical investigation, the free parameters z and k ij will be scanned in the valid ranges to construct the total neutrino mass matrix defined in Eq. (7). After that, the mass eigenstates and the total mixing matrix are calculated numerically with at least 30 digits of precision. Using the relations listed in Eqs. (17) and (18), we reproduce all of the oscillation parameters ∆m 2 ij and s 2 ij then force them satisfying the 3σ allowed data. This will help us to collect the allowed values of z and k ij in evaluating the cLFV branching ratios and (g − 2) µ data. We emphasize that the regions of the parameter space in our numerical investigation are more general than those mentioned in Refs. [36,49].
The Lagrangian for quark masses was discussed previously [38]. Here, we just remind the reader that the Yukawa couplings of the top quark must satisfy the perturbative limit Combining this with the relations in Eqs. (4) and (5) gives a lower bound t β ≥ 0.3, which will be used in the numerical discussion.

B. Higgs bosons
The Higgs potential used here respect the new lepton number defined in Ref. [38], namely where f is a dimensionless parameter, which f ω is the same as that used in previous works.
The minimum conditions of the Higgs potential as well as the identification of the SM-like Higgs were discussed in detailed previously [54,60]. The model always contains a light CP even neutral Higgs boson identified with the SM-like Higgs boson confirmed experimentally.
This Higgs boson gives suppressed contributions to (g − 2) µ hence we will ignore it from now on. The model contains two pairs of singly charged Higgs bosons H ± 1,2 and Goldstone bosons of the gauge bosons W ± and Y ± , which are denoted as G ± W and G ± Y , respectively. The masses of all charged Higgs bosons are [51,60,61] where As a consequence, the parameter f must be positive. In addition, f may be small so that charged Higgs boson masses can be around 1 TeV.

III. ANALYTIC FORMULAS FOR ONE LOOP CONTRIBUTIONS TO ∆a ea AND
CLFV DECAYS e b → e a γ All detailed steps for calculation to derive the couplings that give large one-loop contributions were presented in Ref. [49]. We just collect the final results related with this work.
The condition m e b > m ea is always used to define the one loop form factors c X (ab)R and c X (ba)R introduced in Ref. [62], which are different from our notations by a relative factor m e b .
The relevant Lagrangian of charged gauge bosons is corresponding to the following one-loop form factors: where e = √ 4πα em being the electromagnetic coupling constant, and g = e/s W .

Lagrangian of charged Higgs bosons is
where The one-loop form factors are: where b ≥ a, and The total one-loop contribution to the cLFV and ∆a 331ISS The one-loop contributions from charged gauge bosons to the a ea and the electric dipole moment d ea of the charged lepton e a are [62]: The one-loop contribution to a ea and d ea caused by charged Higgs bosons is [62]: The quantity ∆d ea = d V ea + d H ea is the new one loop contributions predicted to the electric dipole moment of the charged leptons. It equals to zero when our investigation is limited in the case of the Dirac phase δ = π. This zero value of d µ satisfies the current experimental constraint [63] hence we will not consider from now on.
We remind the reader that one loop contributions from neutral Higgs bosons are very suppressed hence they are ignored here. The reason is that the 331ISS model has no new charged leptons, hence the one-loop contributions of any neutral Higgs bosons H 0 to c (ab)R must arise only from the couplings H 0ē a e a derived from the first term of the Yukawa Laragian (6). These couplings have the same Yukawa couplings with the SM- but different mixing factors |c H 0 | ≤ 1 telling the contributions of ρ 0 to the physical state H 0 .
Hence these contributions to a µ have the same form with the one from the SM-like Higgs boson having mass m h 125 GeV [64]. Also, the heavy neutral Higgs will give suppressed one-loop contributions to ∆a µ . The deviation of a µ between predictions by the two models 331ISS and SM are where a SM,W µ = 3.887 × 10 −9 [64] is the SM prediction for the one-loop contribution from R . In this work, ∆a 331ISS µ = ∆a µ will be considered as new physics predicted by the 331ISS and will be used to compare with the experimental data in the following numerical investigation. We note that the discrepancy of a e between experiments and SM is about 2.5 standard deviation [3,20,[65][66][67][68]. In this work we will only pay attention to the ∆a µ which is the very interesting result of 4.2 standard deviation and may be a clear signal of new physics in the near future.
Based on Ref. [62], the branching ratios of the cLFV processes are where . This result is consistent with the formulas given used in Refs. [49,61] for 3-3-1 models.
It is noted that for the gauge boson contributions, we have |c V for m e b > m ea . Similarly, we can estimate that |c H,k (ba)R |/|c H,k (ab)R | 1 for every particular contribution. Anyways, in the general case we cannot ignore c X (ba)R because of the situation that when contributions to c (ab)R have the same order but some of them have opposite signs.
Then the very destructive correlations among particular Higgs contributions in the c (ab)R will result in the same order of both |c (ab)R | and |c (ba)R |. This will happen in the 331ISS model and |c (21)R | ≤ O(10 −13 ) [GeV −2 ], respectively. As a result, we can estimate that the one-loop contributions from two charged Higgs bosons to Br(µ → eγ) are strongly destructive, i.e. c H 1 (12) −c H 2 (12) . Simultaneously, |c H k (12) | ∼ |c H k (22) |, therefore the charged Higgs contributions to ∆a µ must be constructive and satisfy |c H 1 (22) , or they can be destructive but |c H i (22) | |c H j (22) | with i = j. These important properties of charged Higgs boson contributions will be the key point in our numerical investigation to collect data points satisfying the large values of ∆a µ ≥ 10 −9 before considering any cLFV decay constraints. The gauge contributions are suppressed hence we do not discuss qualitatively here, but they are also included in the numerical investigation. We just pay attention to the two key one-loop charged Higgs boson contributions which will affect two other cLFV decays τ → eγ, µγ.
The experimental constraints of the form factors c (ab)R are listed in Table I, where the allowed values of ∆a µ are chosen in the range of 1σ confidence level given in Eq. (1). We Br(τ → eγ) |c ( derive that the allowed regions of the parameter space have the following properties: Normally, our numerical scan gives a relation that |c H k As a result, the huge destructive correlation between charged Higgs contributions to guarantee simultaneously the experimental constraints of Br(µ → eγ) and ∆a µ . Also, the two cLFV decays of τ → eγ, µγ also need smaller but still large destructive charged Higgs contributions to satisfy the upper experimental bounds because some of these particular contributions often satisfy |c H k  For convenience in estimating qualitatively the above properties, we define new important quantities determining the correlations between two charged Higgs contributions in a physical process as follows: The first ratio R X ab shows the relative contribution from the particle X in the loop to the total contribution. The second one shows the relative contributions of both singly charged Higgs bosons. In the 331ISS model, we will see that the relations R W ab , R Y ab R H k ab often happens. The interesting possibility we would like to discuss is that large contributions of The appearance of the gauge singlet X R leads to a possibility that, a new singly charged Higgs bosons h ± 3 ∼ (1, 1, ±1) can be included in the 331ISS model so that they can give one-loop contributions to both ∆a µ and cLFV amplitudes through the following Yukawa interactions: The new contributions to the cLFV decays and ∆a 331ISS Although the contributions of these singly charged Higgs bosons to ∆a µ are normally small and negative, the contributions to the cLFV amplitudes may be significantly large. Consequently, they can affect destructively the total cLFV decay amplitudes. These properties will keep ∆a 331ISS µ reaching the experimental constraint given in Eq. (1), while keeping all other cLFV branching ratios well below the experimental constraints. In this work, we consider the simplest case that h ± 3 does not mix with the other singly charged Higgs bosons in the 331RN, and the mass is another free parameter. All of these properties can be derived easily from the total Higgs potential, hence it will be ignored in this work.
Before discussing on the allowed regions that satisfy all experimental constraints of cLFV decays e b → e a γ as well as (g − 2) µ data, we give some important crude estimation on the allowed regions of parameter space constrained by both large ∆a 331ISS µ ≥ O(10 −9 ) and small Br(e b → e a γ). The way to derive the total mass matrix to calculate numerically the masses and total neutrino mixing matrix U ν were presented in the previous section. We have checked that the input changes of ∆m 2 ij and s 2 ij in the allowed ranges given in Eq. (16) do not change significantly the final results, so we will fix these quantities at their best-fit points.
An exception that the Dirac phase δ = 180 [Deg.] is considered so that the imagine parts of c (ab)R are zeros, leading to a simple case of destruction among the one-loop contributions from charged Higgs bosons.
In the numerical scan, the points in the allowed regions also satisfy simultaneously the following conditions: R ] ∼ 0, which will result in valid regions of the parameter space in which two charged Higgs bosons contributions can cancel each others. Therefore, these regions will contain points which give the very suppressed total contributions to guarantee the small Br(e b → e a γ). We will use this condition in our numerical investigation. These allowed regions will be used to collect the allowed points satisfying the remaining cLFV constraints. The following ranges of the parameter space will be chosen as the necessary conditions of free parameters when scanning to collect allowed points:

A crude numerical scan shows that the condition Re
Numerical values of k ij will be chosen so that they give active neutrino masses and U PMNS consistent with neutrino oscillation data. The value of 5.3 TeV is fixed from the lower bound w obtained from the experimental data of the heavy Z boson mass m Z . But it can be relaxed with larger w without any changes of final conclusions in this work.
Without contributions of the additional singly charged Higgs boson, our numerical investigation shows that the largest values of ∆a µ satisfying all cLFV constraints is ∆a µ ≤ 108.5 × 10 −11 , see an illustration shown in Fig. 1. The corresponding ranges of the free parameters are shown in Table II, where the right panel shows the only contributions from   c (ab)R (a < b) to the decay rates, namely Here the two first lines show the respective minimum and maximum values of the free parameters. The third line shows a particular example of the set of the parameters giving large ∆a µ 108.1 × 10 −11 . The other quantities are shown in Table III, which will be discussed more later. In the left panel of Fig. 1, only Br(τ → µγ) always enhances with increasing ∆a µ . The upper constraint Br(τ → µγ) < 4.4 × 10 −8 gives the largest value of ∆a µ 108.5 × 10 −11 . From the right panel of Fig. 1, we see that |c (23)R | < |c (32)R | in the region predicting large ∆a µ , because the contribution from |c (23)R | to Br(τ → µγ) denoted as Br(23R) is small, namely Br(23R) ≤ 0.2 × 10 −8 with ∆a µ ≥ 108 × 10 −11 . This is in contrast to other normal cases, as we will discuss based on the Table III.   given in Table III show  The above properties are also true for the allowed region of the parameter space given in Table II. They are summarized in Table IV through the quantities defined in Eqs. (39) and (40). We can see that   Table II constraints given in Eqs. (1) and (2).
From the above discussion, we can see that max[∆a µ ] 108.5 × 10 −11 predicted by the 331ISS model comes from the experimental constraints Br(τ → µγ) < 4.4 × 10 −8 , which gets main contributions from c (32)R . On the other hand, Br(τ → eγ) can reach to zero for large ∆a µ ≥ 10 −9 , which is different from the normal behave of these two branching ratios After some numerical checks, we see that the this difference is originated mainly from the following property: each quantity Br(τ → µγ) or Br(τ → eγ) contains only one type of terms with a factor mµ mτ or me    show the respective one-loop contributions from only h ± 3 to ∆a µ and Br(τ → µγ), which are defined as follows: The corresponding benchmark is calculated numerically with 30 digits of precision number.
The numerical values of the free parameters are In this case, the heavy neutrino masses are m n 4 = m n 5 = 137.2 GeV, m n 6 = m n 7 = 4709.4 GeV, m n 8 = m n 9 = 11958 GeV. For simplicity we assume that Y 3 11 = Y 3 12 = Y 3 21 = Y 3 13 = Y 3 31 = 0, therefore the contribution from h 3 does not change the two cLFV decays Br(µ → eγ) 3.93 × 10 −13 and Br(τ → eγ) 1.11 × 10 −8 . They always satisfy the experimental data. The non-zero Yukawa couplings are scanned in the ranges Y 3 ab ∈ [−3.5, 3.5] that satisfy the perturbative limit. This results in the following allowed range of the charged Higgs boson mass 500 GeV ≤ m h 3 ≤ 1158GeV. Numerical values of c (ab)R is shown in Table V. The numerical results shown in Fig. 2 have some interesting properties. In the left panel, the contributions from h ± 3 to ∆a µ are always negative, but much smaller than the total one: 0 < −∆a µ (h ± 3 ) ≤ 1.5 × 10 −10 200 × 10 −11 ∼ ∆a µ . On the other hand, the one-loop contributions c h 3 (32)R and c H 2 (32)R have the same order, but opposite signs. Therefore, the total |c (32)R | is small enough to guarantee that Br(τ → µγ) < 4.4×10 −8 . This is reason why in the right panel, we see that |c (32)R | < |c h 3 (32)R |, i.e. Br(τ → µγ) < Br(τ → µγ)[h 3 ] may happen. More specifically, this property can be seen from a particular numerical illustration presented in Table V. We can see a property that |c (22)R | |c h 3 (22)R | ∼ |c h 3 (32)R | ∼ |c H 2 (32)R | ∼ |c (32)R |, which explains why the contributions from h 3 affect strongly Br(τ → µγ) but weakly ∆a µ .
The allowed regions of parameters allowing ∆a 331ISS µ around the value 200 × 10 −11 can be found easily in the ranges given in Eq. (44). The allowed regions with larger ∆a 331ISS µ are shown in Fig. 3 The heavy neutrino masses are in the following ranges: m n 4 = m n