Decay of standard model-like Higgs boson $h\rightarrow \mu\tau$ in a 3-3-1 model with inverse seesaw neutrino masses

By adding new gauge singlets of neutral leptons, the improved versions of the 3-3-1 models with right-handed neutrinos have been recently introduced in order to explain recent experimental neutrino oscillation data through the inverse seesaw mechanism. We prove that these models predict promising signals of lepton-flavor-violating decays of the standard-model-like Higgs boson $h^0_1\rightarrow \mu\tau,e\tau$, which are suppressed in the original versions. One-loop contributions to these decay amplitudes are introduced in the unitary gauge. Based on a numerical investigation, we find that the branching ratios of the decays $h^0_1\rightarrow\mu\tau,e\tau$ can reach values of $10^{-5}$ in the regions of parameter space satisfying the current experimental data of the decay $\mu\rightarrow e\gamma$. The value of $10^{-4}$ appears when the Yukawa couplings of leptons are close to the perturbative limit. Some interesting properties of these regions of parameter space are also discussed.

I.

INTRODUCTION
Signals of lepton-flavor-violating decays of the standard-model-like Higgs boson (LFVHDs) were investigated at the LHC [1] not very long after its discovery in 2012 [2]. So far, the most stringent limits on the branching ratios (Br) of these decays are Br(h → µτ, eτ ) < O(10 −3 ), from the CMS Collaboration using data collected at a centerof-mass energy of 13 TeV. The sensitivities of the planned colliders for LFVHD searches are predicted to reach the order of 10 −5 [3].
On the theoretical side, model-independent studies showed that the LFVHDs predicted from models beyond the standard model (BSM) are constrained indirectly from experimental data such as lepton-flavor-violating decays of charged leptons (cLFV) [4]. Namely, they are affected most strongly by the recent experimental bound on Br(µ → eγ). Fortunately, large branching ratios of the decays h → µτ, eτ are still allowed up to the order of 10 −4 .
Also, LFVHDs have been widely investigated in many specific BSM models, where the decay rates were indicated to be close to the upcoming sensitivities of colliders, including nonsupersymmetric [5,6] and supersymmetric versions [7]. Among them, the models based on the gauge symmetry SU (3) C ×SU (3) L ×U (1) X (3-3-1) contain rich lepton-flavor-violating (LFV) sources which may result in interesting cLFV phenomenology such as charged lepton decays e i → e j γ [8][9][10][11]. In particular, it was shown that Br(µ → eγ) is large in these models, and hence it must be taken into account to constrain the parameter space. In addition, such rich LFV resources may give large LFVHD rates as promising signals of new physics.
Although the 3-3-1 models were introduced a long time ago [12,13], LFVHDs have been investigated only in the version with heavy neutral leptons assigned as the third components of lepton (anti) triplets, where active neutrino masses come from effective operators [14].
The largest values of LFVHD rates were shown to be O(10 −5 ), originating from heavy neutrinos and charged Higgs bosons [6]. Improved versions consisting of new neutral lepton singlets were recently introduced [8,15]. They are more interesting because they successfully explain the experimental neutrino data through the inverse seesaw (ISS) mechanism. We call them the 331ISS models for short. They predict a large cLFV decay rate of µ → eγ corresponding to recent experimental bounds. They may also contain dark matter candidates [8,15]. These properties make them much more attractive than the original versions of 3-3-1 models with right-handed neutrinos (331RHN) [13]. They predict suppressed LFV decay rates, because all neutrinos including exotic ones are extremely light. Furthermore, loop corrections to the neutrino mass matrix must be taken into account to obtain an active neutrino mass spectrum that explains the experimental data [16]. Hence, LFV signals are an interesting way to distinguish the 331ISS and 331RHN models. More specifically, a simple ISS extension of the SM allows large Br(h → µτ, eτ ) ∼ O (10 −5 ) in the allowed regions satisfying Br(µ → eγ) < 4.2 × 10 −13 [17]. Inspired by this, we will address the following questions in this work: how large is the Br(h → µτ, eτ ) predicted by the 331ISS models under the experimental constraints of the cLFV decays? and, are these branching ratios larger than the values calculated in the simplest ISS extension of the SM? Because these 331 models contain many more particles that contribute to LFV processes through loop corrections, either constructive or destructive correlations among them will strongly affect the allowed regions of the parameter space satisfying the current bound of the decay rate µ → eγ. The most interesting allowed regions will also allow large LFVHD rates, which we will try to look for in this work. Because the discussion on the decay h → eτ is rather similar to the decay h → µτ , we only briefly mention the later.
Our paper is organized as follows. In Sec. II we discuss the necessary ingredients of a 331ISS model for studying LFVHDs and how the ISS mechanism works to generate active neutrino parameters consistent with current experimental data. In Sec. III we present all couplings needed to determine the one-loop amplitudes of the LFVHDs of the SM-like Higgs boson. Sec. IV we show important numerical LFVHD results predicted by the 331ISS model. Section . V contains our conclusions. Finally, the Appendix lists all of the analytic formulas expressing one-loop contributions calculated in the unitary gauge.

II. THE 331ISS MODEL FOR TREE-LEVEL NEUTRINO MASSES
A. The model and neutrino masses from the inverse seesaw mechanism First, we will consider a 331ISS model based on the original 331RHN model given in Ref. [16], where active neutrino masses and oscillations are generated from the ISS mechanism.
The quark sector and SU (3) C representations are irrelevant in this work, and hence they are omitted here. The electric charge operator corresponding to the gauge group SU ( consists of a SU (3) L triplet ψ aL = (ν a , e a , N a ) T L ∼ (3, − 1 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 implies a new righthanded neutrino beyond the SM. The three Higgs triplets are ρ . The necessary vacuum expectation values for generating all tree-level quark masses are ρ = Gauge bosons in this model get masses through the covariant kinetic term of the Higgs bosons, where the covariant derivative for the electroweak symmetry is defined as and T 9 ≡ I 3 √ 6 and 1 √ 6 for (anti)triplets and singlets [18]. It can be identified that where e and s W are, respectively, the electric charge and sine of the Weinberg angle, s 2 W 0.231.
The model includes two pairs of singly charged gauge bosons, denoted as W ± and Y ± , defined as The bosons W ± are identified with the SM ones, leading to v 2 1 +v 2 2 ≡ v 2 = (246GeV) 2 . In the remainder of the text, we will consider in detail the simple case given in Refs. [6,19].
The two global symmetries-namely normal and new lepton numbers denoted respectively, as L and L were introduced. They are related to each other by [16,20]: The detailed values of nonzero lepton numbers L and L are listed in Table I. All tree-level lepton mass terms come from the following Yukawa part: Fields χ η ρ ψ aL e aR  where ijk is the antisymmetric tensor 123 = 1, (ψ aL ) c ≡ ((ν aL ) c , (e aL ) c , (N aL ) c ) T , and h ν is an antisymmetric matrix, h ν ab = −h ν ba . The first term of Eq. (4) generates charged lepton masses m a satisfying h e ab ≡ √ 2δ ab m a /v 1 in order to avoid LFV processes at the tree level.
The second term in Eq. (4) is expanded as follows: where we have used the equality N aL (ν bL ) c = ν bL (N aL ) c ,... The second term on the left-hand side of Eq. (5) contributes a Dirac neutrino mass term −L ν The model can predict a neutrino mass spectrum that is consistent with current neutrino data [21] when loop corrections are included, where all new neutrinos are very light [16]. As a result, they will give suppressed LFV decay rates. Now we consider a 331ISS model as an extension of the above 331RHN model, where three right-handed neutrinos which are gauge singlets, X aR ∼ (1, 0), a = 1, 2, 3 are added. 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 µ X is a 3×3 symmetric matrix and L(X aR ) = L(X aR ) = −1. The last term in Eq. (6) is the only one that violates both L and L, and hence it can be assumed to be small, which is exactly the case in the ISS models. The first term generates mass for heavy neutrinos, resulting in a large Yukawa coupling Y ab with SU (3) L Higgs triplets. In addition, the ISS mechanism allows for large entries in the Dirac mass matrix m D originated from Eq. (4), which is the opposite of the well-known requirement in the 331RHN model.
A four-component (Dirac) spinor n i is defined as n i ≡ (n iL , (n iL ) c ) T = n c i = (n i ) c , where the chiral components are n L,i ≡ P L n i and n R,i ≡ P R n i = (n L,i ) c with chiral operators P L,R = 1±γ 5 2 . Similarly, the definitions for the original neutrino states are ν a ≡ (ν L,a , (ν L,a ) c ) T , ν a ≡ (N L,a , (N L,a ) c ) T , X I ≡ ((X R,I ) c , X R,I ) T , and ν = (ν, N ) T . The relations in Eq. (10) can be written as follows: 2, ..., 9. (11) In general, U ν is written in the form [23] where O is a 3 × 6 null matrix, and U , V , and Ω are 3 × 3, 6 × 6, and 9 × 9 unitary matrices, respectively. Ω can be formally written as where R is a 3 × 6 matrix with the maximal absolute values for all entries |R| satisfying |R| < 1. The matrix U = U PMNS is the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix [24],  (14) and c ab ≡ cos θ ab , s ab ≡ sin θ ab . The Dirac phase δ and Majorana phases α, β are fixed as δ = π, α = β = 0. In the normal hierarchy scheme, the best-fit values of neutrino oscillation parameters are given as [21] ∆m 2 21 = 7.370 × 10 −5 eV 2 , ∆m 2 = 2.50 × 10 −3 eV 2 , s 2 12 = 0.297, s 2 23 = 0.437, s 2 13 = 0.0214, where ∆m 2 21 = m 2 n 2 − m 2 n 1 and ∆m 2 = m 2 n 3 − Hence, the following seesaw relations are valid [23]: where In the model under consideration, the Dirac neutrino mass matrix m D must be antisymmetric. Equivalently, m D has only three independent parameters x 12 , x 13 , and z, where z = √ 2v 1 h ν 23 . In contrast, the matrix m ν in Eq. (17) is symmetric, (m ν ) ij = (m ν ) ji , implying that with i, j = 1, 2, 3. This means that a symmetric matrix M will give a right antisymmetric matrix m D . To fit the neutrino data, there must exist matrices M and m D that satisfy the first equality in Eq. (17). Here we choose M to be symmetric for simplicity. There must exist some sets of z, x 12 , x 13 , and M ij (i ≤ j ≤ 3) that satisfy the six equations From the three equations corresponding to i = j = 1, 2, 3, we can write (M −1 ) ii as three functions of z, x 12 , x 13 , and (M −1 ) ij (i = j).
Inserting them into the three remaining equalities, and taking some intermediate steps, where we exclude the case of x 12 , x 13 = 0. Solving the above three equations leads for two solutions of x 12,13 and a strict relation among (m ν ) ij : Interestingly, the last relation in Eq. (22) allows us to predict possible values of the unknown neutrino mass based on the identification given in Eq. (17). Using the experimental data given in Eq. (15), we derive that m ν 1 = 0 in the normal hierarchy scheme. The Dirac matrix now only depends on z: 8 The above discussion also gives M = diag (10 10 z 2 , 7.029 × 10 10 z 2 , −2.377 × 10 11 z 2 ) for a diagonal M R . In this work, we also consider the simple case where M R is diagonal and all elements are positive. We also assume that |m ν | < µ X |m D | < |M R |. We then derive that heavy neutrino masses are approximately equal to the entries of M R , as given in Eq. (19). However, this approximation is not good for investigating LFVHDs, where a divergent cancellation in the numerical computation is strictly required. Instead, we will use the numerical solutions of heavy neutrino masses as well as the mixing matrix U ν so that a total divergent part vanishes in the final numerical results. This treatment will avoid unphysical contributions originated from divergent parts.
Another parameterization shown in Ref. [8], can be applied to the general cases of nonzero δ as well as both the inverse and normal hierarchy cases of active neutrino masses. With the aim of finding regions with large LFVHDs, we will choose the simple case of m D given in Eq. (23).
For simplicity in the numerical study, we will consider the diagonal matrix M R in the The parameter k will be fixed at small values that result in large LFVHD effects. The total neutrino mass matrix in Eq. (7) depends on only the free parameter z. The heavy neutrino masses and the matrix U ν can be solved numerically, which is not affected by z because |µ X | z.
Using the exact numerical solutions for the neutrino masses and mixing matrix U ν for our investigation, we emphasize that the masses and mixing parameters of active neutrinos derived from the numerical diagonalization of the matrix M ν given in Eq. (7) should satisfy the 3σ constraint of the experimental data. In contrast, neutrino masses and mixing parameters defining the matrix m ν in Eq. (17), which are used to calculate the matrix m D , are considered as free parameters. In other words, the experimental values of neutrino masses and mixing parameters are only used to estimate the allowed ranges of free parameters determining the mass matrix M ν . After that, it is diagonalized numerically to find the neutrino masses as well as the mixing matrix U ν . The mixing parameters will be calculated from the matrix U PMNS , which is related to U ν by the relation (12). Requiring that the expansion of Ω in Eq. (13) and the ISS condition |µ X | > m n 3 are valid, we find that small values of k > 1 are allowed. In particular, we find that if three mixing parameters are fixed at the three respective center values, the two inputs for the active neutrino masses may be outside of (but very close to) the 3σ ranges with k = 5. When k ≥ 5.5, we always find that the input lies within the 3σ ranges of the experimental data that produces the consistent numerical solutions of active neutrino masses. When k ≥ 9, the input corresponding to all center values given in Eq. (15) always produces numerical solutions lying in the 3σ ranges of experimental data.
The LFVHD rates depend strongly on the unitarity of the mixing matrix U ν and heavy neutrino masses. On the other hand, they are weakly affected by the requirement that solutions for active neutrino masses and mixing parameters satisfy the 3σ experimental data. Hence, we will use the matrix m D given in Eq. (23) and k ≥ 5.5 for our numerical investigation. We numerically checked that our choice produces reasonable values for the neutrino data close to the 3σ ranges mentioned above.

B. Higgs and gauge bosons
To study the LFVHD effects, we will choose the simple case of the Higgs potential discussed in Refs. [6,19], namely, where f is a mass parameter and is assumed to be real. The detailed calculations for finding the masses and the mass eigenstates of Higgs bosons were presented in Refs. [6,19], where the minimum condition results in v 1 = v 2 . Here we will only list the part that is involved in LFVHDs.
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 m 2 where t θ = v 2 /w. The relations between the original and mass eigenstates of the charged The neutral scalars are expanded as There are four physical CP-even Higgs bosons h 0 1,2,3,4 and a Goldstone boson of the non-Hermitian gauge boson. The neutral Higgs components relevant for this work are defined where s α = sin α and c α = cos α, and they are defined by There is one neutral CP-even Higgs boson h 0 1 with a mass proportional to the electroweak scale, The decoupling limit t θ and s α 0 [19], resulting in the couplings similar to those predicted by the SM; see Table II. Hence h 0 1 is identified with the SM-like Higgs boson found at the LHC.

A. Couplings
In this section we present Yukawa couplings in terms of U ν and physical neutrino masses.
From this, amplitudes and the LFVHD rate are formulated in terms of physical masses and mixing parameters. The equality derived from Eq. (9) where a, b = 1, 2, 3, and the sum is taken over k = 1, 2, .., 9.
The relevant couplings in the first term of the Lagrangian (4) are The relevant couplings in the second term of the Lagrangian (4) are where the last line is derived following the calculation in Ref. [22]: The first term in Eq. (6) gives the following couplings: where we have used The LFVHD couplings between leptons and charged gauge bosons (W ± , Y ± ) are where , 2, .., 8) are the Gell-mann matrices, and t = g X /g. The charged gauge bosons are W ± µ = By defining a symmetric coefficient λ 0 ij = λ 0 ji , namely, the coupling h 0 1 n i n j derived from Eqs. (32) and (33) is written in the symmetric form gcα 4m W h 0 1 n i λ 0 ij P L + λ 0 * ij P R n j , which gives the right vertex coupling based on the Feynman rules given in Ref. [25]. The Yukawa couplings of charged Higgs bosons are defined by Finally, all of the couplings involved in LFV processes are listed in Table II. The model pre- dicts that the following couplings are zero:

B. Analytic formulas
The effective Lagrangian of the LFVHDs of the SM-like Higgs boson h 0 1 → e ± a e ∓ b is where the scalar factors ∆ (ab)L,R arise from the loop contributions. In the unitary gauge, the one-loop Feynman diagrams contributing to this LFVHD amplitude are shown in Fig. 1. Here The partial width of the decay is with the condition m h 0 1 m a,b . Where m a,b are the masses of muon and tau, respectively.
The on-shell conditions for external particles are p 2 1,2 = m 2 a,b and p 2 [21,26]. The ∆ (ab)L,R can be written as

The corresponding branching ratio is Br
where the analytic forms of ∆ (ab)L,R are shown in the Appendix. They can be calculated using the unitary gauge with the same techniques given in Refs. [6,22]. We have crosschecked this with FORM [27].
The divergence cancellation in the total amplitude (37) is proved analytically in the Appendix, based on the following strict equality: In the model under consideration, the divergent parts coming from the contributions of charged Higgs and heavy gauge bosons are related to both (M ν * M ν ) (a+3)(b+3) and , which are affected by heavy neutrino masses. The cancellation in the total divergent part requires that the physical heavy neutrino masses and U ν must be the exact values. Hence, approximate forms of the heavy neutrino masses and neutrino mixing matrix derived from the ISS mechanism cannot be applied. In contrast, we checked numerically that these formulas are safely used in the usual minimal ISS version extended directly from the SM, because the divergent parts are only involved with the elements (M ν * M ν ) ab = (m * D m T D ) (ab) . Many of the contributions listed in Eq. (37) are suppressed, and hence they can be ignored in our numerical computation. From now on, we just focus on the decay h 0 1 → µτ , and hence the simplified notations ∆ L,R ≡ ∆ (23)L,R will be used. The decay h 0 1 → eτ has similar properties, so we do not need to discuss it more explicitly. We can see that | ∆ L ∆ R | O mµ mτ . In addition, we prove in the Appendix that the following combinations are finite: ∆ 1 , B The regions of parameter space predicting large branching ratios for LFVHDs are affected strongly by the current experimental bound Br(µ → eγ) < 4.2 × 10 −13 [28]. A very good approximate formula for this decay rate in the limit m µ , m e → 0 is [11] Br where G F = g 2 /(4 √ 2m 2 W ) and D R is the one-loop contribution from charged gauge and Higgs boson mediations, R . The analytic forms are Because all charged Higgs bosons couple with heavy neutrinos through the Yukawa coupling matrix h ν ab , this matrix is strongly affected by the upper bound O(10 −13 ) on Br(µ → eγ). In fact, our numerical investigation shows that the allowed regions with light charged Higgs masses are very narrow. The previous investigation in Ref. [8] showed that the 331ISS models predicts a large Br(µ → eγ), where the allowed regions discussed there were chosen such that k ∼ O(10 3 ) and M R ≤ 1 TeV, implying that z ∼ O(1) eV. We checked that our formulas are consistent with these results. In general, the allowed regions are very strict, and satisfy one of the following conditions. First, the regions have a small z and large |M R | and m H ± 2 , implying k 1, including those mainly discussed in Ref. [8]. Second, the regions allow for a large m D and small k, but the strong destructive correlation between the two-loop contributions of charged gauge and Higgs bosons must happen. These regions were also considered in Ref. [8], but they were not given much attention. They are very interesting because they predict large branching ratios for LFVHDs and light particles such as new neutrinos and charged Higgs bosons, which could be found at the LHC and planned colliders [29,30]. Hence, our numerical investigation will focus on this case. IV.

A. Setup parameters
To numerically investigate the LFVHDs of the SM-like Higgs boson, we will use the following well-known experimental parameters [21]: Other parameters can be calculated in terms of the above free ones, namely, Apart from that, the mixing parameter α of the neutral CP-even Higgs was defined in Eq. (28). The Higgs self-coupling λ 2 is determined as [6] In the model under consideration with the quark sector given in Refs. [16,29], only the charged Higgs bosons H ± 2 couple with all SM leptons and quarks. They have been investigated at the LHC in the direct production pp → t(b)H ± , which then decay into two final fermion states [31]. But the specific constraints on them in the framework of the 3-3-1 models have not been mentioned yet, to the best of our knowledge. Instead, the lower bounds on their masses have been discussed recently based on recent data of neutral meson mixing B 0 −B 0 , where a reasonable lower bound of m H ± 2 ≥ 480 GeV was concerned [29]. The values of λ 1,2,12 must satisfy theoretical conditions of unitarity and the Higgs potential must be bounded from below, as mentioned in Ref. [6]. The heavy charged gauge boson mass m Y is related to the recent lower constraint of neutral gauge boson Z in this model.
For the above reasons, the default values of the free parameters chosen for our numerical investigation are as follows. Without loss of generality, the Higgs self-couplings are fixed as λ 1 = 1, λ 12 = −1, which also guarantee that all couplings of the SM-like Higgs boson approach the SM limit when t θ → 0. The default value m Y = 4.5 TeV satisfies all recent constraints [29,32]. The parameter z will be considered in the range of the perturbative limit z < 2 √ π × v 1 617 GeV: in particular, we will fix z = 50, 200, 400, 500, and 600 [GeV]. Finally, the charged Higgs mass m H ± 2 will be investigated mainly in the range of 300 to 5 × 10 4 GeV, where large values of LFVHDs may appear.

B. Numerical results
First, we reproduce the regions mentioned in Ref. [8], where M R was chosen to be from hundreds of GeV to 1 TeV, and the scale of m D (namely, z), was a few GeV, corresponding to k 1. As a result, the respective regions of parameter space always satisfy the experimental bound on Br(µ → eγ) with large enough m H ± 2 . These regions are shown in Fig. 2 with fixed z = 1, 5, 10, 100, and 500 GeV. All allowed regions (i.e, those that satisfy the upper bound Br(µ → eγ) < 4.2 × 10 −13 ) give a small Br(h 0 1 → µτ ) < O(10 −9 ). In general, for larger k we checked numerically that the values of the branching ratio of LFVHDs will decrease significantly, and hence we will not discuss this further.
With small values of k = 5.5 and 9, the dependence of both Br(µ → eγ) and Br(h → µτ ) on m H ± 2 with fixed z are shown in Fig. 3. Most regions of the parameter space are ruled out by the bound on Br(µ → eγ), except for narrow parts where particular contributions from charged Higgs and gauge bosons are destructive. This interesting property of the 331ISS model was indicated previously in Ref. [8]. Furthermore, it predicts allowed regions that give a large Br(h 0 1 → µτ ). In particular, the largest values can reach O(10 −4 ) when k = 5.5 and z = 600 GeV, which is very close to the perturbative limit. In general, the illustrations in two Figs. 2 and 3 suggest that this branching ratio is enhanced significantly for smaller k and larger z, but changes slowly with the change of large m H ± 2 . In contrast, small m H ± 2 plays a very important role in creating allowed regions that predict a large LFVHD. Br(µ → eγ) does not depend on m H ± 2 when it is large enough. Furthermore, the branching ratio decreases with increasing k and it will go below the experimental bound if k is large enough.
The allowed regions in Fig. 3 are shown more explicitly in Fig. 4, corresponding to k = 5.5 and k = 9. Only regions that give a large Br(h 0 1 → µτ ) are mentioned. They are bounded between two black curves representing the constant value of Br(µ → eγ) × 10 13 = 4. Clearly, Br(h 0 1 → µτ ) is sensitive to z and k, while it changes slowly with changing values of m H ± 2 . In contrast, the suppressed Br(µ → eγ) allows narrow regions of the parameter space, where some particular relation between m H ± 2 and k and z is realized. Hence, if these two decay channel are discovered by experiments, depending on their specific values, a relation between heavy neutrino and charged Higgs masses can be determined from the 331ISS framework.

V. CONCLUSION
The 331ISS models seem to be the most interesting among the well-known 331 models because of their rich phenomenology, as indicated in many recent works. This work addressed a more attractive property, namely, the LFVHDs of the SM-like Higgs boson which are being investigated at the LHC. Assuming the absence of the tree-level decays h 0 1 → e a e b and e j → e i γ (j > i), the analytical formulas at the one-loop level to calculate these decay rates in the 331ISS model have been introduced. The divergent cancellation in the total decay amplitudes of h 0 1 → e a e b was shown explicitly. From the numerical investigation, we have indicated that the Br(h 0 1 → µτ ) predicted by the 331ISS model can reach large values of O(10 −5 ). They are even very close to 10 −4 , for example, in the special case with k = 5.5 and z 600 GeV, which is close to the perturbative limit of the lepton Yukawa couplings.
This value is larger than that predicted by the simplest ISS version extended directly from the SM [17]. New charged Higgs bosons may give large contributions to both of the decay rates Br(µ → eγ) and Br(h 0 1 → µτ ), leading to either constructive or destructive correlations with those from the charged gauge bosons. As a by-product, the recent experimental bound on Br(µ → eγ) rules out most of the regions of parameter space with small k and large z, except the narrow regions arising from the destructive correlations between contributions of charged Higgs and gauge bosons. We have shown numerically that only these regions give large Br(h 0 1 → µτ ) > 10 −5 when 400 GeV < z < 600 GeV and k ≤ 9 in the case where the Majorana mass matrix M R is proportional to the identity one. Furthermore, these large values of Br(h 0 1 → µτ ) depend weakly on the masses of the heavy charged gauge bosons, but they require the heavy neutrino mass scale M R and m H ± 2 to be a few TeV, which can be detected at current colliders. Besides, Br(h 0 1 → τ e) has the same result. In conclusion, large branching ratios of the LFV processes like h 0 1 → µτ, eτ will support the 331ISS model and may rule out the original 331RHN model containing only very light exotic neutrinos.
Additionally, many properties of heavy neutrinos and charged Higgs bosons in the 331ISS framework may be determined independently at the SU (3) L scale.
where i = 1, 2. In addition, D = 4 − 2 ≤ 4 is the dimension of the integral, M 0 , M 1 , M 2 are masses of virtual particles in the loop, and µ is an arbitrary mass parameter introduced via dimensional regularization [33]. The external momenta of the final leptons shown in 0 , C 0 , and C 1,2 were shown in Refs. [6,22,34], and hence we will not repeat them here. These functions are used for our numerical investigation. We stress that they were checked numerically to be well consistent with the exact results computed by LoopTooLS [35], as reported in Ref. [36].

The analytic expressions for ∆
Defining ∆ (ab)L,R with i = 4, 6, 9, 10, the analytic expressions for ∆ (2) 0 +m a m b λ L,k * ai λ L,k bi m b + λ R,k * ai λ R,k bi m a B (1) 0 − m 2 a B (2) 0 +m a m b λ R,k * ai λ R,k bi m b + λ L,k * ai λ L,k bi m a B where f 1 = c 2 θ and f 2 = 1/2. The details to derive the expressions in Eq. (A2) are the same as those shown in Refs. [6,22], and hence we do not present them in this work. We note that the scalar functions ∆  U ν * ai U ν bj λ 0 * ij m n j + U ν * (a+3)i λ 0 * ij λ L,1 bj , div ∆ Using the equalities M ν = U ν * M ν U ν † and Eq. (38), we can prove that div ∆ where we have used the antisymmetric property of m D : m T D = −m D . From this, it can be seen that div ∆