Decays $h\to e_ae_b$, $e_b\to e_a\gamma$, and $(g-2)_{e,\mu}$ in a 3-3-1 model with inverse seesaw neutrinos

We will show that the 3-3-1 model with new heavy right handed neutrinos as $SU(3)_L$ singlets can explain simultaneously the lepton flavor violating decays of the SM-like Higgs boson, charged lepton flavor violating decays $e_b\rightarrow e_a\gamma$, and the electron $(g-2)_e$ anomalies under recent experimental data. The discrepancy of $(g-2)_{\mu}$ predicted by the model under consideration and that of the standard model can reach $10^{-9}$. The decay rates of the standard model-like Higgs boson $h\to \tau e, \tau\mu$ can reach the values of $\mathcal{O}(10^{-4})$.

I.
To explain the neutrino oscillation data, the BSMs with the general seesaw(GSS) mechanism also result in LFV decays. But the versions adding only heavy seesaw neutrinos type-I predict suppressed LFV rates that are much smaller than the upcoming experimental sensitivities [20,21]. In contrast, the models with only new inverse seesaw (ISS) neutrinos can predict large LFV rates. In addition, LFVH rates may be large in the regions satisfying constraints of Br(e b → e a γ) [22][23][24][25]. On the other hand, LFVH rates may be smaller when other constraints are considered [26,27]. In the supersymmetric (SUSY) versions of these models with new LFV sources from superparticles, LFVH rates may reach large order of O(10 −5 ) [20,[28][29][30][31][32][33][34][35][36]. LFVH decays were also addressed with other experimental data in many other non-SUSY extensions of the SM . Many BSM predict that the strong constraints of cLFV decay rates Br(e b → e a γ) give small LFVH ones, or suppressed (g − 2) µ .
Unless there is some specific condition of the appearance of very light new bosons, the above cLFV constraints will result in small new one-loop contributions to anomalous magnetic moments (AMMs) of charged leptons (g − 2) ea /2 ≡ a ea , in contrast with recent experimental data. Namely, the 4.2σ deviation between standard model (SM) prediction [68], combined contributions from previous works , and muon experiments [95,96] is This result is slightly inconsistent with the latest one, which calculated the hadronic vacuum polarization for the SM prediction based on the lattice QCD approach, giving a combined value reported in Refs. [77,78,97] closer to the experimental data. This value was shown to fit with other experimental data such as global electroweak fits [98][99][100].
The recent studies of cLFV decays in the regions satisfying the AMM data were done in some specific models such as SUSY with largest Br(h → τ µ) ∼ O(10 −4 ) [102]. Other BSM containing leptoquarks can explain the large ∆a NP µ ∼ O(10 −9 ) [103]. Recent work has discussed an extension of the 3-3-1 with right-handed neutrinos [104][105][106], named the 3-3-1 model with inverse seesaw neutrinos (331ISS) [107], with the aim of giving an explanation of both the (g −2) µ data and the neutrino oscillation data through the ISS mechanism. The model needs new SU (3) L gauge singlets including three neutral leptons X aR and a new singly charged Higgs boson h ± to accommodate all the experimental data of neutrino oscillation, the cLFV bounds in Eq. (1) and the ∆a µ in 1σ deviation given in Eq. (2). Although cLFV and/or LFVH decays were investigated previously with promoting predictions for the 331ISS [108][109][110][111], the AMM data was not included. Our aim in this work is filling this gap. We note that other 3-3-1 models [112][113][114][115] constructed previously can accommodate the (g − 2) µ data only when they are extended such as adding new vector-like fermions, or/and scalars [116][117][118][119][120]. But none of them paid attention to the correlations between LFVH decays and (g − 2) ea anomalies.
Our paper is organized as follows. In Sec. II we discuss the necessary ingredients of a 331ISS model for studying LFVH decays and how the ISS mechanism works to generate active neutrino masses and mixing consistent with current experimental data. In Sec. III we present all couplings needed to determine the one-loop contributions to the LFVH decay amplitudes of the SM-like Higgs boson, cLFV decays, and (g − 2) ea . In section IV, we provide detailed numerical illustrations and discussions. Section V contains our conclusions.
Finally, the appendix lists all of the analytic formulas expressing one-loop contributions to LFVH decay amplitudes calculated in the unitary gauge.

II. THE 331ISS MODEL FOR TREE-LEVEL NEUTRINO MASSES
A. Particle content and lepton masses We summarize the particle content of the 331ISS model in this section. We ignore the quark sector irrelevant in our work, which was discussed previously [121,122]. We also ignore many detailed calculations presented in Ref. [107]. The electric charge operator defined by the gauge group SU 3 ) and a right-handed charged lepton e aR ∼ (1, −1) with a = 1, 2, 3. The 331ISS model contains three neutral leptons X aR ∼ (1, 0), a = 1, 2, 3, and a singly charged Higgs boson σ ± ∼ (1, ±1). There are three Higgs triplets ρ = (ρ + 1 , ρ 0 , ρ . The vacuum expectation values (vev) for generating all tree-level quark masses and leptons are ρ = (0, v 1 √ 2 , 0) T , η = ( v 2 √ 2 , 0, 0) T and χ = (0, 0, w √ 2 ) T . Two neutral Higgs components have zero vevs because of their non-zero generalized lepton numbers [107] corresponding to a new global symmetry In the 331ISS, nine gauge bosons get masses through the covariant kinetic Lagrangian of the Higgs triplets, ., 8, and T 9 ≡ I 3 √ 6 and 1 √ 6 for (anti)triplets and singlets [123]. There are two pairs of singly charged gauge bosons, denoted as W ± and Y ± , defined as with the respective masses m 2 W = g 2 4 (v 2 1 + v 2 2 ) and m 2 Y = g 2 4 (w 2 + v 2 1 ). The breaking pattern of the model is SU (3) L × U (1) X → SU (2) L × U (1) Y → U (1) Q , leading to the matching condition that W ± are the SM gauge bosons. As a consequence, we have where e and s W are, respectively, the electric charge and sine of the Weinberg angle. Similarly to the Two Higgs Doublet Models (2HDM), we use the parameter which leads to v 1 = vc β and v 2 = vs β .
The Yukawa Lagrangian generating lepton masses is: where a, b = 1, 2, 3. The first term generates charged lepton masses as m ea ≡ h e ab v 1 √ 2 δ ab , with the assumption that the flavor states are also physical.
In the basis n L = (ν L , N L , (X R ) c ) T , Lagrangian in Eq. (7) generates a neutrino mass term written in terms of the total 9 × 9 mass matrix consisting of nine 3 × 3 sub-matrices [109], namely (7) is symmetric, and can be considered as a diagonal matrix without loss of generality.
The following approximation solution of U ν is valid for any specific seesaw mechanisms, where R, V are 3 × 6, 3 × 6 matrices, respectively. All entries of R must satisfy |R aI | 1, so that all ISS relations can be derived perturbatively.
The relations between the flavor and mass eigenstates are where n L ≡ (n 1L , n 2L , ..., n 9L ) T , and the Majorana states are n i = (n iL , n iR ) T .
The ISS relations are From experimental data of m ν , we can determine all independent parameters in m D and [109,121]. Namely, the Dirac mass matrix has the antisymmetric form where α 23 ≡ arg[h ν 32 ],m D is an antisymmetric matrix with (m D ) 23 = 1, and is a positive and real parameter. Eq. (13) gives (m ν ) ij = m T D M −1 m D ij for all i, j = 1, 2, 3, leading to six independent equations. Solving three of them with i = j, the non-diagonal entries of M −1 are functions of M −1 ii and x 12,13 . Inserting these functions into the three remaining relations with i = j, we obtain and T we derive that three parameters of the matrix µ X as certain but lengthy functions of (ze iα 23 ), all entries of M R and m ν . While m ν are fixed by experiments, all entries of M R are free parameters. We will fix α 23 = 0, because it is absorbed into the µ X .
In the limit that |R 2 | 1, the heavy neutrino masses can be determined approximately based on Eq. (14), namely We define the reduced matrix M R ≡ zM R , M R ij ≡ k ij , provided that R * 2 = −m D /M R . The matrix M R is always diagonalized by two unitary transformations V L,R [124]: where allk 1,2,3 are always positive andk a 1 so that all ISS relations are valid. Therefore, M R is expressed in terms ofk and V L,R . Then the matrix V in Eq. (14) can be found approximately as follows As a consequence, for any qualitatively estimations we use the approximation that heavy neutrinos masses are m n a+3 = m n a+6 zk a with a=1,2,3; R 1 O 3 ; and We have checked and confirmed that the above approximations give numerical results consistent with those discussed in Ref. [107]. Therefore, these approximate formulas will be used in this work. m ν is chosen as the input with 3σ neutrino oscillation data to fixm D .
The free parameters z 0 andk 1,2,3 , V R will be scanned in the valid ranges to construct the total neutrino mixing matrix U ν defined in Eq. (21). Because which do not depend explicitly on V L , it affects weakly on all relevant processes. We will fix V L = I 3 from now on.
Lagrangian for quark masses were discussed previously [121,122]. Here, we just recall that the Yukawa couplings of the top quark must satisfy the perturbative limit h u 33 < √ 4π, leading to a lower bound of v 2 : v 2 > √ 2mt √ 4π . Combined with the relations in Eqs. (5) and (6), the lower bound of t β is t β ≥ 0.3. The upper bound of t β can be derived from the tau mass, 4π, leading to a rather weak upper bound

B. Higgs bosons
The Higgs potential used here respect the new lepton number defined in Ref. [122], namely where f is a dimensionless parameter, f η,χ are mass dimensional, S = η, ρ, χ. These three trilinear couplings softly break the general lepton number L. For simplicity, we fix f χ = 0 by applying a suitable discrete symmetry. The last line in Eq. (23) contains all additional terms couplings with new charged Higgs singlets compared with Higgs potential considered in previous works [107]. They do not affect the squared mass matrices of both neutral CPodd and CP-even Higgs bosons. The minimum conditions of the Higgs potential as well as the identification of the SM-like Higgs boson have previously been discussed in detailed [33,125], hence we just list the necessary results here. The model contains three pairs of singly charged Higgs bosons h ± 1,2,3 and two Goldstone bosons G ± W,Y of the singly charged gauge bosons W ± and Y ± , respectively. In the limit of f η = 0, the singly charged Higgs The mass of the Higgs singlet σ ≡ h ± 3 is a function of µ 2 s and λ σ S . With f η = 0 considered in this work, the relations between the original and mass eigenstates of the charged Higgs bosons are where These result are consistent with Refs. [123,125,126] in the limits of s α = 0, ±1. The results given in Eqs. (24) and (25) obtained by solving the following 3 × 3 squared mass matrix in the basis (η ± , ρ ± 1 , σ ± ): We will find out that the Higgs masses m h ± 1,2 and the mixing angle α are functions of the Higgs parameters in the Higgs potential.
The model contains five CP-odd neutral scalar components included in the five neutral Three of them are Goldstone bosons of the neutral gauge bosons Z, Z , and X 0 . The two remaining are physical states with masses As a consequence, the parameter f must satisfies f > 0.
Considering the CP-even scalars, there are two sub-matrices 2 × 2 and 3 × 3 for masses of these Higgs bosons in two bases (η 0 2 , χ 0 1 ) and (η 0 1 , ρ 0 1 , χ 0 1 ), namely The matrix M 2 0,2 has one zero value and m 2 h 4 = f t β +λ 13 2 s 2 β v 2 + ω 2 corresponding to one Goldstone boson of X 0 and a heavy neutral Higgs boson h 0 4 with mass at the SU (3) L breaking scale. On the other hand, we see that Det[M 2 0,3 ] = 0 but Det[M 2 0,3 ] v=0 = 0, which implies that there is at least one Higgs boson mass at the electroweak scale that can be identified with the SM-like Higgs boson. In particular, it can be proved that and C h 1 M 2 0,3 C hT 1 ≡ M 2 0,3 satisfying: Therefore, there is a unitary transformation [109,127,128]. Hence h 0 1 is identified with the SM-like Higgs boson found at the LHC, namely h 0 1 ≡ h. For simplicity, we will fix C h 2 = I 3 in this work, and use the relations This assumption leads to a consequence that the m h 0 1 is independent with the Higgs selfcouplings relating with one-loop decays h 0 1 → e a e b , as it will be seen as follows: where non-zero g h 0 1 ij = g hji are In the next section, we will derive all of the remaining couplings giving one-loop contributions of decays mentioned in this work. III.

COUPLINGS AND ANALYTIC FORMULAS
A. Decays e b → e a γ and (g − 2) ea The couplings of charged gauge bosons giving one-loop contributions to LFV amplitudes are: All the calculation steps to derive theses couplings were presented in Ref. [109]. From now on, we always choose that m e b > m ea , equivalently b > a = 1, 2, 3, to define the decays e b → e a γ. One-loop form factors from charged gauge bosons are [129]: where e = √ 4πα em being the electromagnetic coupling constant; and g = e/s W .
The Yukawa couplings of charged Higgs bosons with leptons are defined by where The interactions given in Eqs. (34) and (37) also give tree and loop contributions to the lepton flavor conserved decay µ − → e − ν e ν µ . Regarding the gauge couplings given in Eq.
(34), the couplings of Y ± with active neutrinos are zeros because U ν (c+3)1 = U ν (c+3)2 = 0, the difference of the couplings of W with active neutrinos and charged leptons between the SM and the 331ISS model under consideration is | 1 2 (R 2 R + 2 U ) ab | 1. Regarding the Higgs boson contributions, only λ L,3 ai may give large contributions to the decay amplitude µ − → e − ν e ν µ , because the remaining couplings are always proportional to ab derived from the second term in Lagrangian (7). Based on the well-known formulas of the partial decay width Γ(µ → 3e) at tree level given in the Zee-Babu model [142], the coupling λ L leads to a deviation of the decay width of the decay µ − → e − ν e ν µ between the 331ISS model and the SM as follows: The constraint is derived from the mean life time of muon [132]. The derivation of the formula (39) is summarized as follows. The total amplitude is iM = iM W + iM h ± , where M W and M h ± are the contributions from W and charged Higgs bosons, respectively. In the low energy limit we have Now it can be proved that has an odd number of the gamma matrices in the trace and m e , m νµ , m νe 0, leading to M * W M h ± = 0. In the numerical investigation, we will choose m h ± 3 ≥ z 0 × 10 √ 5 to accommodate the constraint (39). Now we can assume the approximation that Γ 331ISS (µ − → e − ν e ν µ ) . This approximation for calculating the cLFV decay rates is consistent with many works published recently [140,141].
The one-loop form factors are [129]: where b ≥ a, x k,i = m 2 n i /m 2 h ± k , and the one-loop functions F H (x) andF H (x) are The total one-loop contributions to the cLFV amplitude e b → e a γ and ∆a 331ISS The second line of Eq. (42) is derived from the equality that c (ba) The formulas for the contributions to a ea are: One-loop contributions from heavy neutral Higgs bosons are very suppressed, hence they are ignored here. The deviation of a ea between predictions by the two models 331ISS and SM are where a SM µ (W ) = 5g 2 m 2 µ /(96π 2 m 2 W ) is the SM's prediction [130]. In this work, ∆a ea will be considered as new physics (NP) predicted by the 331ISS, used to compare with experimental data in numerical investigations.
The formulas of U ν given in Eq. (21) results in approximate expressions of c (ab)R and c (ba)R with b ≥ a as follows: where equalities in Eq. (22) were used. In addition, we ignore the minor contributions proportional to R † 2 R 2 , and R 2 R † 2 . Because only two terms relating to R 2 Y σ and R † 2 R 2 depend on V L , but give small one-loop contributions to ∆a ea , we fix V L = I 3 without loss of generality.
The above expressions of c (ab)R and c (ba)R given in Eq. (46)  give c (ab)R ∼m Dm * D . As a result, constraints on cLFV decays always exclude the regions of parameter space predicting large (g − 2) e,µ . This conclusion is completely consistent with the numerical results reported in Ref. [107]. In addition, the presence of σ ± and non-zero Yukawa coupling matrix Y σ is necessary to explain 1σ range of (g − 2) µ obtained by experiment. Additionally, the formulas given in Eq. (46) explain explicitly that large ∆a µ also needs large z 0 . And, large t β and non-zero Y σ support more strong destructive correlations to guarantee that (e b → e a ) satisfies the current constraints.
Finally, We emphasize that the (g − 2) e data and LFVH decays were not discussed previously for the 331ISS model. Our numerical investigation showed that large (g − 2) e requires nonzero values of s α , which was not considered in Ref. [107]. In addition, large values of Y σ 22,33,23,32 should be investigated carefully because they may result in too large Br(h → τ µ) that may be excluded by the experimental constraints.
The Yukawa couplings h 0 1 f f is , namely where is a symmetric coefficient λ 0 ij = λ 0 ji corresponding to the Feynman rules given in Ref. [124]. All of the Feynman rules for couplings involved in LFV processes at one-loop level are listed in Table I, where we used s θ = gv 1 /(2m Y ). We focus on the limit of tiny t θ s θ = 0, and the suppressed deviation of the SM-like Higgs mixing mentioned previously [109,127].
Namely, they will be fixed to be zeros in the numerical calculations.
The effective Lagrangian and partial decay width of the decay h 0 1 → e ± a e ∓ b are L LFVH = h 0 1 ∆ (ab)L e a P L e b + ∆ (ab)R e a P R e b + H.c.,  Here V ± = W ± , Y ± ; k, l = 1, 2, 3.
where the analytic forms of ∆ (i)W (ab)L,R and ∆ (i)Y (ab)L,R are shown in the Appendix. There are numbers of tiny one-loop contributions, which we will ignore in the numerical calculations.
They are calculated using the unitary gauge with the same techniques given in Refs. [25,109].
The contributions from diagrams (2), (3), and (5) with Y ± exchanges have suppressed factors The one-loop contributions from diagram (6) are suppressed with heavy singly charged Higgs bosons, which we checked consistently with the result mentioned in Refs. [67,109].

IV. NUMERICAL DISCUSSION
In this work, we will use the neutrino oscillation data given in Refs. [132,139]. The even when Re[c (ab)R ]=0. Therefore, many very complicated relations between parameters must be satisfied to guarantee that all Im and Re parts contributing to these cLFV decays satisfy the experimental constraints. In this work, the limit that δ = 180 [Deg] is fixed for simplicity.
The mixing matrix V R is parameterized using the same formulas given in Eq. (51), V R = f (s r 12 , s r 13 , s r 23 , 0) with |s r ij | ≤ 1. The remaining free parameters are scanned in the following ranges: and m h ± 3 = 40 TeV, so that the decay width of µ − → e − ν e ν µ is consistent with that predicted by the SM. In addition, the collected points satisfy that max|(R 2 R † 2 ) ab | < 10 −3 with all a, b = 1, 2, 3. This constraints also satisfies many other recent experimental results such as electroweak precision tests, cLFV decays [143][144][145][146]. Br(µ → eν e ν µ ) 1, Br(τ → eν e ν τ ) 0.1782, and Br(τ → µν µ ν τ ) 0.1739. We note that the upper bounds of m h ± 1,2 and t β based on the previous work to accommodate large values of ∆a µ . Chosen scanning range of t β also satisfies the perturbative limit mentioned above.
We comment here the results obtained previously in Ref. [107], where large t β ≥ 50 and small values of singly charged Higgs bosons h ± k (k = 1, 2) are required for large (g − 2) µ satisfying 1σ experimental data of (g − 2) µ and all constraints from cLFV decays e b → e a γ. But only the case of s α = 0 and non-zero Y σ 22,33,23,32 was mentioned. Our numerical investigation shows that this case results in small ∆a e which cannot satisfy 1σ range of the experimental data given in Eq. (3). Without σ ± , we obtain two maximal values of ∆a e that ∆a e ≤ 2.5 × 10 −14 and 1.5 × 10 −14 for the NO and IO schemes, respectively. Hence, determining the regions of parameter space giving large ∆a e will be very interesting.
Because of the above reasons, we focus on the regions of parameter space giving large ∆a e that satisfies the 1σ experimental data of (g − 2) e as well as all current constraints of cLFV decay rates Br(e b → e a γ). The investigation shows that the 1σ range of ∆a e ∈ show that Br(h → τ µ) < O(10 −4 ) will lead to Br(h → eµ) < O(10 −6 ), which is still smaller than the planned experimental sensitivity.
From our numerical investigation, we found that the regions allowing 1σ range of ∆a e data and cLFV constraints are very wide. But the regions allowed large (g − 2) µ are difficult to control. This is because of the large number of free parameters in the 331ISS that our numerical code is not still smart enough to collect these points. Because of the special form ofm D that require the non-degenerate matrixk and the strong destructive correlations between the mixing angles s r ab and the entries of Y σ in order to get small Br(e b → e a γ), in the regions allow large ∆a µ . There may exist some certain relations between these parameters for collecting more interesting points allowing large (g − 2) µ at 1σ experimental range. We will determine them in a future work.
Finally, we comment some properties of the current Z boson decay data which may put useful constraints on the parameter space of the 331ISS model. In the limit of v/w → 0, equivalently t θ = 0, the couplings of Z boson with all other SM particles. We can see that all masses of the new heavy neutrinos appearing in the collected points we showed above as the numerical results are much larger than the Z boson masses. Therefore, Z do have not any new tree level decays Z → n I n j with at least a new heavy neutrino n I (I > 3). In addition, all masses of the new heavy particles predicted by the 331ISS models are heavier than the Z boson masses, therefore the invisible decays of the Z boson in this case is the same as that in the SM and the 2HDM discussed in Ref. [67]. We therefore conclude that the current Z boson decay data affects weakly the allowed region of the parameters space we focus on this work.
There is another cLFV decay mode Z → e + a e − b discussed in detailed in 2HDM [67], which is still invisible in the regions predicting large Br(h → e a e b ) and satisfying all the constraints of cLFV decays Br(e b → e a γ). Therefore, this decay channel will not change significantly the allowed regions of parameters discussed in this work. On the other hand, the interesting topic we will focus on is that when the experimental sensitivities are improved, both cLFV decays of µ − → e − ν e ν µ and Z → e a e b may give more significant constraints on those mentioned in this work.

V. CONCLUSION
In this work, we have constructed the analytic formulas for one-loop contributions to the LFV decays of the SM-like Higgs boson h 0 1 → e a e b in the 331ISS model. We also give analytic formulas to explain qualitatively the results of large (g − 2) µ reported previously.
Numerical tests were used to confirm the consistency between the two calculations. We introduced a new parameterization of the heavy neutrino mass matrix to reduce the number of free parameters used to investigate (g − 2) e,µ anomalies, LFV decays e b → e a γ, and h 0 1 → e a e b . Our numerical investigation shows that the model can predict easily the 1σ range of experimental data for (g − 2) e and satisfy simultaneously the cLFV constrains Br(e b → e a γ). But we only obtained the regions of parameter space that give largest values of ∆a µ 10 −9 , which rather smaller than the lower bound of 1σ range reported recently.
The reason is that the recent numerical code used in our investigation only works in the limit of small max[|y σ |] < 0.25. In these regions of the parameter space, the largest values of Br(h 0 1 → τ e) and Br(h 0 1 → τ µ) are order of O(10 −4 ) and 10 −3 , respectively. In addition, large Br(h 0 1 → τ µ) predicts large Br(h 0 1 → µe) ∼ O(10 −5 ), which is close to the recent experimental bounds. The regions with large |Y σ | ab may be more interesting, which is our future work, where many other LFV processes such as Z → e b e a , e b → e c e d e f , and the µ − e conversion in nuclei will be discussed together. The one-loop contributions here are calculated using the notations of Passarino-Veltman (PV) functions [135,136] given in Ref. [27], consistent with LoopTools [137], see a detailed discussion in Refs. [33,138]. The PV functions used in this work are defined as follows: As mentioned in Ref. [27], the two B (1) 1 and C 1 have opposite signs with those introduced in Ref. [33]. They come from the signs of p 1,2 in the internal momenta (k − p 1 ) and (k + p 2 ) shown in Fig.   1, which p 1 has an opposite sign, which is different from the standard notation of k + p 1 defined in LoopTools. The PV-functions used in our formulas are: B (1) 0 (m 2 n i , m 2 n j ), B (1) 1 = B (1) 1 (m 2 W , m 2 n j ), and C 0,1,2 = C 0,1,2 (m 2 W , m 2 n i , m 2 n j ). The analytic expressions ∆ (ab)L,R with i = 4, 6, 9, 10, are (1) 0 (m 2 n i , m 2 n j ), B (1) 1 = B (1) 1 (m 2 Y , m 2 n j ), and C 0,1,2 = C 0,1,2 (m 2 Y , m 2 n i , m 2 n j ).
where B (1) (1) k (m 2 Y , m 2 n i ) (k = 0, 1) and C 0,1,2 = C 0,1,2 (m 2 0 + B where B (2) k = B −λ R,k * ai λ L,k bi m n i C 0 + λ L,k * ai λ L,k bi m a C 1 + λ R,k * ai λ R,k bi m b C 2 , −λ L,k * ai λ R,k bi m n i C 0 + λ R,k * ai λ R,k bi m a C 1 + λ L,k * ai λ L,k bi m b C 2 , h ± k C 0 + m 2 a C 1 + m 2 b C 2 + λ R,k * ai λ R,k bj m b m n j C 2 + λ L,k * ai λ L,k bj m a m n i C 1 + λ 0 ij λ R,k * ai λ L,k bj m n i m n j C 0 + λ R,k * ai λ R,k bj m n i m b (C 0 + C 2 ) + λ L,k * ai λ L,k bj m a m n j (C 0 + C 1 ) + λ L,k * ai λ R,k bj m a m b (C 0 + C 1 + C 2 ) , λ 0 ij λ L,k * ai λ R,k bj B (12) 0 + m 2 h ± k C 0 + m 2 a C 1 + m 2 b C 2 + λ L,k * ai λ L,k bj m b m n j C 2 + λ R,k * ai λ R,k bj m a m n i C 1 + λ 0 * ij λ L,k * ai λ R,k bj m n i m n j C 0 + λ L,k * ai λ L,k bj m n i m b (C 0 + C 2 ) + λ R,k * ai λ R,k bj m a m n j (C 0 + C 1 ) + λ R,k * ai λ L,k bj m a m b (C 0 + C 1 + C 2 ) , where k = 1, 2, 3, B (1) 0 − m 2 a B (2) 0 +m a m b λ L,k * ai λ L,k bi m b + λ R,k * ai λ R,k bi m a −B (1) were shown in Refs. [25,121], and hence we do not present them in this work. We note that the scalar functions ∆ (1)W L,R and ∆ (1,2,3)Y L,R include parts that do not depend on m n i , and therefore they vanish because of the Glashow-Iliopoulos-Maiani mechanism.
The divergent cancellation in the total ∆ L,R is shown as follows.