Heavy singly charged Higgs bosons and inverse seesaw neutrinos as origins of large $(g-2)_{e,\mu}$ in two higgs doublet models

We show that simple extensions of two Higgs doublet models consisting of new heavy neutrinos and a singly charged Higgs boson singlet can successfully explain the experimental data on muon and electron anomalous magnetic moments thanks to large chirally-enhanced one-loop level contributions. These contributions arise from the large couplings of inverse seesaw neutrinos with singly charged Higgs bosons and right-handed charged leptons. The regions of parameter space satisfying the experimental data on $(g-2)_{e,\mu}$ anomalies allow heavy singly charged Higgs boson masses above the TeV scale, provided that heavy neutrino masses are above few hundred GeV, the non-unitary part of the active neutrino mixing matrix must be large enough, two singly charged Higgs bosons are non degenerate, and the mixing between singly charged Higgs bosons must be non-zero.


I. INTRODUCTION
The latest experimental measurement for the anomalous magnetic moment (AMM) of the muon a µ ≡ (g − 2) µ /2 has been reported from Fermilab [1] and is in agreement with the previous experimental result measured by Brookhaven National Laboratory (BNL) E82 [2].
There are many types of different versions of THDM, such as for instance: type I, type II, type X and type Y, which are discussed in Ref. [64], where collider phenomenology will requirements that can be experimentally tested in the near future [58,65]. The conclusion is also true for the Minimal Supersymmetric Model (MSSM) [66][67][68]. Another solution for MSSM requires large SUSY threshold corrections enough to change both sizes and signs of the SUSY electron and muon Yukawa couplings [69], hence allow regions satisfying both (g − 2) e,µ anomalies with large TeV scale of slepton masses, but rather large t β ≥ 70 is needed. As usual solution, there is a variety of extensions of THDM that successfully explain and accommodate the (g − 2) e,µ anomalies by adding new vector-like fermions as SU (2) L multiplets [38][39][40][41][42][43] as well as SU (2) L singlets of neutral and charged Higss bosons [55,57,[71][72][73][74]. Only few THDM models with scalar singlets and fermions can successfully explain the experimental data of both muon and electron anomalous magnetic moments [42,70,72]. Our work here will introduce a more general solution that a wide class of the THDM by adding heavy right handed (RH) neutrino singlets and one singly electrically charged Higgs boson can give sizeable one-loop contributions enough to explain the experimental ranges of values of both muon and electron anomalous magnetic moments in a wide region of parameter space.
This paper is organized as follows. In Sec. II, the THDM extension with new RH neutrinos and singly charged Higgs boson (THDMN R h ± ) will be introduced, where we pay attention to the leptons, gauge bosons, and Higgs sectors, giving all physical states as well the couplings that may give large one-loop contributions to AMM. In Sec. III, the simple inverse seesaw (ISS) and analytic formulas for one-loop contributions to AMM are constructed. Numerical discussion will be shown in details. Our main results will be collected in Sec. IV.

II. REVIEW OF THE THDMN R h ±
In this section we will focus on the analysis of the lepton sector of the THDMN R h ± .
For details of the quark sector of different types of the THDM, the reader is referred to Refs. [75,76]. The leptonic, quark and scalar spectrum with their assignments under the SU (2) L × U (1) Y gauge group are given by [75,76]: Note that we use the old convention for the defining the hypercharge via the original Gell-Mann-Nishijima relation: Here we add six new neutrino singlets to the leptonic spectrum of the model in order to implement the standard ISS mechanism. Furthermore, the model scalar sector is augmented by the inclusion of singly charged Higgs bosons, which allow to generate the new Yukawa interaction (N IR ) c e aR h + thus resulting in sizeable one-loop contributions to AMM.
One of the most important parameters introduced in the THDM is We will also follow the notations s x ≡ sin x, c x ≡ cos x, and t x = s x /c x for any parameter x.
Because the Yukawa couplings of up type quarks and Φ 2 are fixed, we always have a lower bound t β ≥ 0.4. In the THDM type I, all fermions couple with only Φ 2 , hence no upper bounds for t β arise. In the THDM type II and Y, all down type quarks and charged leptons couple with Φ 1 , hence t β ≤ 91. In the type-X, the Yukawa coupling of tau with Φ 1 gives t β ≤ 348.
The THDM type-II and type-X where lepton doublets couple with only Φ 1 are called THDM type-A. The corresponding lepton Yukawa terms are: where a, b = 1, 2, 3 are family indices; k = 1, 2; and I = 1, 2, 3, ...6 are the number of new neutral lepton singlets. A rough condition for the pertubative limit of Y h is |Y h Ia | < √ 4π, and should be smaller for trusty values [77]. The Z 2 charges of N IR and h ± are always chosen to guarantee the Z 2 symmetry of Lagrangian (6). Using the Z 2 charge assignments discussed in Ref. [64], two Higgs triplets Φ 1,2 have different Z 2 charges, therefore N IR couple with only one of them. In addition, we can consider the inclusion of new Z 3 discrete symmetry in the model in order to forbid the quartic scalar interaction Φ † 1 Φ 2 2 + H.c. that generates active neutrino masses at the one-loop level [105]. This guarantees that the Yukawa couplings Y h Ia do not have suppressed constraints arising from the neutrino oscillation experimental data. The detailed assignments of the Z 2 and Z 3 charges are shown in table I. The Z 3 TABLE I. Z 2 × Z 3 charge assignments of the model type-A for the case I (II) with heavy neutrinos couplings with Φ 1 (Φ 2 ), w = e i2π/3 Q aL u aR d aR L a e aR N 1,2,3R N 4,5,6R Φ 1 Φ 2 h + to the experimental value, M h < 126 GeV [106]. This problem can be solved in the model by adding new scalars such as in the THDM [107], and even with the models by adding neutrino gauge singlets [108][109][110]. The reason is that additional scalar couplings to Higgs doublets generating SM-like Higgs boson mass give positive one-loop contributions to the β-functions of the quartic couplings, hence the negative values of these couplings implying the vacuum instability will be relaxed at higher energy scale. This consequence still hold with the appearance of large Yukawa couplings of new fermions with Higgs doublets, which give negative one-loop contributions to the β-functions of the quartic couplings. On the other hand, the upper bounds of these Yukawa couplings are required besides their bounds from the perturbative and unitarity limits, previously discussed for THDM in refs. [111], for example. This situation also applies to the THDMN R h ± model under consideration.
On the other hand, the appearance of the singly charged Higgs boson h ± resulting on the Yukawa coupling matrix Y h , which is irrelevant with any Higgs doublets. Hence Y h is not constrained by the requirement of vacuum stability. Furthermore, the quartic couplings of h ± with Higgs doublets in the Higgs potential will result in the higher energy scale satisfying the vacuum stability.
In the THDM type-I and type-Y, where charged lepton couple with Φ 2 , the first term of Eq. (6) should replace Φ 1 with Φ 2 . We call them the THDM type-B and will discuss later that the qualitative results for AMM do not change. This model type including the model introduced in Ref. [42].
We also emphasize that the models under consideration including the Yukawa part discussed on Ref. [55] focusing on one particular case of Φ k = Φ 2 . As we will show that in our numerical analysis, the regions of the parameter space predicting ∆a e,µ consistent with the experimental data, favor small values of t β , including the range t β ∈ [0. 4,10] which was not mentioned in Ref. [55].
Given that we are working in the basis where the charged lepton mass matrix is diagonal, the leptonic mixing entirely arises from the neutrino sector. This implies that, in order to successfully reproduce the neutrino oscillation experimental data, the neutrino mixing matrix is parameterised in the following form [78]: where V is a 6 × 6 unitary matrix; R is a 3 × 6 matrix satisfying |R aI | < 1 for all a = 1, 2, 3, and I = 1, 2, ..., 6. The 3 × 3 unitary matrix U PMNS is the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix [79]. The condition |R aI | < 1 needed to guarantee that the U ν as a power series of R given in Ref. [78] is convergent and the approximation up to the order R 2 results in seesaw relations well-known in literature. Namely, the total neutrino mass matrix given in Eq. (7) will result in three light masses for active neutrinos and new heavy neutrinos, m n 1,2,3 m n 4,5,... . This form of U ν is consistent with many other well-known parameterizations [102][103][104].
The successful implementation of the ISS mechanism requires to construct the Dirac and Majorana mass matrices in terms of 3 × 3 submatrices as follows [78,80] where 0 3×3 is the 3 × 3 null matrix. Using a new notation M = M R µ −1 X M T R , we have the following ISS relations: Because m D is parameterized in terms of many free parameters, hence it is enough to choose µ X = µ 0 I 3 with µ 0 > 0, M R =M R = M 0 I 3 . We choose a simple form m D = M 0 [78,80,81]. The ISS condition |m ν | |µ X | |m D | M 0 gives √ µ 0mν M 0 0 and meaning that all of six heavy neutrinos have approximately the same mass m n i M 0 for all i > 3. Hence, the above simple forms of M R , and µ X will result in degenerate heavy neutrino masses, leading to small rates of lepton flavor violating decays of charged leptons (cLFV) that satisfy the current experimental constraints [82,83]. Relations in (12) reduce to the following simple forms: wherex ν ≡m satisfying max[|x νa |] 1 for all a = 1, 2, 3.
In our numerical analysis, we will use the best-fit values of the neutrino oscillation data [79] corresponding to the normal order (NO) scheme with m n 1 < m n 2 < m n 3 , namely In our numerical analysis, we have used the following relationŝ x ν = µ −1 0 diag m n 1 , m 2 n 1 + ∆m 2 21 , m 2 n 1 + ∆m 2 21 + ∆m 2 32 , These numerical values of neutrino masses satisfy the cosmological constraint arising from the Planck 2018 experimental data [84]: 3 i=a m na ≤ 0.12 eV. In order to simplify our numerical analysis, we assume m n 1 = 0.01 eV < m n 2 < m n 3 .
The other well-known numerical parameters are given in Ref. [79], namely Also the inverted order (IO) scheme with m n 3 < m n 1 < m n 2 can be considered in the similar way, but the qualitative results are the same with those from NO scheme, so we will not present here.
The non-unitary of the active neutrino mixing matrix I 3 − 1 2 RR † U PMNS is constrained by other phenomenological aspects such as, for instance, electroweak precision tests, cLFV decays [85][86][87], thus leading to the following constraints This constraint is consistent with the data popularly used in recent works discussion on the ISS framework [55,88,114]. The constraint on η may be more strict, depending on particular models. For example in the type III and inverse seesaw models, one has the constraint |η aa | ≤ O(10 −4 ) [36,89]. In our numerical analysis, we will choose the values satisfying |η 33 | ≤ 10 −3 , which are also consistent with the updated constraints on neutrino couplings discussed in Ref. [90,91], including the constraint from lepton universality previously discussed [113].
In the next section, we will discuss the model scalar potential as well as the one-loop contributions to ∆a µ,e .

III. HIGGS BOSONS AND ONE-LOOP CONTRIBUTIONS TO ∆a ea
The Higgs potential satisfying the symmetry mentioned in this work is where the Z 2 soft breaking term µ 2 3 Φ † 2 Φ 1 +H.c. is kept in order to generate non-zero mass for the CP-odd neutral Higgs predicted in this model. The quartic term Φ † 2 Φ 1 2 vanishes because of the Z 3 symmetry, whose assignments for the particle spectrum are given in Table   I. Consequently, the one-loop contributions to active neutrino masses discussed in Ref. [105] do not appear in this case, implying that all entries of Y h are not constrained by this condition. In general, µ 3 , and µ can be complex, while the remaining parameters µ 1,2 , µ h , and λ 1,2, 3,4,8,9,h are real [75,76]. In this work, all of these parameters and vacuum expectation values (vev) are assumed to be real, which corresponds to CP conservation.
It is easy to find the two minimization conditions of the Higgs potential, then inserting them into the Higgs potential in order to find the physical electrically charged scalars as follows: where This result is consistent with the one obtained from the Higgs potential used in Ref [75] after some transformations into a new Higgs basis and parameters.
Namely, the two Higgs doublets Φ i are changed into a new basis as follows: which can be expanded around the minimum as shown below: where v 2 1 + v 2 2 = v 2 = (246 GeV) 2 . In addition, the free parameters of the Higgs potential are transformed as follows µ 2 . . The Higgs potential after the transformation (21) takes the form [75], where the physical charged scalar states and their masses are consistent with those ones given in Eq. (20).
In our numerical calculation, we will use m 2 h + k and the mixing angle ϕ as free parameters.
Three Higgs mass parameters µ 2 2 , µ 2 h , and µ are functions of the remaining parameters. Thus, no perturbative limits on the Higgs selfcouplings are necessary to constrain the dependent functions chosen here.
From the above information we obtain all vertices providing one-loop contributions to the e b → e a γ decay rates as well as to ∆a ea . They are collected from the lepton Yukawa terms given in Eq. (6) respecting all symmetries given in Table (I), namely In order to derive the total neutrino mass matrix in the general ISS form, the leptonic Yukawa terms are rewritten as follows for the ISS mechanism discussed in this work. It is interesting to link this matrix with M D given in Eq. (8), where depending on the Z 2 charges of N IR there are two cases where non-zero couplings with only Φ 1 or Φ 2 correspond to Y N 2 = 0 or Y N 1 = 0, respectively. Namely With this new notation, the Yukawa Lagrangian for the THDM type-A is Then, all relevant couplings are given in the following lepton Yukawa interaction La- where The Feynman diagrams giving one-loop contributions to (g − 2) ea corresponding to the Lagrangian (26) are shown in Fig. 1. We do not list here the couplings of neutral gauge and Higgs bosons because they give suppressed contributions to a NP ea . In particular, the relevant couplings are only the ones with usual charged leptons s 0 e a e a and Z µ e a γ µ e a . The one-loop level contribution arising from the Z gauge boson exchange is the same as the predicted by the SM. The contributions arising from new neutral Higgs bosons are not larger than the one coming from the SM-like Higgs boson since they are suppressed by a factor of the order of O(10 −14 ), because we assume here their masses are at the TeV scale.
We will use the approximation that m 2 n i /m 2 W = 0 with i ≤ 3 and m 2 n i /m 2 Then, the contribution to a ea arising from the W exchange has the form: wheref Because |f V (x W ) + 5 12 | ≤ 5 12 and R * R T aa ≤ O(10 −3 ) 1, Eq. (28) equals to the one-loop level contribution predicted by the SM, see example in Ref. [92]: The one-loop level contribution to ∆a ea arising from the exchange of the electrically charged scalar singlet h ± k is given by [35]: where and the loop functions appearing in Eq. (31) have the forms: And the deviation from the SM is defined as follows: where a (1)SM µ (W ) 3.83 × 10 −9 [92].
Using the approximations m 2 Following Eqs. (13) and (14), we obtain that the one loop level contribution to a ea due to the exchange of h ± 1 is given by In the real part of Eq. (34), the first line corresponds to the chirally-enhanced part proportional to λ L,1 * λ R,1 whereas the second and remaining lines are the parts proportional to λ L,1 * λ L,1 and λ R,1 * λ R,1 , respectively. Now we compare our results given in Eq. (34) with the one-loop contribution due to the exchange of singly electrically charged Higgs bosons in the original versions [64] without N IR and the singlet h ± . Now we assume h ± 1 ≡ H ± in Eq. (20), corresponding to c ϕ = 1, s ϕ = 0. The absence of N IR can be conveniently derived from f H = 0 and f (x 1 ) = 0, thus implying that the only one-loop contribution from h ± 1 only consists of the first term in the third line of the real part given in Eq. (34), which is proportional to −9.19×10 −9 ×m 4 ea t 2 β c 2 ϕ /(24m 2 µ m 2 which yields a small and negative contribution to ∆a NP µ [65]. The dominant contributions to AMM arise from two-loop Barr-Zee type diagrams. The similar conclusion for the THDM type-B where t 2 β → t −2 β , hence has suppressed two-loop contributions to AMM for t β ≥ 0.4. If the mixing between two singly charged Higgs bosons vanishes, namely s α c α = 0, then all of the remaining terms in both Eqs. (34) and (35) are negative, thus not allowing to accommodate the experimental data on muon and electron anomalous magnetic moments.
Using the constraint (19) . Therefore, we will choose a safe upper bound as followŝ To avoid unnecessary independent parameters of Y h Ia without any qualitative AMM results discussed on this work and in order to cancel large one-loop contributions from these Higgs bosons to the cLFV decays e b → e a γ, we assume that The total one-loop level contribution arising from the exchange of two singly charged Higgs bosons is written as where a ea,0 (h ± ) denotes the dominant term of the chirally-enhanced part coming from the second one in the first line of the real part given in Eq. (34)  confirming that the remaining contributions are suppressed were given in Ref. [100]. Finally, a ea,0 (h ± ) = 0 only when the mixing between two SU (2) L doublets and singlets is non-zero s ϕ c ϕ = 0, and their masses are non degenerate m h ± 1 = m h ± 2 . We comment here an important property that a ea,0 (h ± ) in Eq. (39) keeps the same form for both types of THDM A and B, because these models control only the first term of λ R ia in Eq. (27), which is proportional to t β or t −1 β . This key factor controls loop contributions to AMM, which are proportional to suppressed power of t −1 β with large t β in THDM of type-B, including the model type-I. Hence, it is impossible to accommodate the AMM data in the original version of the THDM type I. On the other hand, the THDM type-A has one and two loop contributions to AMM consisting of much enhanced factors of t 2 β and t 4 β , respectively. Therefore, original versions can predict large loop contributions to AMM, provided that the new Higgs bosons are light having low masses of about few hundred GeV. These might be excluded by future collider experiments. Then, the presence of a ea,0 (h ± ) is an alternative way to explain the AMM data. Now we consider the case of f H = t −1 β , which corresponds to the models where the RH neutrino singlets couple with Φ 2 , whereas the charged leptons couple with Φ 1 [58,100]. Now, increasing values of a ea,0 (h ± ) in Eq. (39) require small t β values, thus implying that the scanning range of t β should be chosen from the lower bound t β ≥ 0.4. Numerical illustrations of a ea,0 (h ± ) are shown in Fig. 2 with fixedx ν3 = 10 −3 , i.e., 0 = |(Y N 2 ) ab | < √ 4π. In addition, x 1 Contour plots of ∆a µ (h ± ) × 10 9 and ∆a e (h ± ) × 10 13 as functions of x 1 and x 2 , wherê there are different fixed values of s ϕ , t β , Y d 1 , and Y d 2 shown in the respective panels, namely small t β = 0.5 and large t β = 10. We have checked that a ea,0 (h ± ) a ea (h ± ). Now we estimate the allowed ranges of M 0 , m h ± 1,2 , which are affected from the perturbative limit of the Yukawa coupling matrix Y N 0 defined in Eq. (24) and related with M 0 through Eqs. (11), (8), and (14). We have two different constraints corresponding to f We can choose a more strict upper bound of M 0 ≤x −1/2 ν3 s β v π/2 = 9.7s β TeV for f H = t −1 β andx ν3 = 10 −3 . Therefore forx ν3 ≤ 10 −3 and t β ≥ 0. In the last discussion we will focus on the allowed regions consisting of masses of heavy neutrinos and singly charged Higgs bosons below few TeV so that they can be detected by future colliders. The allowed regions are defined as they result in the two values of ∆a µ (h ± ) and ∆a e (h ± ) both satisfying the experimental data of the muon and electron anomalous magnetic moments within the 1 σ level, and all perturbative limits of the Yukawa couplings The region of parameter space used to scan is chosen as follows: Here we fix m n 1 = 0.01 eV corresponding to the NO scheme used in our numerical analysis. Smaller values of m n 1 will result in small allowed ranges of Y d 1,2 because of the perturbative limit affecting the relation given in Eq. (37). The scanned range of x ν3 satisfies the non-unitary constraint given in Eq. (19). The numerical results confirm that |a ea,0 (h ± )/a ea (h ± )| 1, namely 0.995 < |a e,0 (h ± )/a e (h ± )| < 1.005 and 1 ≤ |a µ,0 (h ± )/a µ (h ± )| ≤ 1.03. Hence the discussion about correlations between different contributions in Eq. (35) will not be shown. The correlations between important free parameters vs. ∆a µ (h ± ) are shown in Fig. 3. As mentioned above, large t β favors small ∆a µ (h ± ), leading to an upper bound on t β in the allowed regions, see the left-panel of Fig. 3. The scanned ranges in Eq. (41) allow all experimental ranges of ∆a e,µ . In addition, the dependence between ∆a e and ∆a µ is not interesting. We know instead that ∆a e (h ± ) ∆a e,0 (h ± ) is a function of the Yukawa coupling Y d 1 . Hence, the dependence of ∆a e (h ± ) on ∆a µ (h ± ) can be seen from the dependence of the Yukawa coupling Y d 1 in the right panel, where it is bounded in a more restrictive range than the one given in Eq. (41), see Table II, where other allowed ranges are also listed.
As we mentioned above, the dominant contributions to ∆a ea are ∆a ea,0 given in Eq. (39).
This property can be seen in Fig. 4, showing the dependence of the ratio on ∆a e (h ± ). The vertical width of the allowed region is controlled by both 1σ ranges of ∆a e (h ± ) and ∆a µ (h ± ). Also, the right panel shows the linear dependence of the allowed  [94][95][96][97][98]. Because of the sizeable mixing angle ∼ √ x ν,3 between ISS and active neutrinos ν aL , the main production channel of heavy neutrinos n I (I=4,..., 9) with mass M 0 at the LHC is via the Drell Yan anihilation process ud → n I e + a mediated by the W gauge boson in the s channel. Then the decay channel of n I can be n I → e − a W + , n a Z, n a h, where h is the standard model-like Higgs boson. The ILC can produce heavy neutrino in the processes e + e − →n a n I through the exchange of virtual W and Z bosons in the t and s-channels, respectively. The model under consideration also predicts the production channel of a heavy neutrino pair e + e − →n I n I through the virtual exchange of h ± k . In addition, the singly charged Higgs bosons in the model under consideration can be searched in a proton-proton collider through the processes , where the Yukawa couplings Y h give an important contribution [115]. The ILC can produce two singly charged Higgs bosons e + e − → h + k h − l through the n I exchange in the t-channel. Studying these processes are beyond the scope of this work, but will be investigated in more detail elsewhere. Because of the non-vanishing mixing between h ± and the singly charged components of the Higgs doublets, another decay into a CP-odd neutral Higgs boson A, such as h ± k → AW ± , can occur [112]. The correlations relating the masses with ∆a µ (h ± ) are shown in Fig. 5. We can see that GeV. The allowed regions of large m h ± 1,2 at TeV scale corresponding to (g − 2) µ experimental data may be tested indirectly through the process µ + µ − → γ * → hγ at multi-TeV muon colliders [99].
Finally, the correlations showing significant dependence of free parameters and t β are given in Fig. 6, where the allowed regions with large t β ≥ 20 require both conditions of large mixing |s 2ϕ | = 1 and largex ν3 . We obtain a small allowed range of t β ≤ 10 that was missed in Ref. [55]. Our result is consistent with the discussion corresponding to 3-3-1 models given in Ref. [101], where the THDM is embedded.
therefore M 0 may be small with small c β equivalently with large t β . We always have M 0 < 1.6 TeV for t β = 0.4 andx ν,3 = 10 −3 . Although t β ≥ 0.4 is always kept, smallerx ν,3 ≤ 10 −3 may allow large M 0 which also allow large m h ± 1,2 . On the other hand, from Eq. (39), |a ea | ∼ t β (x ν3 ) 1/2 |s 2ϕ Y d a |, hence the allowed regions are easily satisfied for small t β and large values of other parameters. Our numerical analysis shows that large x ν3 ≥ 10 −4 allows all t β ≤ 50 and all singly charged Higgs masses at the TeV scale. Contour plots of ∆a µ (h ± ) × 10 9 and ∆a e (h ± ) × 10 13 forx ν,3 = 5 × 10 −4 with small t β = 0.5 and large t β = 20 are shown in Fig. 7. x 1 Regarding the numerical analysis in the scanning range (41) with f H = −t β , the allowed ranges are tighter than the scanning regions as shown below In addition small M 0 ≤ 30 GeV is also allowed. The numerical results of the correlations be- Fig. 8, where the allowed ranges of Y d 2 correspond to a rather narrow allowed range of Y d 1 , see Eq. (43). This implies that the phenomenology of the singly charged Higgs boson at colliders related with these two couplings will have some certain relations that should be experimentally verified.
We have checked numerically that although lower bounds for allowed ranges of s ϕ and x ν3 are tiny, but never vanishes. In addition, small allowed values ofx ν3 near lower bounds require both large t β → 50 and |s 2ϕ | = 2|c ϕ s ϕ | → 1, implying the existence of ∆a e a,0 (h ± ).
Also, small allowed values of |s 2ϕ | near the lower bound require both largex ν3 → 10 −3 and t β . Illustrations for these comments are shown in Fig. 9. We see also that, the relation shown in Fig. 4 is also true with f H = −t β .
Finally, we have comments regarding the (g − 2) e data [26]: ∆a NP e = −(8.7 ± 3.6) × 10 −13 .  The allowed regions of parameter space corresponding to this data can be derived from the above described numerical analysis. Namely, excepting Y d 1 all allowed ranges of free parameters are kept unchanged to guarantee consistency with the experimental (g − 2) µ data. On the other hand, Y d 1 is changed into new values such that Y d 1 → (−8.7±3.6) (4.8±3) × Y d 2 and exclude too large values violating the perturbative limit. The two models under consideration predict allowed regions of parameter space that are different from those discussed in Ref. [59,112] for the THDM model where fermions couple to two Higgs doublets with the aligned assumption of the two respective Yukawa couplings Y 2f = ζ f Y 1f with f = u, d, , ν, where ν denote new right handed neutrinos ν aR . The large contributions to (g − 2) µ comes from the two-loop Barr-Zee contributions with the necessary condition of very light neutral CP-odd Higgs boson mass, m A ≤ 60 GeV. On the other hand, large and negative sign of (g − 2) e satisfying the 1σ experimental data comes from one-loop contribution with the condition that ζ ζ ν < 0. The model in Ref. [59] does not include the case Y 1ν = 0 and Y 2ν = 0, corresponding to f H = t −1 β mentioned in our work. In addition, the region of parameter space with ζ ν = ζ = 0 corresponding to f H = −t β is excluded by the 1σ range of (g − 2) e data.
In contrast, our models always assume that Y 1u = ζ = 0, and only one of the two Yukawa coupling matrices Y 1ν or Y 2ν being non-zero. The Yukawa couplings Y d a between only gauge singlets N a(b+3) e aR h + give main one-loop contributions to both ∆a e,µ , leading to nearly linear relations of these two quantities. Finally, our models predict regions of parameter space that successfully accommodate the experimental data on both (g − 2) e,µ anomalies without the requirement of small m A and rather light m h ± k ∼ O(10 2 ) GeV. Therefore, our models will be a another solution for the (g − 2) e,µ anomalies if a light CP odd scalar A is excluded by future experiments.

IV. CONCLUSION
In this work we have shown that the appearance of heavy ISS neutrinos and singly charged The presence of ∆a ea,0 is an alternative way to explain the AMM data if light mass ranges are excluded by future collider searches. This solution will enlarge the allowed regions of parameter spaces of the THDM which can simultaneously explain the (g − 2) e,µ data thanks to the one loop exchange of electrically charged Higgs bosons with masses within the LHC reach. The existence of ∆a ea,0 also yieds the following consequences: i) the non-zero mixing s 2ϕ = 0, ii) the non-unitary parameterx ν,3 has a lower boundx ν,3 ≥ 4 × 10 −7 for f H = t −1 β andx ν,3 ≥ O(10 −10 ) for f H = −t β , iii) and lower bounds of new heavy neutrino masses are of the order of O(100) GeV. Tiny but non vanishing values ofx ν,3 ∼ 10 −10 require very large t β and |s 2ϕ | → 1. Despite the large number of parameters, the model is economical with a small amount of BSM fields, much lower than the corresponding to several models considered in the literature. Besides that, apart from explaining the (g − 2) anomalies, the model considered in this paper can feature interesting collider signatures mainly related with electrically charged scalar and heavy neutrino production at the LHC, that can be useful to test that theory at colliders.