Contribution of heavy neutrinos to decay of standard-model-like Higgs boson h→ μτ in a 3-3-1 model with additional gauge singlets

Abstract In the framework of the improved version of the 3-3-1 models with right-handed neutrinos, which is added to the Majorana neutrinos as new gauge singlets, the recent experimental neutrino oscillation data is completely explained through the inverse seesaw mechanism. We show that the major contributions to Br(μ → eγ) are derived from corrections at 1-loop order of heavy neutrinos and bosons. But, these contributions are sometime mutual destructive, creating regions of parametric spaces where the experimental limits of Br(μ → eγ) are satisfied. In these regions, we find that Br(τ → μγ) can achieve values of 10−10 and Br(τ → eγ) may even reach values of 10−9 very close to the upper bound of the current experimental limits. Those are ideal areas to study lepton-flavor-violating decays of the standard modellike Higgs boson (h1). We also pointed out that the contributions of heavy neutrinos play an important role to change Br(h1 → μτ), this is presented through different forms of mass mixing matrices (MR) of heavy neutrinos. When MR ∼ diag(1, 1, 1), Br(h1 → μτ) can get a greater value than the cases MR ∼ diag(1, 2, 3) and MR ∼ diag(3, 2, 1) and the largest that Br(h1 → μτ) can reach is very close 10−3.

I.

INTRODUCTION
The current experimental data has demonstrated that the neutrinos are masses and flavor oscillations [1]. This lead to a consequence that lepton-flavor-violating decays of charged leptons (cLFV) must be existence and it is strongly dependent on the flavor neutrinos oscillation. The cLFV is concerned in both theory and experiment. On the theoretical side, the processes l a → l b γ are loop induced, we therefore pay attention both neutrino and boson contributions, with special attention to the latter because the former is included active and exotic neutrinos, which can be very good solved through mixing matrix [2]. From the experimental side, branching ratio (Br) of cLFV have upper bound as given in Ref. [3].
Br(µ → eγ) < 4.2 × 10 −13 , Br(τ → eγ) < 3.3 × 10 −8 , After the discovery of the Higgs boson in 2012 [4,5], lepton-flavor-violating decays of the standard model-like Higgs boson (LFVHDs) are getting more attention in the models beyond the standard model (BSM). The parameter space regions predicted from BSM for the large signal of LFVHDs is limited directly from both the experimental data and theory of cLFV [5][6][7]. Branching ratio of LFVHDs, such as (h 0 1 → µτ, τ e), are stringent limited by the CMS Collaboration using data collected at a center-of-mass energy of 13TeV, as given Br(h 0 1 → µτ, τ e) ≤ O(10 −3 ). Some published results show that Br(h 0 1 → µτ ) can reach values of O(10 −4 ) in supersymmetric and non-supersymmetric models [8,9]. In fact, the main contributions to Br(h 0 1 → µτ ) come from heavy neutrinos. If those contributions are minor or destructive, the Br(h 0 1 → µτ ) in a model is only about O(10 −9 ) [10]. With the addition of heavy neutrinos, the models have many interesting features, for example, besides creating large lepton flavor violating sources, it can also explain the masses and mixing of neutrinos through the inverse seesaw (ISS) mechanism [11][12][13][14].
We recall that the 3-3-1 models with multiple sources of lepton flavor violating couplings have been introduced long ago [15,16]. With the β parameter, the general 3-3-1 model is separated into different layers, highlighting the properties of the neutral Higgs through contributions to rare decays that can be confirmed by experimental data [17][18][19][20][21][22][23][24][25]. However, LFVHDs have been mentioned only in the version with heavy neutral leptons assigned as the third components of lepton (anti)triplets, where active neutrino masses are generated by effective operators [26,27]. An other way of investigating LFVHDs is derived from the main contributions of gauge and Higgs bosons. Unfortunately, these works can only show that the largest values of Br(h 0 1 → µτ ) is O(10 −5 ) [28,29]. Recently, the 3-3-1 model with right-handed neutrinos (331RHN) is added heavy neutrinos which are gauge singlets (331ISS) has shown very well results when investigating LFVHDs [30][31][32]. As a result, the above model has predicted the parameter space regions, where satisfying the experimental upper limit of the Br(µ → eγ) and Br(h 0 1 → µτ ) can reach value of O(10 −5 ) [32]. In contrast, 331ISS model still has some questions to be solved, such as: in the parameter space regions satisfying the experimental limits of Br(µ → eγ), are Br(τ → eγ) and Br(τ → µγ) excluded? What are the contributions of neutrinos, gauge and Higgs bosons to the LFVHDs? How does the parameterization of neutrinos mixing matrices (both active and exotic neutrinos) affect LFVHDs? In this work, we will solve those problems.
The paper is organized as follows. In the next section, we review the model and give masses spectrum of gauge and Higgs bosons. We then show the masses spectrum of the neutrinos through the inverse seesaw mechanism in Section III. We calculate the Feynman rules and analytic formulas for cLFV and LFVHDs in Section IV. Numerical results are discussed in Section.V. Conclusions are in Section VI. Finally, we provide Appendix A, B to calculate and exclude divergence in the amplitude of LFVHDs.

A. Particle content
We now consider 331ISS model, which is structured from the original 331RHN model as given in Ref. [16] and additional heavy Majorana neutrinos. The electric charge operator corresponding to the electroweak group SU (3) L ⊗U (1) X is Q = T 3 +βT 8 +X, where β = − 1 where U L and D αL for α = 1, 2 are three up-and down-type quark components in the flavor basis, while N c aL ∼ = N aR are right-handed neutrinos added in bottom of lepton triplets.
The scalar sector consists of a triplet field χ, which provides the masses to the new heavy fermions, and two triplets ρ and η, which give masses to the SM fermions at the electroweak scale. These scalar fields are assigned to the following (SU (3)C; SU (3)L; U (1)X) represen- There are two triplets (η, χ) have the same quantum numbers and different from a remain.
Neutral components of scalar triplets are shown relevantly with real and pseudo scalars as: The electroweak symmetry breaking (EWSB) mechanism follows where the vacuum expectation values (VEVs) satisfy the hierarchy ω u, v as done in Refs. [19,33]. In order to generate heavy neutrino masses at tree-level and arise mixing angles, we use the ISS mechanism. Thus, the 331RHN model is extended, where three right-handed neutrinos which are singlets of SU (3) L , F aR ∼ (1, 1, 0), a = 1; 2; 3 are added [32]. This model exits two global symmetries, one noted L and L are the normal and new lepton numbers, respectively. They are related to each other by L = 4 √ 3 T 8 + L [16,34]. Requiring L must be softly broken, one add a Lagrangian terms, which is relevant to F a fields. The general Lagrangian Yukawa relates to leptons and heavy neutrinos is given follows: Two of first term in Eq.(5) generate masses for original charged leptons and neutrinos. The next term describes mixing between N a and F a , and the fouth term generates masses for Majorana neutrinos F a .
To use the simple Higgs spectra, we will choose the case of the Higgs potential discussed in Refs. [28,35], namely, where λ 1 , λ 2 , λ 3 are the Higgs self-coupling constants, f is a mass parameter and is imposed as real.

B. Gauge and Higgs bosons
In the model under consider, we denote g and g as the coupling constants of the electroweak symmetry (SU (3) L ⊗ U (1) X ). Then, the relation between coupling constants and sine of the Weinberg angle following: Gauge bosons in this model get masses through the covariant kinetic term of the Higgs bosons, The model comprises two pairs of singly charged gauge bosons, denoted as W ± and Y ± , defined as The bosons W ± as first line in Eq. (8) are identified with the SM ones, leading to In the remainder of the text, we will consider in detail the simple case [28,29,35]. Under these imposing conditions, we get h 0 1 mixed with the three original states. From the potential Higgs was given in Eq. (6), one can find the masses and the mass eigenstates of Higgs bosons. There are two pair of charged Higgs H ± 1,2 and Goldstone bosons of W ± and Y ± , which are denoted as G ± W and G ± Y , respectively. The relations between the original and physics states of the charged Higgs bosons are : with masses With components of scalar fields are constructed as Eq.(4), the model contains four physical CP-even Higgs bosons h 0 1;2;3;4 and a Goldstone boson of the non-Hermitian gauge boson (X 0 ). The mixing of neutral Higgs h 0 4 and Goldstone of boson X 0 depends on γ angle, (t γ ≡ tan γ = u ω ). Three CP-even Higgs h 0 1;2;3 mutually mix and relate to their original components as: 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 and is identified with SM-like Higgs boson.
The remaining two neutral Higgs (h 0 2 , h 0 3 ) in Eq. (11) have masses on the electroweak symmetry breaking scale (m h 0 2,3 ∼ ω), which are outside the range of LFVHDs so they are not given here.

III. NEUTRINOS MASSES AND ISS MECHANISM
We now consider the Yukawa Lagrangian in Eq. The second term in Eq. (5) is expanded by: From the last term of Eq. (14), using antisymmetric properties of h ν ab matrix and equality N aL (ν bL ) c = ν bL (N aL ) c , we can contribute a Dirac neutrino mass term −L ν mass = ν L m D N R + H.c., with basic, The third term in Eq.(5) generates mass for heavy neutrinos, this consequence comes from the large value of Yukawa coupling Y ab . To describe mixing N a and F a , (M R ) ab = Y ab ω √ 2 is introduced. The last term in Eq.(5) violates both L and L, and hence µ F can be assumed to be small, in the scale of ISS models.
In the basis In the normal seesaw form, M ν can be written: To obtain masses eigenvalue and physics states of neutrinos, one can diagonal M ν by 9 × 9 matrix U ν : The relations between the flavor and mass eigenstates are P L n i = n iL = U ν * ij n jL , and P R n i = n iR = U ν ij n jR , i, j = 1, 2, ..., 9.
In general, U ν is written in the form [36] The matrix U = U PMNS is the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix [37,38],  (22) and c ab ≡ cos θ ab , s ab ≡ sin θ ab . The Dirac phase (δ) and Majorana phases (σ 1 , σ 2 ) are fixed as δ = π, σ 1 = σ 2 = 0. In the normal hierarchy scheme, the best-fit values of neutrino oscillation parameters which satisfied the 3σ allowed values are given as [1] s 2 12 = 0.32, s 2 23 = 0.551, s 2 13 = 0.0216, where ∆m 2 21 = m 2 n 2 − m 2 n 1 and ∆m 2 32 = m 2 n 3 − m 2 n 2 . Hence, the following seesaw relations are valid [36]: where In framework of 331RHN model adding F a as flavor singlets, the Dirac neutrino mass matrix m D must be antisymmetric. From results in Ref. [30], with the aim of giving regions of parameter space with large LFVHDs, m D can be chosen in form to suit for both the inverse and normal hierarchy cases of active neutrino masses as: Therefore, m D has only three independent parameters x 13 , x 23 , and = √ 2vh ν 23 . In general, the matrix m ν in Eq. (25) is symmetric, (m ν ) ij = (m ν ) ji , their components are given by: We can calculate in detail, From Eq.(30), we have two solutions x 13 , x 23 and one equation, which express the relation of components of matrix m ν : Based on experimental data of neutrinos oscillation in Eq.(23), the matrix m D is parameterized and only depends on .
It should be emphasized that m D has a form like Eq.(32) and differ from Ref. [32]. This is caused Eq.(30) can have many different solutions, but the use of Eq.(32) is suited very well to investigate LFVHDs.

IV. COUPLINGS AND ANALYTIC FORMULAS
In this section, we will calculate amplitudes and branching ratios of the LFVHDs in term of M ν and physical neutrino masses. With this aim, all vertices are presented in term of physical masses and mixing parameters. From relation in Eq. (17), one can derive equation as follow: Here a, b = 1, 2, 3 and m n k is mass of neutrino n k , with k run taken over 1, 2, ..., 9. It is interesting that, the relation in Eq.(33) lead to represent Yukawa couplings in term of M ν and physical neutrino masses.
We then pay attention the relevant couplings of LFVHDs. These couplings are derived by Lagrangian Yukawa, Lagrangian kinetics of lepton (or scalar) fields and Higgs potential.
From the first term in Eq.(5), we can give couplings between leptons and Higgs boson as follow: The relevant couplings in the second term of the Lagrangian in Eq. (5) are The couplings get contributions of M R matrix given by: Neutrinos interact with gauge bosons based on the kinetic terms of the leptons. When we only concern with the couplings of the charged gauge bosons, the results are.
To calculate h 0 1 n i n i coupling, we use results in Eqs. (36,38) as given above. Furthermore, we can define symmetry coefficient λ 0 ij = λ 0 ji as Ref. [32], the result therefore obtained.
This result is coincidence with the Feynman rules given in Ref. [39]. In such way, the h 0 1 n k n j coupling can be written in the symmetric form h 0 1 n k n j ∼ h 0 1 n k (λ 0 kj P L +λ * 0 jk P R )n j . For brevity, we also define the coefficients related to the interaction of charged Higgs and fermions as follows: The couplings related to LFVHDs are given in Tab.(I). Especially, based on the characteristics of this model, some couplings of h 0 1 such as h 0  The regions of parameter space predicting large branching ratios (Br) for LFVHDs are affected strongly by the current experimental bound of Br(l a → l b γ), with (a > b) [40]. Therefore, we will simultaneously investigate the LFV decay of charged leptons and LFVHDs. In the limit m a,b → 0, where m a,b are denoted for the masses of charged leptons l a,b , respectively, we can derive the result, which is a very good approximate formula for branching ratio of cLFV as given in Ref. [2] Br(l a → l b γ) The analytic forms are represented as: , , Some quantities are defined as below. Especially, the F (x) function is derived from the characteristics of PV functions in this model.
To investigate the LFVHDs of the SM-like Higgs boson h 0 1 → l ± a l ∓ b , we use scalar factors ∆ (ab)L and ∆ (ab)R , which was first obtained in Ref. [32]. Therefore, the effective Lagrangian of these decays is L eff The partial width of h 0 With conditions p 2 where the analytic forms of ∆ Based on the finite terms ∆ 1,L,R , ∆ 2,L,R , ∆ 3,L,R (∆ i , i = 1, 3 -for short), we can investigate the change in total amplitude of h 0 1 → l a l b with the mass of the heavy neutrinos and other parameters of the model. Higgs m H ± 1 . Therefore, the dependent parameters are given follows.
The heavy charged gauge boson mass m Y is related to the lower constraint of neutral gauge boson Z in 3-3-1 models, which have also been mentioned in Refs. [43,44]. To satisfy those constraints, we choose the default value m Y = 4.5 TeV. The values of the Higgs self-couplings must satisfy theoretical conditions of unitarity and the Higgs potential must be bounded from below, which also guarantee that all couplings of the SM-like Higgs boson approach the SM limit when v ω. For the above reasons, the Higgs self-couplings are fixed as Based on recent data of neutral meson mixing B 0 − B 0 [17], we can choose the lower bound of m H ± 1 ≥ 500GeV. This is also consistent with Ref. [32]. Characteristic for the scale of the matrix m D is the parameter as shown in Eq. (32), considered in the range of the perturbative limit, = √ 2vh ν 23 ≤ 617 GeV. In the calculations below, we fix the values for to be: 100, 200, 400, 500 and 600 GeV. To represent masses of heavy neutrinos (F a ), we parameterize the matrix M R in the form of a diagonal. In particular, the hierarchy of a diagonal matrix M R can yield large results for the LFVHDs.

B. Numerical results of cLFV
In this section, we numerically investigate of l a → l b γ decays with a > b use a diagonal and non-hierarchical M R matrix. That means M R ∼ diag(1, 1, 1). We overhaul regions mentioned in Ref. [32], where was chosen to be from a hundred of GeV to 600GeV, and In a similar way, we can investigate the contributions of gauge and Higgs boson to Br(τ → eγ) and Br(τ → µγ) according to change of m H ± 1 and find out the regions of parametter space which comply with current experimental limits Br(τ → eγ) < 3.3 × 10 −8 and Br(τ → µγ) < 4.4 × 10 −8 [3]. However, the parameter space is only really meaningful when all the experimental limits are satisfied. For the above reasons, we will examine Br(τ → eγ) and Br(τ → µγ) in narrow space regions allowed to satisfy the experimental limit of Br(µ → eγ).
We choose k = 9 and fix = 200, 400, 600GeV, the range of m H ± 1 is from 500GeV to 10TeV,  Br(la →lb γ)  . We illustrate how Br(l a → l b γ) change with m H ± 1 , in the case k = 9 and the fixed values of = 200, 400, 600GeV, corresponding to the plots in the first row of Fig.4. Here, we obtain narrow parameter spaces that satisfy the experimental limits of the Br(µ → eγ), respectively, shown in the second row. The expected space (colorless) is between the two curves 4.2, the remain rules out of the experimental limit (green). It should be emphasized that, for other fixed values of within the limits of the perturbation theory, we can also investigate in the same way.
In each allowed narrow space, where the Br(µ → eγ) is within the experimental limits, Br(τ → eγ) and Br(τ → µγ) also satisfy the upper bound limits of the experiment. In particular, thee values of Br(τ → eγ) can reach as high as 10 −9 and Br(τ → µγ) is about 10 −10 , close to the accuracy found in today's large accelerators. These results are shown in Tab  The l a → l b γ processes have also been studied previously in the context of the 331 models such as Ref. [45]. According to the result, the parameter space areas satisfy the experimental limit (comply with Refs. [46,47]) of the l a → l b γ are given. However, it has two restrictions: i) the limit of µ → eγ is not tight (2.4×10 −12 ), ii) has not shown the region of the parameter space suitable for all l a → l b γ decay. These restrictions have been overcome as shown in Tab.II. This is very interesting result given in the framework of this model and a suggestion for the verification of physical effects in the model from current experimental data.

C. Numerical results of LFVHD
We consider the narrow spatial regions where the Br(l a → l b γ) approach the experimental upper limit, this may be predicting large LFVHD. Therefore, we will investigate the contributions of ∆ i,R,L , i = 1, 3 in Eq. (45) to Br(h 0 1 → l a l b ) and then, we continue to examine Br(h 0 1 → l a l b ) in the narrow spaces mentioned above. This is done both in the case of hierarchical and non-hierarchical M R .
All contributions to Br(h 0 1 → µτ ) in case M R = 9 diag(1, 1, 1) are presented in Fig.6. Here, we chose the fixed values of = 100, 200, 400, 500, 600GeV and m H ± 1 in the range 500GeV to 100TeV. As the result, Br(h 0 1 → µτ ) increases with value of and changes very slowly with the change of large m H ± 1 . The maximum value that Br(h 0 1 → µτ ) can reach is about 0.71 × 10 −3 . This result is very close upper limit of current experimental data Ref. [3].
We use free parameters derived from M R and m D matrices to give the results above.
It is necessary to emphasize the difference with Ref. [32] in parameterizing the matrix m D , we parameterize the matrix m D in the same form as Eq. (32). Using that consequence and including all contributions, especially heavy neutrinos, we have shown that Br(h 0 1 → µτ ) is close to 10 −3 .
The change rule of ∆ i,R,L , i = 1, 3 in Fig.9 produces the survey results of Br(h 0 1 → µτ ) as shown in Fig.10. The largest value of Br(h 0 1 → µτ ) as presented in the right panel of Fig.10 is about 0.160 × 10 −3 . This value is approximately to corresponding ones in case M R = 9 × diag(3, 2, 1), but also smaller when M R = 9 × diag(1, 1, 1). In fact, when M R is chosen in different diagonal form, the heavy neutrinos (N a , F a ) have different masses. This is caused that the contributions of charged Higgs and gauge boson will be destructive interference at the different m H ± 1 . These are consequences that the narrow regions of parameter space where satisfy the experimental limits of Br(l a → l b γ) have different ranges of m H ± 1 as shown in right panels of Fig.6, Fig.8, Fig.10. However, the values of Br(h 0 1 → µτ ) are only really meaningful when considered in narrow spaces that satisfy the experimental limits of Br(l a → l b γ) . These allowed spaces are confined to the two curves 4.2 (black) on the right part of Fig.6, Fig.8  Performing numerical investigation, we point out that ∆ 1,2 are main contributions to Br(h 0 1 → µτ ) while ∆ 3 is ignored because it is very small compared to ∆ 1,2 . All of these contributions are less than 10 −3 in the selected parameter space of this model. In this appendix we use Passarino-Veltman (PV) functions [28,48] for representing all analytic formulas of one-loop contributions to LFVHDs defined in Eq. (43). We also use notations for one-loop integral of PV functions, such as The analytic expressions for ∆ Donations with the participation of Y ± -boson + λ R,s * ai λ R,s bj m b m n j C 2 − λ L,s * ai λ L,s bj m a m n i C 1 + λ 0 ij λ R,s * ai λ L,s bj m n i m n j C 0 + λ R,s * ai λ R,s bj m n i m b (C 0 + C 2 ) + λ L,s * ai λ L,s bj m a m n j (C 0 − C 1 ) + λ L,s * ai λ R,s bj m a m b (C 0 − C 1 + C 2 ) , ∆ + λ L,s * ai λ L,s bj m b m n j C 2 − λ R,s * ai λ R,s bj m a m n i C 1 + λ 0 * ij λ L,s * ai λ R,s bj m n i m n j C 0 + λ L,s * ai λ L,s bj m n i m b (C 0 + C 2 ) + λ R,s * ai λ R,s bj m a m n j (C 0 − C 1 ) + λ R,s * ai λ L,s bj m a m b (C 0 − C 1 + C 2 ) , ∆ (7)Hs L = g 2 λ ± Hs f s 16π 2 m 2 W 9 i=1 −λ R,s * ai λ L,s bi m n i C 0 − λ L,s * ai λ L,s bi m a C 1 + λ R,s * ai λ R,s bi m b C 2 , ∆