Low energy phenomena of the lepton sector in an $A_4$ symmetry model with heavy inverse seesaw neutrinos

An extension of the two Higgs doublet model including inverse seesaw neutrinos and neutral Higgs bosons was constructed based on the $A_4$ symmetry in order to explain the recent neutrino oscillation data. This model can distinguish two well-known normal and inverted order schemes of neutrino data once both the effective masses $m_{\beta}$ in tritium beta decays and $\langle m\rangle$ in the neutrinoless double beta decay are observed. The lepton flavor violating decays of the charged leptons $e_b\rightarrow e_a\gamma$, $\mu\rightarrow3e$, the Standard model-like Higgs boson decays $h\rightarrow e_be_a$, and the $\mu$-e conversions in some nuclei are generated from loop corrections. The experimental data of the branching ratio Br$(\mu\rightarrow e\gamma, 3e)$ predict that the upper bounds of Br$(\tau \rightarrow \mu\gamma,e\gamma)$ and Br$(h\rightarrow e_{a}e_b)$ are much smaller than the planned experimental sensitivities. In contrast, the $\mu$-e conversions are the promising signals for experiments.


I. INTRODUCTION
Observation of neutrino oscillation requires that models beyond the Standard Model (BSM) must be considered to explain both properties of the tiny active neutrino masses and the structure of the lepton mixing matrix U PMNS named as the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix. The simple tri-bi maximal (TB) form of U PMNS was introduced in Refs. [1][2][3][4], which can be explained theoretically as the consequence of discrete symmetries such as A 4 group [5][6][7][8]. The TB form implies exact zero value of a mixing angle θ 13 defined in the standard form of U PMNS [9] which is inconsistent with non-zero but small θ 13 pointed out by experiments so that this form must be modified, see a recent review in Ref. [10]. Various modifications were carried out in order to looking for models as simple as possible . Many of them are A 4 models generating active neutrino masses based on the wellknown standard seesaw (SS) mechanism [20][21][22][23], some of them based on the inverse seesaw (ISS) [25].
In the A 4 models containing heavy neutrinos to generate active neutrinos through the SS or ISS mechanisms, the lepton flavor violating (LFV) couplings of neutrinos will give loop corrections to many LFV processes such as the decays of charged leptons (cLFV) e b → e a γ, µ → 3e, decays of the Standard Model-like (SM-like) Higgs boson (LFVHD) h → e + a e − b , and the µ-e conversions in nuclei. For the standard SS models, these corrections to the LFVHD are very suppressed [33], therefore only the branching ratio (Br) of the decays µ → eγ may reach experimental sensitivities. In contrast, using the so-called Casas-Ibarra parametrization [34] and the simple diagonal form of heavy neutrino mass matrix to determine the mixing parameters and neutrino masses, many ISS models predict large LFV corrections from heavy neutrinos to LFVHD [35][36][37][38]. Namely, Br(h → τ µ, τ e) can reach the order of O(10 −5 ) under the very small experimental constraint Br(µ → eγ) < O(10 −13 ).
In general, these sensitivities are still larger than the upper bounds predicted by the models containing only loop contributions to the LFVHD.
The most strict constraints from experiments for cLFV decays are Br(µ → eγ) < 4.2 × 10 −13 and Br(µ + → e + e + e − ) < 10 −12 [46], with the respective planned sensitivities will be 6 × 10 −14 and O(10 −16 ) [47]. These two constraints often result in much smaller values of Br(τ → µγ, eγ) than the recent experimental upper bounds, or even the planned sensitivities of O(10 −9 ). These cLFV decays and the µ-e conversions in nuclei were also investigated in the standard SS models [48,49] and ISS model [50], including the minimal supersymmetric (MSSM) versions [51,52]. Depending on the specific structures of the Higgs and lepton sectors of different discrete symmetric models, the allowed regions satisfying all cLFV experimental bounds will predict different possibility to observe the LFVHD in the future experiments.
Based on the SS models containing two SU (2) L Higgs doublets transforming as A 4 singlets [11,26], in this work we will introduce a non-supersymmetric A 4 × Z 3 × Z 11 model with the ISS mechanism (A 4 ISS) to generate active neutrino masses and U PMNS enough to explain the recent oscillation data and allow low mass scale of heavy neutrinos. In addition, the LFV signals will be discussed as other promoting channels to constrain the parameter space. Our work also pay attention to the LFVHD, which is often ignored in discrete symmetric models because the SM-like Higgs boson is difficult to realize in complicated Higgs potentials of both SUSY and non-SUSY versions. To keep the Yukawa term unchanged in the original SUSY versions, the discrete symmetries Z 3 and Z 11 are introduced to exclude unwanted terms appearing in the non-SUSY version. The model also consists of additional SU (2) L neutral Higgs singlets (flavon) enough to generate the active neutrino mass matrix corresponding to U PMNS close to the TB form in a special degenerate limit of the two independent parameters in the heavy neutrino mass matrix. Using a deviate parameter˜ to relaxing this limit will result in the real form of U PMNS . The observable parameters defined by the standard form of U PMNS will be determined based on the recent work [10]. Using this in the numerical investigation, we collected allowed regions of the parameter space to study low energy observable quantities such as effective neutrino masses related with the neutrinoless double beta (0νββ) and Tritium beta decay, the cLFV decays e b → e a γ, µ → 3e, and h → e a e b .
We note that the decay h → e a e b were rarely discussed in previously in discrete symmetric models, including the SUSY versions. The numerical results of the allowed regions will be compared with previous works as well as the current and future experimental constraints. In contrast with original SS models [11,26], in this work all allowed values of the Dirac phase δ are considered and the parameter˜ is solved exactly in this A 4 ISS model. Hence the allowed regions of parameters will be determined more exhausted, leading different predictions for the low energy observable quantities. In addition, different from the ISS models with the Casas-Ibarra parameterisation [34], the masses and the particular form of the total lepton mixing matrix of all neutrinos originated from the A 4 × Z 3 × Z 11 symmetric breaking will lead to new predictions for the LFV signals.
The A 4 ISS contains two SU (2) L Higgs doublets and only flavons, therefore it inherits many properties of the well-known two Higgs doublet models (2HDM) type I and II, see a review in Ref. [53]. Based on many discussions on the 2HDMs, we will constrain many important parameters affecting strongly on the signals of LFV processes predicted by the A 4 ISS model. Namely, the most important parameters are the ratio of the two vacuum expectation values (vev) of the two neutral components in the Higgs doublets t β , the charged Higgs mass, and the parameter s δ defining the deviation of the SM-like Higgs boson couplings between the A 4 ISS and the SM. Based on the experimental results from LHC searches and precision electroweak test, the recent constraints on the parameter spaces of the 2HDM related with the A 4 ISS model were discussed in detailed in Refs. [54][55][56]. In the future project of energy collision of 100 TeV, promising signals of heavy higgs bosons with masses at O(10) TeV were mentioned [57]. In direct signal of heavy Higgs bosons predicted by the 2HDM may also appear in the future e + e − colliders [58], where large allowed t β corresponds to the alignment limit s δ → 0. The predictions on cLFV signals may depend strongly on t β , leading to another channel to determine which SU (2) L Higgs doublet generates quark masses, i.e. the information to distinguish the type I and II of the 2HDM. This work is organized as follows. In section II, we introduce the A 4 ISS mechanism, constructing the analytic formulas of active neutrino masses and mixing parameters as functions of free parameters. We also give out the allowed regions of parameter space satisfying the recent neutrino oscillation data. The predictions of the effective neutrino masses corresponding to the two decays, neutrinoless double beta and Tritium beta are also discussed. In section III, important properties of the SM-like Higgs and charged Higgs bosons are summarized. Analytic formulas and numerical results relating to the LFV processes are separated into two sections IV and V. Finally, the summary of our new results is given in section VI.
There are four appendices presenting more details on the product rules of the A 4 symmetry, the full Higgs potential, the one-loop formulas contributing to the LFV decay amplitudes.

A. The particle content and lepton masses
The non-Abelian A 4 is a group of even permutations of 4 objects and has 4!/2 = 12 elements. The model has three one-dimension (1, 1 , 1 ) and one three-dimensional (3) irreducible representations. The important properties of this group and its representations needed for model construction were reviewed in the appendix A. The transformations for leptons and scalars under the total symmetry SU (2) L × U (1) Y × A 4 × Z 3 × Z 11 × U (1) L as well as their VEVs of the A 4 ISS model is shown in Table I. Here L is the normal lepton  Lepton where the electric charge operator is well-known as Q = T 3 + Y /2.
Considering the effective operators up to five dimension (dim.) needed to generate masses of lepton, the Yukawa Lagrangian respecting the total symmetry consists of two parts. In particular, the first part is renormalizabe as follows While the second part consisting of all effective operators of five dim., including all terms breaking the lepton number L relevant with neutrino masses, is Here Λ is the cut-off scale of the model under consideration,h u = iσ 2 h * u , the charge conjugation of the neutral leptons N R and X R are N c L ≡ P L N c = (N R ) c = C(N R ) T and X c L = (X R ) c , respectively. After spontaneous symmetry breaking, the charged lepton mass matrix comes out diagonal with m e = yev T v d Λ , m µ = yµv T v d Λ , and m τ = yτ v T v d Λ . Correspondingly, the Higgs doublet h d plays a similar role to the SM Higgs doublet in generating charged lepton masses. Only the last term in Eq. (3) breaks the lepton number L with two units, giving neutrino mass term µ X ≡ (λ X v u v d /Λ), which can be small so that the ISS mechanism can work. The non-renormalizable terms generating charged lepton masses can be seen as the tree level mass originated the exchange of the heavy vector-like leptons transforming as E L,R = (E 1 , E 2 , E 3 ) T L,R ∼ (1, −2) (3,1,1,1) with the total symmetry (SU (2) L , U (1) Y ) (A 4 ,Z 3 ,Z 11 ,L) .
The renormalizable Lagrangian relating with E L,R that respect the total symmetry is where we assume that Λ v S , v T , u, u so that all of the above Higgs bosons do not contribute significantly to very heavy vector-like lepton masses E 1,2,3 . Therefore, m Ea ≡ Λ for all a = 1, 2, 3. The non-renormalizable terms are reduced form the diagrams given in Fig.1. For example in the regions where the vector-like lepton momenta are much smaller than their masses we the effective term i ye The Lagrangian for neutrino mass is where m D , M R , and µ X are 3 × 3 Dirac and Majorana neutrino mass matrices have the following forms and The tilde notation inx implies a complex parameter, and x ≡ |x| > 0. Without loss of generality, the real and positive M 0 is assumed in Eq. (8) and 0 ≤ φ 1 < 2π. The effective neutrino mass matrix m ν is then obtained by the ISS relations [59], which is a specific of the general SS framework [60,61]: where m ν ≡ U * 3νm ν U † 3ν ,m ν ≡ diag (m n 1 , m n 2 , m n 3 ) relating with the mixing matrix U 3ν and the masses m na (a = 1, 2, 3) of the three active neutrinos, leading to the following form of the U PMNS , where U L,e is defined by the relation U † L,e m D U R,e = diag(m e , m µ , m τ ). The standard form of the U PMNS is the unitary matrix defined as follows [9] U where s ij ≡ sin θ ij , c ij ≡ cos θ ij , i, j = 1, 2, 3 (i < j), 0 < θ ij < 180 [Deg.] and 0 < δ ≤ 720 [Deg.]. In this work, U L,e = U R,e = I 3 , hence where the phases ψ 1,2,3 are absorbed into the charged lepton states. The standard form U PDG PMNS provides the experimental quantities s ij and δ, α 21 and α 31 . Let us remind that the TB framework of the lepton mixing matrix was predicted in an A 4 model with the standard SS mechanism. It also happens in this model in the degenerate condition thatρ =κ, corresponding to˜ = 0. As a result, the neutrino mixing matrix U 3ν given by Eq. (9) have the TB form: corresponding to s 13 = 0. This is in contrast with the experimental data of neutrino, which is divided into two cases of normal (NO) and inverted (IO) schemes. We consider the real case, where all θ ij is in the 3σ ranges of the experimental data. We assume a solution that the deviation from the TB data arises from only the condition that ρ =κ in the matrix M R given in Eq. (7), equivalently˜ = 0. Then, the mixing matrix U 3ν in Eq. (10) is calculated by writing it as which can be identified with the well-known form given in Eq. (12). We then have, where U TB is given in Eq. (14) and The form of m ν results in the form of U 1 as follows: where c θ ≡ cos θ, s θ ≡ sin θ, and φ 0 are real. The matrix U 1 is found by diagonalizing the matrix m † ν m ν , leading to the total the neutrino mixing matrix given in Eq. (21). It can be identified with the standard form given in Eq. (12) by the following relation [10,62] , cos δ = cos 2θ 13 cos 2θ 23 They were found by the requirement that | (U 0 PMNS ) ij | = | (U TB U 1 ) ij | for all i, j = 1, 2, 3. These equations show that apart from θ, φ 0 , other parameters s 12 and the Dirac phase can be written in terms of s 13 and s 23 . The formulas of cos φ 0 and cos θ in (22) result in the following well-known relations [10]: sin φ 0 = − sin 2θ 23 sin δ, and Therefore, using tan 2 φ 0 + 1 = 1/ cos 2 φ 0 and cos φ 0 / cos δ < 0 to formulate cos φ 0 as follows Based on Eq. (22), s 2 23 is written in term of the following function of δ and s 13 , From now on, s 23 and φ 0 are investigated as the above functions of s 13 and δ. The requirements | cos φ 0 |, | sin φ 0 | ≤ 1 always satisfy under the recent 3σ data allowing only very small s 2 13 . In contrast, 3σ allowed range of s 2 23 gives an upper bound on δ, more strict than that from the 3σ experimental data as given in Fig. 2. The new upper bound of δ for both NO and IO schemes is Apart from relations in Eq. (22), the matrix U TB U 1 casts into the well-known standard form: where the unphysical phases ψ i are absorbed into the charged lepton states. In contrast, the two phases γ 2,3 will contribute to the Majorana phases α ij . By identifying ij , ψ i and γ 2,3 are found as follows: by identifying that e ix = e iy is equivalent with x = y without the unnecessary repeated term k2π, k ∈ Z. Hence, the non-zero contribution to the Majorana phase is consistent with the result mentioned in Ref. [10].
From now on, three parameters θ 12 , θ 23 , and φ 0 will be written as functions of s 13 and δ based on the relations given in Eqs. (22), (23), (24), and (25). Taking these functions for diagonalisation m ν by requiring that [U T 1 m ν U 1 ] 13 = 0 will lead to s 2θ (m ν ) 33 13 , and (m ν ) 33 are given in Eq. (19). The result is leading to the following exact solution of˜ , Here, we choose˜ guarantees that lim t 2θ →0˜ = 0, consistent with our assumption that θ = 0 results in the TB form of U PMNS . We note that formula of˜ is more general than that given in the original SS models [11,26], where˜ was assumed to be small for the approximate solutions. The Majorana phases defined in Eq. (17) are also formulated based on the phases of U T 1 m ν U 1 . Then, the more precise form of U PMNS defined in Eq. (12) is where γ 3 is given in Eq. (28) contributes to the phase α 31 , apart from those come from the phases of the neutrino masses in the diagonal matrix U T 1 m ν U 1 . At this step, all of the active neutrino masses and mixing parameters can be formulated in terms of the five independent parameters s 13 , δ, m 0 φ 1 , and κ. Until now, we do not now the orders of m 0 and κ, hence it is necessary to take some numerical estimation to figure out these orders.
As a result, ∆m 2 32 is written as a function of ∆m 2 21 andκ = κe iφ 1 for both NO and IO schemes. The experimental data of s ij , ∆m 2 21 , ∆m 2 32 , and δ will give constraints on κ and φ 1 . In the TB limit, the right-hand sides of Eqs. in (33) depend only on κ and φ 1 . The condition ∆m 2 21 > 0 requires κ < 1.6 being useful for estimating the deviation of (ρ −κ) to looking for a real mixing angle s 13 .
For the first estimation of κ and φ 1 with˜ = 0, we find numerically that ∆m 2 21 /m 2 0 depend weakly on s 2 13 . Therefore, it is fixed at the best-fit point when ∆m 2 21 /m 2 0 is plotted as a function of κ and φ 1 in the range 0 ≤ φ 1  Adding a requirement that the formula of ∆m 2 21 /m 2 0 given in Eq. (33) must be positive, values of κ > 1.9 are excluded, see the numerical illustrations in Fig. 3. The allowed regions of {κ, φ 1 } depend rather complicatedly on δ. In Fig. 3, the contour plot of ∆m 2 32 /m 2 0 are also shown, where the two yellow and non-color regions supporting the respective IO and NO schemes are separated by the constant curve ∆m 2 32 /m 2 0 = 0. It can be estimated that, the allowed regions for the NO scheme always require that 1.4 × 180/π < φ 1 < 4.
Now, from numerical illustrations shown in Fig. 4, the allowed regions (blue) are more constrained. We also add here the contour plots of m 0 as a function of φ 1 and κ given in Eq. (33) and ∆m 2 21 is fixed at the best-fit point. The numerical investigation above shows that for the NO scheme, κ must satisfies 0.4 < κ < 1.6. And the order of m 0 is around 10 −10 eV. In addition, every allowed φ 1 corresponds to a very narrow allowed range of κ. ]. In addition, the allowed ranges of κ and m 0 are similar to those from the NO scheme, hence will be used to scan for determining more general allowed regions in the following numerical investigation.

B. Low energy observables
As we discussed above, the neutrino mass matrix in Eq. (19) depends on the set of five free real parameters (s 13 , δ, φ 1 , κ, m 0 ) used to scan numerically. By fixing s 13 , δ, and ∆m 2 21 , we have estimated the reasonable ranges of all parameters m 0 , κ, and φ 1 . They will be   used to scan all five independent parameters (s 13 , δ, m 0 , κ, φ 1 ) to collect all allowed points satisfying the 3σ experimental data in the general case. For both NO and IO schemes, the unknown parameters get random values in the following ranges: and the two parameters s 13 and δ runs over 3σ allowed range of the experimental data. We stress that wider ranges of m 0 and κ were investigated before choosing the best ranges for the following numerical illustration. The parameter spaces (κ, , φ 0 , φ 1 , φ 2 = arg[˜ ]) and their correlations are respectively plotted in Fig. 5, where the red and the blue patterns represent the allowed regions predicted by the NO and the IO schemes, respectively. Hereafter, we continue using these conventions if there is no more explanation. We note that the upper bound values of m 0 for the typical NO (IO) scheme was estimated from our numerical investigation. We can find that the allowed regions of the model's parameters are separated completely for the two schemes. In particularly, the ranges 0.35 < κ < 1.5, 0.1 < < 2.3, ] support the NO scheme while the IO implies the allowed values of the parameters as 0.3 < κ < 1.3, 0 < < 3,  (15) and (16).
The light neutrino masses predicted by the model are respectively plotted in Fig. 6 as function of the light neutrino mass scale m 0 for the NO (left panel) and the IO (right panel).
There the red, blue, green plots represent for m 1 , m 2 , m 3 , respectively. We can recognize that the neutrino masses are strong hierarchy with small values of m 0 and they can be quasi  We all know that neutrino oscillation experiments could not determine the absolute scale of active neutrino masses. Instead, it can be measured by non-oscillation neutrino experiment such as Tritium beta decay [63], neutrinoless double beta decay (0νββ) [64], or by cosmological and astrophysical observations [65]. As the model consequences, we would like to study the effective neutrino masses associate with 0νββ (| m |) and beta decay (m β ) which are defined as | m | = In Eq. (36), we can use directly U 3ν given in Eq. (17) instead of U PMNS defined in Eq. (31).
The predictions of the allowed regions presenting the relation between effective masses | m | and m β with the lightest active neutrino mass, m lt , are respectively plotted in Fig. 7 and Fig. 8 (left panel). Whereas, the correlation between m β and | m | is shown in the right panel of the Fig. 8.   The consequence is that the IO scheme is excluded for m β < 0.05 eV, or | m | < 0.015 eV. In contrast, the NO is excluded for m β < 0.01 eV or | m | < 0.005 eV. In addition, the sharped allowed regions in the right panel of Fig. 8 implies the very strict relations between m β and | m | for both schemes. On the other word, the model is very predictive. Because once both of these quantities are observed, the reality of the model is immediately confirmed.
If it survive, based on the two separated regions corresponding the NO and IO schemes, the model can figure out only one of the two schemes survive.

III. THE HIGGS SECTOR
The same as the SM, the covariant derivative for local SU (2) L ⊗ U (1) Y symmetry is Here T a (a = 1, 2, 3) are the generators of the SU (2) L symmetry. The kinetic terms of the two Higgs doublets generating gauge boson masses are: The SM charged gauge bosons are W ± ≡ (W 1 ∓ iW 2 )/ √ 2. The same relation known in two Higgs doublet models (2HDMs) for the two vevs, where s 2 W = 0.231. The model also contains the SM Z boson satisfying W 3 A detailed analytical calculation on the mass and mixing of the Higgs in the 2HDMs was done previously, for example in Ref. [69]. Similarly to the very important parameter t β defined by the ratio of the two neutral Higgs VEVs in the MSSM and the 2HDM, t β in the model under consideration is defined as follows: It is noted that from now on we will use the following conventions for any angles x = δ, α, β, 2α, 2β, and (β − α): t x ≡ tan(x), s x ≡ sin(x), c x ≡ cos(x). Because the Yukawa couplings of charged leptons with h d in Eq. (3) is fixed, t β defined in this model is consistent with Refs. [51,70], but equivalent to 1/t β defined in other MSSM and 2HDMs discussed previously.
For the Higgs spectrum, this model must contain at least one SM-like Higgs boson observed by the LHC. It is one of those included in the squared mass matrix of the CP-even Higgs boson, which is a 10×10 matrix consisting of a large number of the independent Higgs self-couplings. In this work, we will choose a simple case of the realistic Higgs potential, which is summarized in the appendix B. Here we will focus on the identification of SM-like Higgs boson. In particular, a simple regime is chosen that these Higgs doublets decouple to other Higgs singlets, and a soft term generate masses for CP-odd neutral Higgs and results in heavy charged Higgs bosons, the same as those appear in the well-known 2HDM [53].
The model contains two singly charged Higgs bosons ϕ ± and two massless states G ± W which are Goldstone bosons eaten by gauge bosons W ± . The masses and relations between the original and physical states of the charged Higg components are: Regarding CP-odd neutral Higgs components, the squared mass matrix has a zero determi- s β c β = m 2 ϕ − λ 4 v 2 , implying that µ 2 12 > 0. For the CP-even neutral Higgs corresponding to the simple case we assume above, the total squared mass matrix will separate into two submatrices, namely a 2 × 2 and a 8 × 8 ones. The 8 × 8 matrix gives eight physical heavy Higgs bosons with masses depending on heavy VEVs v S and v T , see appendix B for the details. In the original basis (S u , S d ) T given in Eq. (1), the 2 × 2 matrix containing the SM-like Higgs boson has the following form, where α is defined based on the following relation Because the SM-like Higgs boson h were found experimentally at LHC, we require that the its couplings with normal fermions and gauge bosons have small deviations with those predicted by the SM, see the Feynman rules given in table II. Here all quarks multiplets are assumed to be singlets under all A 4 and other discrete symmetries listed in table I. The Yukawa Lagrangian of quark is Correspondingly, the top quark mass is m t = v d y t , leading to the perturbative limit m t /v d = y t < √ 4π. As a result, we have The model now treats like the 2HDM type-I, which the recent constraints on the model parameters were given in Ref. [55,57]. The typical constraints on the free parameters are |s δ | ≤ 0.05 and t β ≤ 3.3.
There is another assignment of quark sector corresponds to the case of the 2HDM type-II and the MSSM: all upper quarks only couple with h u , leading to the condition t β ≥ 0.3 for the perturbative limit. The recent constraints of the model parameters were given in Refs. [54,57], where the typical ranges of free parameters are |s δ | < 0.008 and 0.2 ≤ t β ≤ 5 after casting into the model under consideration. The Lagrangian contains Yukawa couplings of quarks with charged Higgs boson is: where V is the Caibibbo-Cobayashi-Maskawa matrix, m qa with q = u, d and i = 1, 2, 3 is the mass of the quark q a .
We note that t β defined in this work is equivalent with 1/t β defined in Ref. [54,55] for both type-I and II of the 2HDM. Our numerical calculation does not depend on the assignments of quark, except the µ − e coversion rate, which was discussed for the MSSM, including the contributions from the charged Higgs bosons in the 2HDM type II. But the phenomenology investigated in this work strongly depend on t β , leading to another channel to determine which SU (2) L Higgs doublet generates quark masses.
From the recent experimental data, the coupling of the SM-like boson with normal charged leptons, namely he a e a given in Table II, must be consistent with the coupling predicted from in the SM, leading to a consequence that s α −c β . A small deviation from the SM couplings correspond to a small parameter δ satisfying equivalently δ ≡ π 2 + α − β 0, which is known in the 2HDMs. The constraints of δ, t β , and charged Higgs bosons comes from the electroweak and Higgs precision measurements discussed recently for the 2HDMs [54,55,58], In conclusion, in the numerical investigation, the free parameters in the Higgs sector we will use are β, δ, m h ≡ m h , m ϕ , λ 2 , and λ 4 . The SM-like Higgs mass m h ≡ m h 125.09 GeV was determined experimentally from LHC. Three parameters α and λ 1,3 are determined from three Eqs. (48), (44), and (45). We are interesting in the regions allowing large range of t β , hence we will fix λ 4 = 0, and m A = m H 2 = m ϕ , which were shown from previous discussions [55][56][57]. Based on the general formulas of the three parameters λ 1,2,3 given in Ref. [69], they can be written as follows: The analytic formulas of these three Higgs self couplings will be used to calculate the coupling hϕ + ϕ − that contributes to the LFVHD.  2) and (3). The corresponding mass matrix is diagonalized by a total 9 × 9 unitary mixing matrix U ν determined as follows: where n = (n 1 , n 2 , ..., n 9 ) T and n i = (n iL , n iR ) T (i = 1, 2, .., 9). The matrix U ν is parameterized as follows [34,71,72], where O is the 3 × 3 matrix with all elements being zeros, V 6 is the unitary matrix, exp(x) ≡ ∞ i=0 x i /i!, and R is the 3 × 6 matrix satisfying the ISS condition that max|R ai | 1 for all a = 1, 2, 3 and i = 1, 2, .., 9. Accordingly, the ISS relations given in Eq. (9) are obtained from the expansion U ν up to the order O(R 2 ). Correspondingly, the sub-matrix R and heavy neutrino masses in this case are where .., M 6 ) presents six heavy neutrino masses.
Before determining Feynman rules to calculate Br of LFV decays, we remark some requirements to guarantee the ISS relations given in Eqs. (9) and (53). In principle, in order to consistent between the approximation solution from the ISS relations and the exact one obtained directly from calculating numerically the total neutrino mass matrix given in Eq. (5), the condition max(|R ai |) 1 must satisfy, equivalently max |m D (M T R ) −1 | ai 1, leading to a consequence that the two scales v u f and M 0 given in (7) must satisfies v u f /M 0 1.
Direct numerical calculation to derive masses and mixing of the active neutrino from the total neutrino mass M ν , it is necessary to know the input as the following set of free parameters {m d , M 0 , µ X ,κ,ρ}. We limit our investigation in the allowed regions of the parameter space indicated in section II B. Namely, two parametersκ andρ reduce to two allowed values of κ and φ 1 . The parameter µ X is considered as a function of m 0 from Eq. (20), then allowed values of m 0 can be determined from the numerical investigation.
These allowed values will be used as the inputs. The remaining parameters m D ≡ vf , and M 0 do not appear in the approximate solution of active neutrino data using ISS relation (10) because it is absorbed in m 0 . Instead, the constraint of m 0 is used to determine µ X . In contrast, determining M ν requires all of the three parameters. In numerical investigation, the exact solution is found by using the software mathematica 11 to find the mixing matrix U ν relating with the eigenvectors of the matrix M ν † M ν and the neutrino masses.

B. Feynman rules for LFV decays
In the Yukawa Lagrangian parts (2) and (3), couplings relating to LFV decays are where the third term results in the coupling hn i n j for the ISS mechanism discussed in detail in Ref. [35,73]. We can prove that which was introduced firstly in Ref. [35]. Using the definition λ ij ≡ D ij m n i + D * ij m n j needed to write the right Feynman rules for Majorana neutral lepton in terms of Dirac spinors [33,36,72], the coupling hn i n j is written in the symmetric form as follows: In the unitary gauge, the Feynman rules for vertex couplings relating with the decay processes in this work are given in table III, consistent with those mentioned in the 2HDMs [53].
In the limit m 2 a /m 2 b 1, m a,b being masses of charged leptons e, µ and τ , the Brs of the cLFV decays e b → e a γ is determined as follows [74,75],  [26], and consistent with [74,75]. The two loop contributions mentioned in the 2HDM type-X model given in Ref. [76] do not appear in our model.
The analytic formulas base on the analytical results for non-supersymmetric contributions to LFV decays in the 2HDM presented in Refs. [48,51,52], and results Refs. [77][78][79] which general formulas can be used for the 2HDM. The µ-e conversion rate for the ISS model were given in Ref. [50]. The partial decay width µ → 3e is taken from Ref. [51], namely only the non-suppersymmetric contributions are collected. The decay rate is where Γ µ = α 2 W m 5 µ /(384πm 4 W ), α W = g 2 /(4π), and the loop functions are listed in the appendix C. In the limit of the ISS model, the Eq. (56) is consistent with that given in Refs. [48,49]. Apart from the gauge boson and Higgs contributions taken from Refs. [48,49], the Higgs contributions were checked with the results given in Refs. [80,81].

E. µ − e conversion in nuclei
Based on the results given in Refs. [48,51,70], we collect all one-loop non-supersymmetric contributions to the µ − e conversion rates, see the detailed formulas listed in appendix C.
In the model under consideration the µ − e conversion rate R J µ→e in a nuclei J consisting of Z protons and N neutrons is where where Q q and I 3 are the electric charge and iso spin of the quark q = u a , d a . To determine V X q with X = L, R, we just consider the limit of the 2HDM type II so that the couplings of the charged Higgs bosons with all quarks in our model are the same as those given in Ref. [51]. Well-known values of Z ef f , F p , and Γ capt corresponding to various nuclei are given in table IV [48,82]. The form factors F X γ,µe , F X Z,µe , G X γ,µe , and F µeqq box are collected in   from SM is small with |s δ | ≤ 0.05 defined in Eq. (48), hence we can ignore this change in our numerical investigation. In notations constructed in [83], the ∆ (ba)L,R can be written as where ∆ In this section we will apply the allowed regions of the parameter space mentioned above to investigate the LFV decays. First, we discuss the independent parameters needed to determine numerically the masses and mixing matrix of all neutrinos from the total mass matrix given in Eq (5). The five independent parameters s 13 , δ, φ 1 , κ and m 0 will be constrained from the experimental data, as we have presented precisely. In exact numerical calculation using the direct neutrino mass matrix (5), more unknown parameters are m D ≡ vf , t β , and M 0 .
Apart from the free parameters mentioned above, the are unknown parameter relating with the Higgs sector. In particularly, contributions of charged Higgs bosons ϕ ± depends on m ϕ , the mixing angle δ between the CP-even neutral Higgs bosons, and four Higgs-self couplings λ 1,2,3,4 in coupling factor hϕ + ϕ − . But all of them are not independent and we can choose the new independent parameters are m h , m ϕ , and δ as we discussed previously. The perturbative constraints for these Higgs self-couplings of the model are [54,69,84]: We fix m h = 125.09 GeV. The Dirac mass scale m D = vf < 174 × √ 4π 616 GeV.
The heavy neutrino masses are originated from the A 4 breaking scale M 0 hence they can be very large. This situation is completely different from that discussed previously in Refs. [36], where heavy neutrino masses are bounded from above because of the Casas-Ibarra parametrization [34] specializing the structure of Dirac mass term m D , resulting in the perturbative limit of the Yukawa coupling. In the numerical investigation, we require M 0 /(m D s β ) ≥ 10 1, necessary to obtain the consistent ISS relations given in Eqs. (9) and (53).
On the other hand, in the case of λ 4 = 0 and s δ → 0, the values of t β is relaxed to a large values of t β ∼ 50. In the study on the parameter space corresponding more heavy masses of the Higgs bosons [57], the large t β values are still allowed. For the case of the 2HDM type-I, based on the recent results given in Ref. [55], the allowed regions of parameters are: In the numerical investigation for studying the LFV phenomenology using the allowed values of the set {s 13 , δ, m 0 , κ, φ 1 } that obtained from scanning the ranges given in Eq.
To looking for regions of the parameter space allows large Br(h 0 1 → τ µ, τ e) and satisfy the constraints Br(µ → eγ) < 4.2 × 10 −13 and Br(µ → 3e) < 10 −12 , we need to scan s δ in the range |s δ | ≤ 0.05 (0.008) for model type-I (II) and adding a requirement given in Eq. (61). We note that the small ranges of m D GeV were scanned but they give suppressed values of Br(h 0 1 → e a e b ). In the numerical investigation, we just collect points satisfying max[Br(h 0 1 → µτ, eτ ) ≥ 10 −11 . To discuss on the µ-e conversion rates predicted by the 2HDM type-II only the following allowed range of t β is considered: 0.3 ≤ t β ≤ 50. The most interesting allowed range 0.2 ≤ t β ≤ 5 given in Ref. [54] will also be remarked.

LFVHD and cLFV processes
The correlation between Br(h 0 1 → e a e b ) and Br(µ → eγ) are plotted in the consideration, the A 4 symmetry results in a very strict structure of the total neutrino mass matrix, which is the origin of the very strict relations between the cLFV and LFVHD decays.
We continue with the numerical results for the cLFV decays Br(e b → X) with {e b , X} = Regarding to the µ−e conversion corresponding to the 2HDM type-II, with 0.3 ≤ t β ≤ 50.
The correlations between Br(µ → eγ) and the four µ-e conversions rates are shown precisely in Fig. 13. All µ-e conversion rates are constrained strictly by the data of the decay µ → γ. They decrease with smaller Br(µ → γ). The allowed regions corresponding to the three nuclei Ti, Au, and Pb are nearly the same, while the allowed region for Al is more narrow. Finally, we remind that the above allowed regions of parameters satisfy the recent experimental bound (µ → 3e) < 10 −12 . If this channel is not observed with planned sensitivity of 10 −16 , the upper bounds of cLFV decays and or µ-e conversion rates will be more suppressed, see illustration in Fig. 15   The reality of the model and many interesting predictions on the observable quantities such as m β , | m |, cLFV decay rates, and µ-e conversion rates in nuclei will be confirmed or ruled out by the upcoming experiments.
Only rep. 3 is used for generating the neutrino mass matrix. In the mentioned basis, T is complex and T * = T in general so the complex conjugate representation r * of a representation r (r = 1 , 1 , 3) is not the same as r, although they are all real reps.. It is determined by the following rules [94,95]: with the Higgs self-couplings. Because of the very large number of the Higgs self-couplings, the VEV pattern assumed in this work is easily guaranteed. Therefore, the lengthy and unnecessary minimal conditions will not be presented here.
In general, the squared mass matrix of the CP-even Higgs bosons are 10 × 10 matrix, where the main contribution to the SM-like Higgs boson arises from the two Higgs doublets h u and h d . Hence, for simplicity in studying the LFV decay of the SM-lik Higgs boson, we will choose the regime that these Higgs doublets decouple to other Higgs singlets, namely To generate non-zero masses for CP-odd neutral and charged Higgs bosons we adopt a soft term breaking Z 3 × Z 11 in the Higg potential, namely µ 2 ud (iσ 2 h u ) T h d + h.c.. included in the second line of the Higgs potential (B1). It results in that the masses and eigenstates of all Higgs bosons arising from two Higgs doublet h u and h u are the same as those well-known in the 2HDM. The one-loop three-point PV functions [96], called C− functions, which specific definitions were given in Ref. [26]. For cLFV decay processes e b → e a γ, where are the masses of charged leptons m a,b satisfy m 2 a,b /m 2 W 1 and m 2 a,b /m 2 ϕ 1, the C-functions are where t = M 2 1 /M 2 2 . The value t = 1 gives C 0 = −1/(2M 2 2 ), C 1 = 1/(6M 2 2 ), and C 11 = −1/(12M 2 2 ). Contributions from W and ϕ ± bosons to the Br(e b → e a γ) defined as D (ba)L,R given in Eq. (55) are calculated based on the general form given in Ref. [75], where C (ba)L,R = U ν * bi U ν ai 2(C 12 + C 22 − C 1 )m 2 W + m 2 b (C 11 + C 12 + C 1 ) +m 2 n i (C 0 + C 1 + 2C 2 + C 12 + C 22 ) , U ν * bi U ν ai 2(C 11 + C 12 − C 2 )m 2 W + m 2 a (C 12 + C 22 + C 2 ) +m 2 n i (C 0 + 2C 1 + C 2 + C 11 + C 12 ) with C 0,a,ab = C 0,a,ab (m n i , m W , m W ), and with C 0,a,ab = C 0,a,ab (m n i , m ϕ , m ϕ ). In the limit m 2 = 0, C W (ba)L,R is consistent with that given in [97][98][99]. Also, with t ϕ,i = m 2 n i m 2 ϕ , C ϕ (ba)L,R have consistent forms with Ref. [99]. Loop functions relating with only gauge bosons are [48,49]: x − y (4 + xy) The loop functions relating with both charged gauge and Higgs bosons are included in Refs. [51,52]. Here, we use the results given in Ref. [51] and the analytic functions given in Ref. [80] to cast the contributions to the Higgs and gauge bosons into the analytic functions summarized in the following.