One loop contributions to neutral Higgs decay h-->mu tau

In calculating one-loop contributions to amplitudes of the lepton-flavor violating decays of the neutral Higgses (LFVHD) to different flavor charged leptons, the analytic expressions can be written in term of Passarino-Veltman functions. Then, they can be computed numerically by Looptools [1]. Another approach is using suitable analytic expressions established for just this particular case. We compare numerical results obtained from Looptools and those computed by different expressions that have been applied recently. Then we derive the preferable ones that are applicable for the large ranges of free parameters introduced in the models beyond the Standard Model (SM). For illustration, the LFVHD in a simple model, which has been discussed recently, will be investigated more precisely.


I. INTRODUCTION
In 2012, the detection of the standard model (SM) Higgs was a milestone of particle physics [2]. Experimental searches for lepton-flavor violating decays of neutral Higgses (LFVHD) have seen important new improvements recently [3,4]. Motivated by this, studying for signals of LFVHD at colliders in forthcoming years has been being given attention [5].
The one-loop contributions to LFVHD in SUSY models are usually formulated by the mass insertion method [8], which will lose accuracy when the new (SUSY) scale is not much larger than the electroweak scale. There is another way that is applied to both SUSY and non-SUSY models, in which one-loop contributions are written in terms of the PV-functions [8,13,14] before they are numerically investigated using well-known computation packages [24].
Recently, there have been efforts to find convenient analytic expressions used for calculating one-loop contributions to LFVHD in non-SUSY models, without using numerical packages. The reason is that it is advantageous for studying models with simple loops contributing to LFVHD. In addition, it is extremely useful for qualitative estimation of particular loop contributions before making concrete investigations. The mass insertion method may not work well because these models predict new relevant scales rather close to the electroweak scale. Refs. [11,17,20] did not exhaustively solve integrals in analytic formulas, although Ref. [11] did cross-check them numerically with those built from PV-functions.
According to our experience, these expressions will not work well if the loops contain two small masses of virtual particles like active neutrinos or MeV masses of the exotic neutrinos in (inverse) seesaw models. Ref. [9] tried to find final analytic forms solving all integrations, but they are only valid in very special cases, e.g., when masses of new particles are much heavier than the SM-like Higgs mass. Refs. [12,21] used directly an expression containing C 0 functions-the simplest scalar integral in the set of one-loop-three point PV-functions.
But there were no analytic formulas for C 0 introduced. It was calculated using one set of fixed values of internal masses and external momenta. Ref. [16] used some particular assumptions for evaluating approximate analytic forms of one-loop contributions to LFVHD in a radiative neutrino mass model. In an effort to investigate LFVHD in a 3-3-1 model [18], an analytic expression for C 0 was introduced, needing only reasonable conditions of very small masses of normal charged leptons e, µ and τ , corresponding to the approximate zero on-shell momenta. This result was derived as a particular case of the general expression given in [25]. It is then very easy to deduce all other one-loop-three point PV-functions, which are well-known as C-functions. One of our main purposes is proving numerically the very consistency of LoopTools [1] and these analytic expressions. We then compare them with other formulas that have been used recently. In particular, we will study formulas of C-functions introduced in [9] and [16] with the aim of finding regions of parameter space that these two expressions are still valid.
We would like to stress that, in calculating LFVHD at the one-loop level, the key problem in constructing simple analytical forms of PV-functions is that the external momentum of the SM-like Higgs bosons cannot be taken to zero. It seems more dangerous when decays of heavy neutral Higgs bosons in models beyond the SM are considered. Many analytic forms of the one-loop-two-point PV-functions, denoted B-functions, are available and very consistent with LoopTools [25,26], such as expressions given in (A4). Therefore, we will not consider them in this work. But for a C-function, even two external momenta relating with charged leptons can be taken to zero, the momentum of the external Higgs cannot be ignored when the loop contains at least one heavy virtual particle. Hence, the analytic expressions for three zero external momenta, e.g., those given in [27], cannot be applied to calculating LFVHD in general.
Another new result of this work is that we will use the analytic formula of C 0 to reinvestigate the LFVHD in a lepton-flavored dark matter model introduced in [9]. In particular, we will focus on the ranges of small masses of sleptons and neutral Majorana leptons that cannot apply to expressions used in [9].
Our paper is arranged as follows. Section II will check the consistency between numerical results given by Looptools and analytical formulas introduced in [18], concentrating on the C 0,1,2 -functions. The other scalar factors of tensor integrals are easily derived by reduction procedures. We then use the analytic formulas to compare with those used in some recent studies. Section III will restudy the LFVHD in the model given in [9] and discuss on the possibilities of detection new particles predicted by this models at LHC and future colliders.
The final section is our conclusion. The three appendices list the PV-functions discussed in this work.

II. PV-FUNCTIONS FOR CALCULATING LFVHD
A. LoopTools versus analytic expressions in [18] The PV-functions relating with one-loop contributions to LFVHD are two-and threepoint functions. Conventions for external momenta are shown in Fig. 1.
FIG. 1: Notation and directions of external momenta. The left panel is from LoopTools, where k 1 = p 1 and k 2 = p 1 + p 2 , while the right panel is from [18].
We use a prime to distinguish notation between LoopTools and analytic forms, i.e., PV functions computed by LoopTool are C ′ -functions with external momenta p ′ 1 , p ′ 2 , and p ′ 3 . The notation for analytic forms is unchanged, namely C-functions defined in Eqs. (A6) and (A9). Relations between the two expressions for the same one-loop-three-point functions are as follows. The last external momentum satisfies the on-shell condition p 2 3 = p ′2 3 = m 2 h , where m h is the mass of some neutral Higgs boson, including the SM-like Higgs boson.
For the scalar function C 0 we have where the corresponding notation is p ′ and m 3 → M 2 . This result is proved by changing the variable q = −k + p 1 between the two notations. Another proof is given in Appendix B.
We need only one tensor integral C µ for calculating LFVHD in the 't Hooft Feynman gauge. The standard definition for it according to LoopTools is where the inverses of Feynman propagators in the loops are denoted as The standard definitions for the analytic expressions are listed in Appendix A, namely 2 . Now we try to find the relation between C ′µ and C µ that have relations among momenta, (2) does not affect the integral value. Therefore, we get a new expression for C ′µ : where D ′ 0 = k 2 − m 2 2 , D ′ 1 = (k − p 1 ) 2 − m 2 1 , and D ′ 2 = (k + p 2 ) 2 − m 2 3 . Fixing m 2 = M 0 , m 1 = M 1 and m 3 = M 2 and comparing Eq. (4) with Eq. (3), we obtain an important equality As a result, the relations between scalar functions are Now we will check numerically the consistency between LoopTools and analytic expressions based on the three equalities shown in (1) and (6), where C ′ 0,1,2 (M 1 , M 0 , M 2 ) is computed by LoopTools; and C 0 , C ′ 1 = C 1 − C 2 − C 0 , and C 2 are computed using C 0,1,2 given in Appendix A.
In many models, the case of M 1 = M 2 often occurs in LFVHD calculations, where M 1 (M 2 ) may be masses of charged Higgs bosons; fermions including active neutrinos, exotic leptons or quarks; charged gauge bosons and their Goldstone bosons. Therefore, we will check with M 1 = M 2 in the two following cases: (i) all M 0,1,2 are new particles from beyond the SM, (ii) at least one of these masses is an active neutrinos mass. On the other hand some models, such as the one introduced in Ref. [9], and contain three different internal masses. Therefore we will consider two more cases: In order to estimate the discrepancy between the analytical results and LoopTools, we define the relative error as follows: The following results are presented only for the real parts of the functions C 0 , C 1 , C 2 in comparison with the corresponding ones in LoopTools. Although the imaginary parts of these functions are not shown in this paper, it is noticeable that they are in perfect agreement with LoopTools as well (with relative errors all smaller than O(10 −6 %)). The figures in this section will be plotted on a logarithmic scale. illustrate the cases of the C 1 and C 2 functions. Again, we find the same conclusion as the case of the C 0 -function.
To finish the comparison with LoopTools, we emphasize that the analytic results mentioned here can be successfully applied to calculating one-loop contributions to LFVHD of  In the next subsection, we discuss other analytic forms used for calculating LFVHD at the one-loop level. B.

Discussion of other expressions for C-functions
1. Comparison with results in Refs. [12,21] In this subsection we would like to compare the numerical results of analytic forms of C-functions in [18] with other recent expressions. Refs. [12,21] used directly a formula containing C 0 -functions but with only one set of fixed values of internal masses and external momenta. The relevant loops are defined in the formula (1) 0 = 1.941. As a result, g 1 (λ, 650GeV) = −(0.26 + 0.12λ), being consistent with the value given in [12,21].
2. Approximation of C 0 -function in Refs. [9] and [16] Now we consider special cases used in Ref. [9], where the notation of the C 0 -function is the same as that in [18], in particular C 0 (0, 0, m 0 , m 1 , m 2 ) ≡ C 0 (M 0 , M 1 , M 2 ). Apart from the approximation p 2 µ , p 2 τ ≃ 0, calculation in [9] assumed a very special limit where M 2 0,1,2 ≫ m 2 h = (125.1GeV) 2 . In our notation, the C 0 -function derived from [9] is as follows where r i ≡ M 2 i /M 2 0 , i = 1, 2. The relative difference between |C 0 | and |C ′′ 0 | is defined by the quantity |δC 0 | given in Eq. (7). The two cases of r 1 = r 2 and r 2 ≫ r 1 = 1 are shown in Fig. 5. Here we also include an approximate function for C 0 (m N , m V , m V ) used in Ref. [16]. The precise formula is collected in Appendix C. We just show the relative difference with the main analytic formula in the right panels of Fig. 5 with the dotted curves. Fig. 5 shows that two analytic forms in [16] and [18] are more consistent than the expression in [9], if internal masses are few hundred GeV. If all internal masses are as large as TeV scale or more, all three results are well consistent with the relative differences being smaller than 0.1%.
3. Approximation of C 2 -function in Ref. [16] The LFVHD was also investigated in Ref. [16] with some special conditions. In the light of today's experimental data, though the analytic expressions of one-loop contributions from diagrams with W ± mediations may give large errors compared with LoopTools, they are still applicable to diagrams with new particle mediations such as new heavy charged scalars, gauge bosons and fermions in models beyond the SM [9,17,18]. Comparing two analyses for particular diagrams in the 't Hooft Feynman gauge, e.g., Ref. [16], diagram 1 a), we can derive an approximate formula for the C 2 function denoted C ′′′ 2 . It is listed in Appendix C.
The m V now can be considered as the mass of some new particle playing the role of W ± bosons in the loops. All of the assumptions given in [16] are still valid in this case, specifically Similarly, the analytic expression for The comparison is shown in Fig. 6. We can see that for seesaw models with loops containing W ± bosons and neutral exotic leptons, the expressions in [16] are not good, with the relative discrepancy around 5%. In the models with top quarks in the loop, the relative discrepancy is better, with values smaller than 1%.
We obtain the same conclusion for other C-functions where the analytic formulas are shown in Appendix C. In general, these formulas are very well consistent with all internal masses larger than 300 GeV. The details of the model are given in [9], we collect here only the ingredients relating to the LFVHD.
The Lagrangian containing all LFVHD couplings is whereφ ℓ = iσ 2 φ * ℓ , H is the SM Higgs doublet, and V 2HDM , which is the same as the Higgs potential of two Higgs doublet models (2HDM), can be found in [9]. The slepton doublet φ ℓ and the SM Higgs doublet have the same U(1) Y charge; therefore φ ℓ can be regarded as the second Higgs doublet in the 2HDM, except that the neutral component has zero vacuum expectation value (VEV). This is similar to the case of extension the SM Higgs sector, which is one of the necessary conditions to get large Br(h → µτ ) without any inconsistencies with the experimental constraint of LFV of charged lepton decays [7,21].
Neutrino masses in this model are originally from radiative corrections [9]. We will ignore contributions from active neutrinos to LFVHD because they are much smaller than those from new LFV couplings [18]. Only charged sleptons and N involve as LFV mediators.
After symmetry breaking, these new sleptons get mass from the two mass terms of φ ℓ and φ e , as well as the part coming from the trilinear coupling µvφ e φ + ℓ / √ 2 + h.c. They are all new physics beyond the SM, leading to new free parameters of the model. The original and mass eigenstates (φ ± ℓ , φ ± e ) and (ẽ ± 1 ,ẽ ± 2 ) relate to each other through the following relations where and ∆m 2 φ ≡ m 2 φ ℓ − m 2 φe . The masses ofẽ 1,2 are The mixing angle θ can be also read as Eq. (14) implies that m 2 e 2 , θ, µ, and m 2 e 1 are not independent of each other, i.e., one of them must be treated as a function of the remaining ones. The strict decoupling condition is µ = 0 and θ = 0, ±π/2. For convenience, it is enough to assume that µ > 0, 0 ≤ sin θ ≤ 1 √ 2 , leading to the consequence that mẽ 2 ≥ mẽ 1 . The signs of µ and sin 2θ will be commented on if needed. The LFVHD amplitude contains only the functions C 0 (M, mẽ i , mẽ j ) with i, j = 1, 2 for two sleptonsẽ 1 andẽ 2 . Interestingly, the partial decay width of this decay is proportional to the following part [9], New variables are defined as x 1,2 ≡ m 2 e 1,2 /M 2 . The Br(h → µτ ) was estimated from the rate R τ ≡ Br(h → µτ )/Br(τ → µγ) which is proportional to G(x 1 , x 2 )/(F (x 1 ) − F (x 2 )) or (G(x 1 ) + G(x 2 ))/(F (x 1 ) − F (x 2 )) in the decoupling or maximal mixing limit. They are denoted by a common function r(x −1 1 , x −2 2 ). According to [9], under the constraint of Br(τ → µγ) < 4.4 × 10 −8 , the value of r(x −1 1 , x −2 2 ) should be large enough to explain the current experimental value of Br(h → µτ ). In particular, the preferable regions of parameter space are as follows. In the maximal limit, masses of the two sleptons should be degenerate. In the next section we will use many analytic expressions constructed in [9] to discuss the more interesting aspects of LFVHD; in particularly we will pay attention to the regions of parameter space with few hundred GeV masses of new particles.

B. New results for LFVHD
First, we consider the function G(x 1 , x 2 ) defined in (9), where G(x 1 ) ≡ G(x 1 , x 1 ) and G(x 2 ) ≡ G(x 2 , x 2 ), and Our numerical investigation shows that the difference between the results produced from the two analytic expressions in [9] and [18] does significantly increase with small masses of M and m 2 e 1,2 . Especially if all of them are around 300 GeV, not far away from m h = 125.1 GeV.
In the following investigation, we will use the formulas for Br(h → µτ ), Br(τ → µγ) and the deviation c γ of the hγγ coupling that are established in [9], except that the G(x 1 , x 2 )function is replaced with the accurate C 0 -function mentioned above.
Regrading the estimation of Br(h → µτ ) with the ratio R τ , which is defined as Br(h → µτ ) ≡ R τ × Br(τ → µγ) [9], seems not very good, for the following reasons. First, even with very large R τ , a tiny value of Br(h → µτ ) may correspond to a very small Br(τ → µγ), and vice versa. Second, it does not show the allowed regions satisfying the bound Br(τ → µγ) < 4.8 × 10 −8 , because this constraint may rule out the regions with large Br(h → µτ ). We will give a numerical discussion of these points after reviewing the formulas required from [9].
where the two particular forms for this BR are θ → 0 in the decoupling limit and θ → π/4 in the maximal mixing limit. The parameters µ, y Rτ , y Lµ and M are introduced in the Lagrangian (10). The Yukawa couplings can be fixed as |y Rτ y * Lµ | = 1 because they are independent of charged slepton masses. In contrast, the parameter µ affects the masses of sleptons through Eq. (13), implying that it affects G(x 1 , x 2 ). Hence we believe that µ and G(x 1 , x 2 ) do not independently affect Br(h → µτ ), as discussed in [9]. In addition, the increasing µ, corresponding to decreasing x −1 1,2 , changes all values of G(x 1 ), G(x 2 ), and G(x 1 , x 2 ). So the BR in (17) depends complicatedly on µ.
The relation between Br(τ → µγ) and Br(h → µτ ) is given by Eq. (18) contains two specific limits of decoupling and maximal mixing, which are separately considered in [9]. Here we use the ratio 1/R τ instead of R τ . Recall that these ratios cancel all Yukawa couplings appearing in both expressions of the branching ratios. If we consider simultaneously both (17) and (18), the discussion in Ref. [9] for large LFVHD is illustrated in another way, as shown in Fig. 7, where sin θ = 0.1 for the decoupling limit. The three quantities Br(h → µτ ), Br(τ → µγ) and r(x −1 1 , x −1 2 ) are represented in the same figure. We emphasize that in this investigation the µ parameter is expressed as a function of mẽ i and θ, given by (14). 1 We can see that the constraint of Br(τ → µγ) seems to favour small Br(h → µτ ) in the region with degenerate slepton masses. Fig. 7 also shows two interesting points: i) a large r(x −1 1 , x −1 1 ) does not always correspond to a large Br(h → µτ ) when the experimental constraint of Br(τ → µγ) is considered, and ii) the allowed regions with large Br(h → µτ ) are sensitive to the variation of M but seem not sensitive to the change of r(x −1 1 , x −1 2 ). Furthermore, an increasing M enhances Br(h → µτ ), but causes Br(τ → µγ) to decrease. As a result, the allowed region expands wider. These conclusions are, in general, different from those indicated in [9]. An illustration of the first point is that with x 1 ≃ x 2 , increasing r will cause a decrease in the value of Br(h → µτ ); see the lower-left and upper-right of all panels in Fig. 7. For the second point, enhancement of Br(h → µτ ) can be explained by considering formula (17). With large M > 1 TeV and x 2 ≫ x 1 ≃ 1, the approximate expressions (9) are applicable and very convenient for qualitative estimation. Fig. 7 suggests two allowed regions giving large Br(h → µτ ): (i) x 1 = x 2 ≫ 1, and (ii) x 1 → 1 while x 2 → ∞ assuming that x 2 > x 1 . To understand, using (14), with v ≃ 0.25 TeV, m 2 e i = M 2 x i , we get a formula for LFVHD derived from (15) and (17): Assuming that M and θ are fixed, in the first region with x 1 → x 2 , the total factor relating to x i goes to zero in the limit, then the BR of LFVHD will go to zero too. So the large LFVHD is not caused by this degeneration between slepton masses. Instead the enhancement of BR of LFVHD originates from the large product |M sin 2θ/(1TeV)| 2 . Comparing the upper and lower panels of Fig. 7, the effect of M is clearly illustrated, while the effect of sin θ can be seen in the left and right panels, where sin θ = 0.1 and 1/ √ 2, respectively.
In the second region, where and Both of them can be arbitrary large if x 2 is not constrained. Combining with the factor of |M sin 2θ/(1TeV)|, it is easy to derive that the value of Br(h → µτ ) will take arbitrary values with large M and nonzero sin 2θ. In contrast, in this case the Br(τ → µγ) is very suppressed, which can be explained as follows. From Eq. (18), or the precise form of the partial decay width of the LFV process τ → µγ shown in [9], we can see that with F 2 (x) given in (19). It is easy to prove that lim x 1 →1 F 2 (x 1 ) = 0 and lim x 2 →∞ F 2 (x 2 ) = 0.
Hence, if x 2 or M is large enough, the Br(τ → µγ) always satisfies the experimental bounds.
Because x 2 = m 2 e 2 /M 2 and with relation (14), a large m 2 e 2 will correspond to a large µ, leading to a very narrow allowed region with large µ values, as we will show below. In this case, the experimental data such as LFVHD and hγγ coupling will give information on the parameters of the model.
In the next section, we focus just on the small values of M below 1 TeV, where N can be detected by experiments and addressed with dark matter candidates [9]. Another reason is that small values of M can be accurately investigated using the analytic expressions we mentioned above.
In this investigation, we will combine two constraints of BR(τ → µγ) and hγγ coupling to estimate Br(h → µτ ). The free parameters chosen are M, θ, µ, and mẽ 1 , apart from fixed Yukawa couplings. As an illustration, the parameter M will be fixed in two cases: small M = 300 GeV and large M = 1 TeV. The value of θ will be chosen the same as the assumption from Ref. [9], with the two limits for decoupling sin θ = 0.1 and maximal mixing The loop-induced coupling of decay h → γγ is nonnegative and constrained by [9] 0 ≤ δ γ ≡ δ cγ c SM,γ = 1 48 × 0.81 × (µv) 2 m 2 e 1 (m 2 e 1 + √ 2µv sin 2θ) < 0.20, where we have used m 2 e 2 − m 2 e 1 sin 2θ = √ 2µv. In contrast to Ref. [9], where mẽ 2 is ignored, here we include this mass in (21). The interesting consequence is that δ γ is always positive, unlike the conclusion about the two allowed regions indicated in previous work. It is easy to see that the hγγ coupling deviation gives an upper bound on µ.
Illustrations are shown in Figs From above discussion, one can conclude that if M ≤ 1 TeV the most interesting region giving large BR(h → µτ ) corresponds to the degeneration of M and the lighter slepton, i.e., M = mẽ 1 ≪ mẽ 2 or x 1 = 1 ≪ x 2 . Now we will focus on this special case. As mentioned above, because lim x 1 →1 F 2 (x 1 ) = 0 and lim x 2 →∞ F 2 (x 2 ) = 0, resulting in very suppressed Br(τ → µγ), the constraint now comes from the Higgs coupling c γ . Illustrations are drawn in Fig. 10. Both of the coupling and decoupling limits can explain the experimental LFVHD value of 5×10 −3 when M = mẽ 1 ≥ 400 GeV. And the constraint from hγγ coupling deviation gives lower bound on these masses. The parameter µ can be determined rather strictly from Because the masses of the DM candidate N and mẽ 1 , especially M ≃ mẽ 1 , are consistent with values discussed in [9], the electroweak scale of M and mẽ 1 can give the correct relic density of DM. And this conclusion does not depend on µ. But in order to satisfy both the condition of large Br(h → µτ ) and c γ constraint, we indicate that the value of µ should not be larger than 2 TeV, which is smaller than the values used in [9] for investigating direct searches of DM, µ ≥ 5 TeV. It can be seen that the DM-nucleon scattering rate decreases with decreasing values of µ [9]. In particular, the DM-nucleon scattering is generated by radiative corrections via mediations of the virtual photon and Higgs boson. The contribution from photon mediation is significant only for lighter sleptonsẽ ± 1 , but is insensitive to µ. The scattering rate will be much smaller than the current LUX sensitivity if only photon mediation is considered and M is in range O(100 GeV) [9,30]. The contribution from the Higgs mediation can be presented by the effective coupling λ hN (0) being proportional to µ [9]. Because the scattering rate can reach close to the current LUX sensitivity for µ ≥ 5 TeV, the smaller value of µ results in a smaller scattering rate. Hence, it is harder to detect DM from DM-nucleon scattering with small µ indicated by our investigation.
The region of parameter space we discussed above, with degeneration of masses of M and mẽ 1 , is an especially interesting explanation for gamma ray peak being internal bremsstrahlung in DM annihilation through a charged t-channel mediatorẽ ± i [30]. This parameter region is also the most promoting region for finding signals of Majorana DM from planned XENON1T [31] and LUXZEPLIN [32] experiments [30].

C. Productions of new particles at colliders
As we have shown, at least the masses of theẽ ± 1 and lepton-flavored Majorana DM N can be smaller than 1 TeV. Therefore, they may be detected at current running energies of the LHC or near future e + e − colliders such as the International Linear Collider (ILC) [33,34] and the Compact Linear Collider (CLIC) [35,36]. Interestingly, the sleptons defined in the model under consideration couple to SM particles in a very similar way to the sleptons in supersymmetric (SUSY) extensions of the SM. In addition, all new particles in both kinds of models are odd under Z 2 symmetries. Hence the lightest neutral particle N is a DM candidate, and plays the same role as the lightest neutralino in the SUSY models with Rparity conservation. In general, the new particle spectrum in the considering model can be seen as a simplified version of the superpartner spectrum, which has been hunted by LHC [37], especially searches for slepton productions [38][39][40][41]. For SUSY models, current channels of experimental searches are pp →ll → (ℓ χ 0 1 )(ℓ χ 0 1 ), where ℓ ± denotes an SM lepton state: e, µ, τ . In notation of the model under consideration, these channels correspond to processes The slepton productions at the LHC happen via virtual gauge bosons, i.e., pp → γ * , Z * → ℓ 0 ℓ 0 * , ℓ + ℓ − ; or pp → W ± * → ℓ ± ℓ 0 , because the gauge bosons are always lighter than the new particles. Based on couplings of new particles at final states, we see that the signals of detection of new particles in both kinds of models, namely SUSY and the model studied in our work, are of the same order. Couplings of SUSY particles are given in detail in many textbooks, e.g., [42]. The relevant couplings of sleptons and N predicted by the model under consideration are collected in Table I Coupling Vertex slepton production from gluon fusion channel gg → h 1 0 →ẽ +ẽ− relating to a top quark loop. The above discussion shows that searches for SUSY sleptons can be applied for sleptons in the model under consideration. Current lower bounds of sleptons are a few hundred GeV, which do not exclude the light sleptons and N in the region of parameter space we discussed above. And they may be detected at the LHC [43][44][45]. For example, with the condition of very small differences between masses of stau and the lightest neutralino (not larger than 1 GeV), Ref. [43] suggested that the expected number of staus may be several hundred at 8 and 14TeV LHC run with light masses larger than 450 GeV.
The sleptons and N can also be searched for in future e + e − colliders with collision energies up to 3 TeV, such as the ILC and CLIC. Predictions for signals of sleptons and DM were indicated in SUSY models [46][47][48][49] and models with lepton-flavored DM [45,50]. Similarly to the LHC, slepton productions will be searched through s channels of e + e − → Z * , γ * → e + il − j , φ 0 ℓ φ 0 * ℓ . In contrast to the LHC, where quarks and gluons do not couple toẽ i and N, there are additional t(u) channels through the exchange of N, leading to enhancements of slepton production at the ILC. In addition, e + e − colliders give a direct channel of DM search e + e − → N Nγ, corresponding to a signal of a mono-photon plus missing energy. Another indirect search is the one-loop contribution, where sleptons and N run in loops, to lepton pair production, e + e − → l + i l − i , and multiflavor lepton final state e + e − → l + i l − j (i = j) [45]. Recent bounds of new particle masses obtained from the LEP are a few hundred GeV [45,51,52], which is consistent with results obtained from the LHC.
In summary, the lower constraints of masses of sleptons and N under recent experimental results are a few hundred GeV. The parameter region of the lepton-flavor DM we discussed in this work is still valid, and is very predictive for many future projects of (in)direct searches for these new particles.

IV. CONCLUSION
The one-loop contributions to LFV decays of neutral Higgs bosons are now very interesting in many models beyond the SM, where many new particles may inherit masses that are not far from the electroweak scale. In some models, even the top quarks can play the role of LFV mediations in the loop. These one-loop contributions can be conveniently written in terms of the one-loop-three point C-functions, before taking any approximations for more precise forms used for numerical investigations. We have shown that numerical results obtained from the analytical forms of the C-functions introduced in Ref. [18] are in great agreement with those evaluated by LoopTools. This conclusion is true for all ranges of mass values in the loops, even with loops containing active neutrino masses smaller than a few eV. We have compared this with the two other analytic approximations given in [9] and [16]. We have found that the latter two expressions are still safe with all masses in the loops large than 1 TeV for the case of studying LFVHD of the SM-like Higgs boson. But they fail with masses in the loops below a few hundred GeV. Furthermore, they can not be applied for LFVHD of new heavy neutral Higgs bosons appearing in many models beyond the SM.
However, the results in [18] still work very well.
Based on the above conclusions, the analytic formulas of C-functions given in [18] have been used to reinvestigate the LFVHD mentioned in Ref. [9], focused on the regions of small masses of Majorana dark matter M and slepton masses mẽ 1,2 . We stress that these regions could not be accurate with the approximation used in previous works. We found many interesting results that are not metioned in [9]. In particular, large Br(h → µτ ) depends strongly on M, namely it enhances with increasing values of M. Even when constraints of both Br(τ → µγ) and hγγ coupling deviation are included, the LFVHD can be arbitrary large with very large M if the following condition is satisfied: M = mẽ 1 ≪ mẽ 2 . In the case of M below 300 GeV, the large Br of LFVHD near the recent experimental report occurs only in the region having M = mẽ 1 . The BR of LFVHD is constrained by the hγγ coupling deviation, where the largest value is of order 10 −3 . With M around 1 TeV, the LFVHD constraint from experiment leads to the consequence that µ should be smaller than few TeV. The parameter region discussed in this work can be tested by the LHC and the ILC in coming years. 2 and (p 1 + p 2 ) 2 = m 2 h . In this work, with m 1 and m 2 are the respective masses of muon and tauon, and m h is the SM-like Higgs mass. The tensor integrals are where B (i) 0 and C 0,1,2 are PV-functions. It is well known that C 0,1,2 are finite, while B are divergent. We define ∆ ǫ ≡ 1 ǫ + ln 4π − γ E + ln µ 2 , where γ E is the Euler constant. The divergent parts of the B-functions can be determined as Div[B x k ln 1 − 1 where x k , (k = 1, 2) are solutions of the equation The C 0 -function was given in [18] consistent with that discussed in [28], namely where Li 2 (z) is the di-logarithm function, x 1,2 are solutions of Eq. (A5), and x 0,3 are given as In the limit p 2 1 , p 2 2 → 0 the C 1,2 -functions are The parameterization of C ′ 0 is chosen as 1 . From the equalities (p ′ 1 +p ′ 2 ) 2 = m 2 h and 2p ′ 1 .p ′ 2 = (p ′ 1 +p ′ 2 ) 2 −p ′ 1 −p ′ 2 = m 2 h − p ′ 1 − p ′ 2 , we get Comparing with C 0 shown in [18], namely Appendix C: Analytic approximation from [16] Here we list the needed approximation The approximation in some special cases are