Analysis of the $h, H, A \to \tau \mu$ decays induced from SUSY loops within the Mass Insertion Approximation

In this paper we study the lepton favor violating decay channels of the neutral Higgs bosons of the Minimal Supersymmetric Standard Model into a lepton and an anti-lepton of different flavor. We work in the context of the most general flavor mixing scenario in the slepton sector, in contrast to the minimal flavor violation assumption more frequently used. Our analytic computation is a one-loop diagrammatic one, but in contrast to the full one-loop computation which is usually referred to the physical slepton mass basis, we use here instead the Mass Insertion Approximation (MIA) which uses the electroweak interaction slepton basis and treats perturbatively the mass insertions changing slepton flavor. By performing an expansion in powers of the external momenta in the relevant form factors, we will be able to separate explicitly in the analytic results the leading non-decoupling (constant at asymptotically large sparticle masses) and the next to leading decoupling contributions (decreasing with the sparticle masses). Our final aim is to provide a set of simple analytic formulas for the form factors and the associated effective vertices, that we think may be very useful for future phenomenological studies of the lepton flavor violating Higgs boson decays, and for their comparison with data. The accuracy of the numerical results obtained with the MIA are also analyzed and discussed here in comparison with the full one-loop results. Our most optimistic numerical estimates for the three neutral Higgs boson decays channels into $\tau$ and $\mu$ leptons, searching for their maximum rates that are allowed by present constraints from $\tau \to \mu \gamma$ data and beyond Standard Model Higgs boson searches at the LHC, are also included.

1 Introduction framework in [34] and a full-one loop diagrammatic computation in the physical SUSY particle basis was done in [37]. The issue of non-decoupling of the heavy SUSY particles in the LFVHD within this same MSSM context with NMFV was addressed numerically in [45].
In the present work we re-explore the LFV leptonic decays of the three neutral Higgs bosons, h, H, A → l klm (k = m), within the context of the MSSM with NMFV, and calculate their partial widths at the one-loop level with general slepton flavor mixings. These mixings are parametrized by means of a complete set of slepton flavor mixing dimensionless parameters, δ AB mk with AB = LL, RR, LR, RL, and flavor indices m, k = 1, 2, 3, with m = k. These parameters take into account, in a model-independent way and without any assumption on their particular origin, all the possible flavor mixings among the SUSY partners of the leptons with either left-handed or right-handed chirality. The novelty of this new computation is that we use a different technique, the so-called Mass Insertion Approximation (MIA) [49][50][51], that works with sleptons in the electroweak basis instead of the physical basis of the full one-loop computation, and treats the off-diagonal in flavor entries of the slepton squared mass matrices ∆ AB mk perturbatively, i.e., by means of mass insertions inside the propagators of the electroweak interaction sleptons eigenstates, instead of performing the exact diagonalization of the mass basis involved in the full one-loop computation. Recent studies have additionally shown that the MIA results can alternatively be also obtained if one expands properly the starting expressions in the mass basis [52,53]. The main advantage of using the MIA for the one-loop computation of the Γ(H x → l klm ) partial widths (H x = h, H, A) is clear: it provides very simple analytic formulas for the form factors involved which after a proper expansion, to be valid in the case of heavy sparticle masses of our interest here, say m SUSY O(1 TeV), can be recast in simple LFV effective vertices V eff Hxlml k , and these in turn are very useful for a simplified phenomenological study of the LFVHD rates in terms of the generic δ AB mk 's and their comparison with data. In this work, by applying the MIA at the first (linear) order in the off-diagonal mass insertions ∆ AB mk (m = k), we will compute analytically all the relevant diagrams that contribute at the one-loop level to the LFV partial widths Γ(h, H, A → l klm ). Furthermore, the MIA will also allow us to perform an analytic expansion of the involved form factors in powers of the external momenta and, in consequence, we will be able to capture analytically for the first time both the leading non-decoupling contributions of O((m h,H,A /m SUSY ) 0 ), i.e., those that go to a constant value in the asymptotically large SUSY masses limit, and the next-to-leading decoupling contributions of O(m 2 h,H,A /m 2 SUSY ), which are numerically much smaller than the leading ones but they turn out to play an important role for some of the studied cases of the flavor mixings. A few comments and estimates will also be done for the next-to-leading decoupling contributions of O(M 2 W /m 2 SUSY ), which are numerically very tiny. In this work we will also include a numerical computation of the LFVHD rates with the MIA for the case of most interest at present, h, H, A → τμ, which will be systematically compared with the full one-loop results to be able to conclude on the goodness of this approximation, the MIA, and its range of applicability. Finally we will conclude with simple analytic formulas for the useful LFV effective vertices, V eff Hxτ µ , and with a numerical estimate of the maximum expected BR(h, H, A → τμ) rates that are allowed by the present experimental constraints from τ → µγ [54] and by the ATLAS and CMS searches for neutral Higgs bosons beyond the SM [55,56]. This numerical study will be performed in terms of the most relevant model parameters, emphasizing which flavor mixing parameters will be most efficiently tested at future colliders.
The paper is organized as follows. In section 2 we summarize the relevant aspects of the MSSM with general sfermion flavor mixing and present the chosen scenarios for our numerical estimates. Section 3 contains our analytic computation of the LFVHD widths within the MIA. We select and compute the relevant one-loop diagrams and derive the form factors for LFVHD, their proper expansions, and the corresponding effective vertices. Section 4 contains all our numerical results for BR(h, H, A → τμ) and the comparison with the full one-loop predictions. The conclusions are summarized in section 5. The technicalities, including the relevant Feynman rules for the interaction vertices in the MSSM with NMFV, the analytic expressions of the form factors for each diagram, and the proper expansions of the loop integrals are collected in the appendices A, B, and C, respectively. 2 The MSSM with general flavor mixing in the charged slepton and sneutrino sectors In order to describe the MSSM with general sfermion mixing, the relevant model pieces are the superpotential and the soft SUSY-breaking Lagrangian. The superpotential of the MSSM in terms of the relevant superfields is given by: where the Yukawa couplings Y u,d,e are 3 × 3 matrices in flavor space. All indices, including the flavor ones, have been omitted in eq. (1) for simplicity. The relevant soft SUSY-breaking MSSM Lagrangian for generic sfermion mixing is: h = − sin α φ 0 1 + cos α φ 0 2 , A = − sin β χ 0 1 + cos β χ 0 2 , and use m A and tan β as input model parameters of the MSSM Higgs sector.
Since we are interested here in the Lepton Flavor Violating Higgs decays of these three neutral MSSM Higgs bosons, H x → l klm with H x = h, H, A, we will focus on sfermion mixing in the slepton sector and we will ignore the possible sfermion mixing in the squark sector. Furthermore, we will work within a general flavor mixing context at the low energies, i.e., without assuming any high-energy hypothesis for the generation of the relevant soft-breaking terms producing this slepton flavor mixing. Therefore, we will work within a NMFV framework which goes beyond the more frequently used MFV hypothesis in which the sfermion mixing is always induced by the Yukawa couplings.
The most general hypothesis for flavor mixing among sleptons assumes a mass matrix in the interaction basis that is not diagonal in the space of flavor, both for charged sleptons and sneutrinos. In the charged slepton sector the mass matrix is 6 × 6, since there are six electroweak interaction eigenstates,l L,R with l = e, µ, τ . For the sneutrinos the mass matrix is 3 × 3, since within the MSSM there are only three electroweak interaction eigenstates,ν L with ν = ν e , ν µ , ν τ .
The non-diagonal 6 × 6 squared mass matrix of sleptons when expressed in the electroweak interaction basis, that we order here as (ẽ L ,μ L ,τ L ,ẽ R ,μ R ,τ R ), is written in terms of left-and right-handed blocks M 2 l AB (A, B = L, R), which are non-diagonal 3 × 3 matrices, as follows: where: with flavor indices i, j = 1, 2, 3 running by the three generations, respectively; and (m l 1 , m l 2 , m l 3 ) = (m e , m µ , m τ ) are the lepton masses. It is worth recalling that the non diagonality in flavor comes exclusively from the soft SUSY-breaking parameters, that could be non vanishing for i = j. Specifically: the masses mL ij for the slepton SU (2) doublets, (ν LilLi ), the masses mR ij for the slepton SU (2) singlets, (l Ri ), and the trilinear couplings, A l ij . In the sneutrino sector there is a 3 × 3 squared mass matrix that, when expressed in the (ν eL ,ν µL ,ν τ L ) electroweak interaction basis, is given by: where As a consequence of the SU (2) L gauge invariance, the same soft masses mL ij enter in both the slepton and sneutrino LL mass matrices. It should be noted that if the neutrino masses and neutrino flavor mixings (oscillations) were taken into account, the soft SUSY-breaking parameters in the sneutrino sector would differ from the corresponding ones for the charged slepton sector by a rotation with the PMNS matrix. This would be somehow equivalent to what happens in the squark sector where the soft masses for the squarks of down type and those of up type differ by a rotation given by the CKM matrix. However, due to the smallness of the neutrino masses, we do not expect large effects from the inclusion of neutrino masses in the present computation and consequently we will neglect them in this work, as it is usually done in the context of the MSSM.
The general flavor mixing in the slepton sector is introduced via the non-diagonal terms in the soft breaking slepton mass matrices and trilinear coupling matrices, and these are defined here in a model-independent way in terms of a set of dimensionless parameters δ AB ij (A, B = L, R; i, j = 1, 2, 3, i = j), where L and R denote the "left-" and "right-handed" SUSY partners of the corresponding leptonic degrees of freedom, and i, j indices run over the three generations. We assume here that the δ AB ij 's provide the unique origin of LFV processes with potentially measurable rates. Specifically, we define: Some comments are in order regarding our parametrization above. First, for simplicity, in all this work we are assuming that all δ AB ij parameters are real, hence, hermiticity of the squared mass matrices implies δ AB ij = δ BA ji . Second, the diagonal entries in eq. (12) have been normalized as usually done in the literature, namely, by factorizing out the corresponding lepton Yukawa coupling: Third, it should be noted that the choice in eqs. (11), (12), and (13) is to normalize the non-diagonal in flavor entries with respect to the geometric mean of the corresponding diagonal squared soft masses. Thus, the non-diagonal LL and RR terms, with m = k, are normalized as: and However, in the case of sfermion mixing of LR (and RL) type, and taking into account that the origin of these off-diagonal mass entries is intrinsically connected to the value of the soft trilinear couplings, having dimension of mass, we find more appropriate for the purpose of this work, dealing with very large SUSY masses, to normalize them alternatively as follows: and similarly, This implies an obvious relation between δ LR mk andδ LR mk which should be kept in mind: and similarly for the RL case. Besides, if one wishes to relate the previous electroweak interaction basis and the physical mass basis one must perform the corresponding rotations:  where Rl and Rν are the respective 6 × 6 and 3 × 3 unitary rotating matrices that provide the diagonal mass-squared matrices as follows, Regarding the particle interactions that are involved in the present computation of the LFV Higgs decay widths, Γ(H x → l klm ) with k, m = 1, 2, 3, k = m, and H x = h, H, A, we have collected all the relevant Feynman rules in Appendix A, including all the needed insertions, vertices, and propagators, which we have expressed in the proper basis here. Concretely, we work with the mass basis for the external particles, H x , l k , andl m , and with the electroweak interaction basis for the internal sparticles in the loops, which from now on will be shortly denoted by:l L,R i (i = 1, 2, 3), ν i (i = 1, 2, 3),W ± ,W 3 ,B,H ± , andH 1,2 . This choice of basis is the most convenient one for the computation in the MIA, in contrast to the full one-loop computation where the physical mass basis is also usually set for the internal sleptons, sneutrinos, charginos, and neutralinos:l α (α = 1, ..., 6), ν α (α = 1, 2, 3),χ ± i (i = 1, 2), andχ 0 i (i = 1, ..., 4). Finally, to close this section of model specifications, we shortly summarize next the heavy SUSY scenarios that we work with for the estimates of this research. In order to simplify our numerical analysis, and to reduce the number of independent parameters, we define here three simplified SUSY scenarios, where the relevant parameters with mass dimensions are related to a single SUSY mass scale, m SUSY : • Equal masses scenario. In this scenario we choose the simplest case with all the relevant parameters involved set to be equal: • GUT approximation scenario. In this scenario we set an approximate GUT relation for the gaugino masses: And, for simplicity, we also relate the soft parameters and the µ parameter to a common scale by choosing: where a is a constant coefficient that we will fix to two different values for comparison, namely, a = 3 4 and 4 3 . • Generic scenario. In this scenario we wish to explore the non-equal mass generic case. Thus, we set different values for all the mass parameters involved. For the purpose of this work, the particular values of each parameter is not much relevant, but the important feature here is setting all of them to be heavy by a common m SUSY scale. Concretely, we take: For the first two scenarios that are defined above, we use a short notation for the common soft masses, namely, mL for mL = mL 1 = mL 2 = mL 3 , etc. For simplicity, in all the three scenarios we have also assumed a vanishing soft-trilinear coupling for the first generation in the charged slepton sector, i.e., A e = 0. Concerning the soft masses of the squark sector, they are indeed irrelevant for LFV processes. However, since we want to identify the discovered scalar boson with the lightest MSSM Higgs boson, we set these parameters to values which give a prediction of m h compatible with the LHC data in the mass range of 125 GeV ± 3 GeV, and fix them to the particular values mQ = mŨ = mD = A t = A b = 5 TeV in the three scenarios described just above. Besides, as already said, the other MSSM input parameters to be set in the numerical analysis are m A and tan β. Finally, regarding the δ AB ij parameters they will be taken in the conservative interval, |δ AB ij | < 1, since we wish to keep our computation in the perturbative regime. This computation will be reported in the next section.

Analytic results of the LFVHD widths in the MIA
Here we present our analytic computation of the partial widths for the LFVHD, Γ(H x → l klm ) with k, m = 1, 2, 3, k = m, and H x = h, H, A. These can be written with full generality in terms of the two form factors F (x) L,R involved in the decay amplitude of this H x (p 1 ) → l k (−p 2 )l m (p 3 ) process, as follows: Figure 1: Full one-loop diagrams for H x → l klm decays in the MSSM mass basis.
where p 1 is the ingoing Higgs boson momentum, −p 2 the outgoing momentum of the lepton l k , and p 3 the outgoing momentum of the antileptonl m , with p 1 = p 3 − p 2 . We focus here on the H x → l klm channel, but due to the fact that we work with real parameters, the predictions for the CP -conjugate channel H x → l mlk will be equal. The present computation of Γ(H x → l klm ) is performed by taking into account the following assumptions and considerations: 1) The amplitude is evaluated at the one-loop level, 2) only loops containing sleptons and sneutrinos contribute since they are the only particles propagating the LFV by means of the ∆ AB mk entries with m = k, 3) the particle content assumed here is that of the MSSM, 4) the external particles h, H, A and l k ,l m are expressed in the physical mass basis, 5) the internal loop sparticles are expressed in the electroweak interaction basis, and 6) we use the Mass Insertion Approximation [49][50][51] to describe the propagation of slepton mixing changing flavor, and work in the linear approximation for each insertion ∆ AB mk , with AB = LL, RR, LR, RL, and m = k, i.e, considering one single insertion at a time.
In order to estimate the goodness of the MIA that we use here, we have systematically compared all our results with the full-one loop results which were firstly computed in [37]. In this case, all the particles involved in the H x → l mlk decay, both external and internal to the loops, are usually expressed in the mass basis. For completeness, we display in figure 1 the eight one-loop diagrams that contribute to the full one-loop result. For our posterior numerical analysis and comparison with our computation in the MIA, we have also implemented in our code the full one-loop formulas for each of these eight diagrams contributing to F (x) L,R , which we take from [37]. From now on, we will use the labels (i), with i = 1, ..., 8 associated to each of these diagrams according to figure 1, in the comparison of the full versus MIA results.
Next, we present our computation of the form factors F L,R within the MIA. The results are presented in the following simple form, where the contribution from each single insertion is explicitly separated. In order to extract all the relevant contributions in the MIA to each of these form factors, we have selected and computed, in a systematic way, all the diagrams that dominate the decay rates in the kinematic region of our interest here, namely, for very heavy internal sparticle masses as compared to the external particle masses: m SUSY m Hx , m l k , m lm . The set of contributing one-loop diagrams in the MIA are displayed in figures 2, 3, 4, and 5, for each case with a non-vanishing insertion, ∆ LL mk , ∆ LR mk , ∆ RL mk , and ∆ RR mk , correspondingly. The labels assigned to these diagrams refer explicitly to the particular class of full one-loop diagram they should be compared with. For instance, the contributions from the MIA diagrams with labels (1a), (1b), when added, should be compared with the full diagram (1), the ones with labels (3a), (3b), when added, should be compared with (3), etc.
It should be noted that, in the scenarios that we are working with, all the sparticle masses are considered to be heavy by means of a unique common SUSY mass scale, generically called here m SUSY . In each of the three considered scenarios, the particular relation between each soft mass and m SUSY varies, but in all scenarios the sparticle masses grow linearly with the common m SUSY scale. Saying it in different words, we are integrating to one-loop order all the internal SUSY particles, considering all of them very heavy, and without keeping any of them at low energies with a fixed mass. This feature allows us to classify the various contributions from the loop diagrams in the MIA into two categories, depending on their behavior in the asymptotic limit m SUSY → ∞: 1) Contributions that go to a constant, which will be called from now on non-decoupling contributions, and 2) contributions that go to zero, which will be called from now on decoupling contributions. Furthermore, among these later we will distinguish between the dominant decoupling contributions, which decrease with m SUSY as powers of (m Hx /m SUSY ), and the subdominant decoupling contributions, which decrease with m SUSY as powers of (m EW /m SUSY ), with m EW being any of the other electroweak masses involved, namely, M W , M Z , m l k , and m lm . Here we will not include these subdominant decoupling contributions. For instance, diagrams of type (2), that would be classified as (2a), (2b), with a H xνLkνLm vertex and one insertion of ∆ LL mk type into one of the two sneutrino internal propagators, would be one of these cases, leading to contributions that are subdominant and decoupling by powers of (m EW /m SUSY ), and consequently we have not included them into our selected diagrams. Although for all the cases studied in this work, we have checked that these corrections are not relevant numerically from a phenomenological point of view, in some specific cases in which important cancellations among the leading nondecoupling contributions occur, we have found that they may play some important role in order to obtain a better convergence between the full and the MIA results. This will be commented later in our numerical analysis.
The analytic results of the form factors in eq. (32), F (x)AB L,R with AB = LL, LR, RL, RR, from all the diagrams in figures 2, 3, 4, and 5 are collected in Appendix B. The contributions from each diagram are presented separately and expressed in terms of the relevant one-loop functions, C 0 , C 2 , D 0 , andD 0 , which are given in Appendix C. Some comments are in order regarding these analytic results. First of all, it is immediate to learn that within the MIA each diagram by itself is ultraviolet finite, since the contributing loop functions, C 0 , C 2 , D 0 , andD 0 , are all convergent. This is in contrast to the full one-loop computation, where there are some diagrams that are ultraviolet divergent, more specifically, all diagrams in figure 1 except (2) and (6) and, of course, the total full one-loop result given by the sum of the eight diagrams is ultraviolet convergent [37]. Second, according to our previously explained classification into non-decoupling and decoupling contributions, we can already conclude from these analytic results which particular terms will dominate. In particular, by selecting just the contributions from the loop functions at zero external momenta, we are capturing all the non-decoupling terms, and we can already conclude that these only appear in the LL and RR form factors but not in the LR and RL ones.
By considering zero external momenta in eqs. (57) through (64), neglecting m µ , and after some algebraic simplifications due to the symmetry properties of the loop functions, we obtain for the case of our main interest here, H x → τμ with k = 3 and m = 2, the following simple results for the non-decoupling (ND) part of the form factors, which is by far the dominant part: Figure 3: Relevant one-loop diagrams within the Mass Insertion Approximation for H x → l klm decays in the MSSM electroweak interaction basis for the internal SUSY particles, with one insertion changing flavor given by × = ∆ LR mk .
It is interesting to notice that only the loop function D 0 at zero external momenta is involved in these simple expressions for the ND parts. The definition of this D 0 for arbitrary masses is given in eq. (79). It is clear that if one considers all mass parameters to be asymptotically heavy, the two functions in eqs. (33) and (34) tend to a constant value, meaning that the integration out of the heavy SUSY particles does leave as a remnant a non-vanishing value of the Γ(H x → τμ) partial widths that is constant with m SUSY if either δ LL 23 or δ RR 23 are non vanishing. We also wish to emphasize that for the particular choice of µ = m L 3 there is an important cancellation in the RR form factor between the two contributing terms in eq. (34), leading to a vanishing of the ND part in this case. It is also worth mentioning that the above analytic results at zero external momenta of eqs. (33) and (34) are in agreement with previous results obtained in the alternative framework of the effective Lagrangian approach [34].
On the other hand, the above simple expressions also tell us about how large can be this constant value as a function of the other relevant parameters, namely, tan β and m A . Indeed, these two dependencies are fully contained in the factor inside the big squared parenthesis, which can be easily derived using eq. (54) and setting s α and c α in terms of the input parameters m A and tan β, . This simple exercise gives the relevant dependence with m A and tan β in the two form factors above. We find that for the case of δ LL 23 and δ RR 23 mixings, and for generic masses, the modulo of the form factors go at large tan β as: σ By collecting all findings together, we can summarize this general power counting with all the relevant factors in the case of δ LL 23 and δ RR 23 as follows: which show, on the one hand, the expected decoupling behavior with m A in the lightest Higgs boson h case, recovering the well known feature of vanishing LFVHD rates within a SM Higgs-like scenario, and, on the other hand, the also well known feature of the enhanced heavy H and A LFVHD rates at large tan β, which grow as (tan β) 4 .
In the case of δ LR 23 and δ RL 23 mixings, the effective form factors, as we have said, decouple with the large sparticle masses, since the potential non-decoupling terms coming from the evaluation of the loop functions at zero external momenta vanish when adding the two relevant diagrams: (6c) and (8m) in the LR case and (6d) and (8m) in the RL one. In these two cases, the leading contribution then comes from the decoupling (D) terms of O(m 2 Hx /m 2 SUSY ) in the C 0 loop functions expansions: and It is remarkable that these results above for the LR and RL cases are not dependent on the lepton masses nor on µ. We also see the factorized dependence in tan β, this time inside σ (x) 1 . Thus, we can summarize the general power counting with all the relevant factors in the case ofδ LR 23 as follows: and similarly forδ RL 23 , by interchanging L by R in the formulas above. Finally, to finish this section we find illustrative to include the analytic results in the simplest scenario where all soft mass parameters are equal, i.e, the Equal masses scenario with just one SUSY scale: m SUSY = m S . In this case the formulas can be greatly simplified and they could be useful both as a reference benchmark scenario to compare with and to perform an easy phenomenological analysis. First, the form factors are expressed as: Then, by using the formulas of the loop functions in eq. (82) of Appendix C, we have found the results collected at the end of Appendix B, where we explicit the contributions from each diagram. After adding the contributions from all the diagrams, the total results of the form factors, which can be interpreted as effective LFV interaction vertices, are the following: .
We can see clearly in the total results above how relevant are the strong cancellations that occur in this Equal masses scenario among the various contributing diagrams. In fact, the behavior of the RR case at large m S changes qualitatively with respect to a generic scenario with heavy sparticles, since we find in contrast a decoupling behavior, with the form factor going as (m 2 Hx /m 2 S ), due to an exact cancellation of the non-decoupling terms in this particular case. Regarding the LL case, we find again non decoupling, and for the LR and RL cases we find decoupling as in the generic case.
Interestingly, if we keep just the leading non-decoupling terms and neglect m µ in the previous formulas for the Equal masses scenario, we are left with only one relevant form factor,F (x)LL L , and therefore the total effect of the heavy SUSY particles can be summarized in terms of a very simple effective LFV vertex given by (−igV eff Hxτ µ P L ) with: It should be noted that this is valid for all tan β values. We have further checked that when the large tan β limit is considered in this eq. (51), we find agreement with the results found from the full one-loop computation in [37]. Concretely, using eq. (35) for the lightest Higgs boson vertex we find the expected decoupling behavior in the large m A M Z limit going as (M Z /m A ) 2 , which then makes this h boson to resemble as the SM Higgs boson. We also get agreement with the expected (tan β) 2 enhanced LFV vertex [37] in the case of the H and A Higgs bosons:  (64) and also the complete one-loop formulas of [37]. The masses of the three neutral MSSM Higgs bosons, with two-loop corrections included, and their corresponding total decay widths are computed by means of the code FeynHiggs [58]. We have explicitly checked that all the numerical results for BR(H → τμ) are nearly equal to those of BR(A → τμ) and, for shortness, we will show in this section only the latter.
We start the presentation of the numerical results with the most general scenario considered along this work, the Generic scenario, in which all the SUSY mass parameters are different. We show in figure 6 the contributions of the dominant diagrams and the total one to BR(h → τμ) (left panels) and BR(A → τμ) (right panels) as functions of m SUSY , within this scenario with m A = 800 GeV and tan β = 40, for δ LL 23 = 0.5 (upper panels), δ RR 23 = 0.5 (middle panels), andδ LR 23 = 0.5 (lower panels). In each case, the other flavor changing deltas are set to zero. Since the results for δ RL 23 = 0.5 are identical to those ofδ LR 23 = 0.5, they are not shown here. The BR(h → τμ) for the LL case is displayed on the upper left panel of figure 6. First of all, we can see that each diagram contribution and the total prediction present the expected nondecoupling behavior with m SUSY , with a very good agreement between the full one-loop results and the MIA ones. The agreement is found for each diagram contribution and for the total result. It should be noted that although the non-decoupling behavior of the partial width manifests in that it goes to a constant value at large m SUSY , in the plots we see however a slight increase of the branching ratios due to the slight decrease of the total width with m SUSY . Regarding the dominant contributions, they come from diagrams 1 and 4, and we have found that they nearly cancel between each other. The rest of subdominant diagrams (3, 5, 6, and 8) are indeed important, since the remnant contributions of diagrams 1 and 4 interfere negatively with their contributions and fall down the total contribution below the diagram 3 one, what is the lowest one. Therefore, it is clear that there is in the LL case a strong cancellation among diagrams of the BR(h → τμ) that reduces the rates around three orders of magnitude, from the dominant contributions (diagrams 1 and 4) to the total one. This strong cancellation does not occur for BR(A → τμ) as we can observe on the right panel of figure 6. The dominant contribution to this process in the LL case comes  from the diagram 4, followed by far by the diagram 8. There is indeed a small negative interference between these two diagrams, resulting in total contributions slightly lower than the diagram 4 ones.
The non-decoupling behavior with m SUSY of all the contributions is also manifest, and the results in the MIA are very close to the full one-loop ones again. Now we move our attention to the RR case. The dominant and total contributions to BR(h → τμ) are depicted on the middle left panel of figure 6, in which we can see that the diagram 6 is the dominant one, followed by the diagram 5 and secondly by the diagram 8. In this case there is again a very strong cancellation among the contributions of these three diagrams, and the surviving contribution comes from the diagram 7, which reproduces very well the total result for BR(h → τμ). As in the previous cases, the agreement between the full one-loop calculation and the MIA one is very good, and all the contributions to the LFVHD partial width show a constant behavior as m SUSY grows. On the other hand, we observe on the middle right panel that for BR(A → τμ) the dominant contribution comes from the diagram 8, reproducing extremely well the total result. In this RR case there cannot be any class of interference among diagrams because the rest of contributions are at the most two orders of magnitude smaller than dominant one. All of them present also the expected non-decoupling behavior with m SUSY . We obtain again a very good agreement between the full and the MIA calculations. As we have already said in the LL case, it happens also here for the RR case that the slight increase of both branching ratios, BR(h → τμ) and BR(A → τμ), with m SUSY has not its origin in the LFV Higgs partial decay widths, since they are constant as m SUSY grows, but it is due to a small reduction of the total decay widths with m SUSY .
To end up with this Generic scenario, the results of the h → τμ and A → τμ rates for the LR case are displayed on the lower panels of figure 6. Both LFVHD rates can be understood by means of the contributions from the most relevant diagrams that are diagrams 6 and 8 in this case. The dominant non-decoupling terms, constant with m SUSY , of these two diagrams are identical in the MIA but with opposite sign. Thus, they exactly cancel and the remaining dominant decoupling contributions in the branching ratios are proportional to (m Hx /m SUSY ) 4 , what explains the final decoupling behavior of these rates with m SUSY observed in the plots. A good agreement between the full result and the MIA one is again achieved.
As main conclusions of the figure 6, we could say that in the Generic scenario the MIA approximates very well the full one-loop results, diagram by diagram and the total contributions. The LFVHD rates present a clear non-decoupling behavior with m SUSY if δ LL 23 or δ RR 23 is the responsible for the flavor mixing, whilst these rates have a strong decoupling behavior with the SUSY mass scale when theδ LR 23 orδ RL 23 is connected. In figure 7 the results for BR(h → τμ) and BR(A → τμ) as functions of m SUSY are displayed in the GUT approximation scenario with µ = 3/4 m SUSY (left panels) and µ = 4/3 m SUSY (right panels), for δ LL 23 = 0.5 (upper panels), δ RR 23 = 0.5 (middle panels), andδ LR 23 = 0.5 (lower panels). In both scenarios we have set m A = 800 GeV and tan β = 40. The first conclusion from this figure is that the LFVHD rates in this GUT approximation scenario show again a non-decoupling behavior with m SUSY in the LL and RR cases and a decoupling behavior with m SUSY in the LR case, as in the Generic scenario. We also see that the MIA works very well in all the cases LL, RR, and LR cases, reproducing accurately the results of the full one-loop computation at large m SUSY . The only exception is the prediction of BR(h → τμ), where we have found some discrepancies between the MIA and the full results in the RR case and also a little one in the LR case, being these differences larger for µ = 3/4 m SUSY than for µ = 4/3 m SUSY . We have also detected that these discrepancies are due to the fact that, in the light Higgs boson case, the missing decoupling terms in our MIA computation of the form factors of O(M 2 W /m 2 SUSY ) compete with the leading decoupling terms of O(m 2 h /m 2 SUSY ) and, for some particular cases in which there are strong cancellations among the dominant non-decoupling contributions, they may play some important role in order to obtain a better convergence between the MIA and the full results. We have also checked that this divergence      appears more pronounced where there is some degree of degeneracy among the mass parameters, as it happens partially in the GUT approximation scenario and totally in the Equal masses one.     Indeed, for this latter scenario with RR mixing, we have checked that, by means of an explicit analytic computation in the MIA of these decoupling O(M 2 W /m 2 SUSY ) contributions from the most relevant additional diagrams (see at the end of Appendix B), we achieve a better convergence between the MIA and the full results. However, we believe that is not worth including those extra terms in our general estimates here, since they are numerically extremely tiny and therefore irrelevant for the associated phenomenology.
Next we study in figure 8    close to 1, since they are still perturbative (remember eqs. (16) and (17)). Again, the observed small discrepancies in BR(h → τμ) between the MIA and the full results are due to the missing subdominant decoupling contributions of O(M 2 W /m 2 SUSY ) in our MIA calculation. The dependence of the LFVHD rates as functions of tan β is depicted in figure 9 within the Equal masses scenario with m SUSY = 5 TeV, m A = 800 GeV, and δ XY 23 = 0.5, with XY = LL, RR, LR, in each case. The full/MIA agreement in the LL and LR cases is very accurate for both LFVHD branching ratios, BR(h → τμ) and BR(A → τμ), while there is an appreciable disagreement for the RR predictions, of up to two orders of magnitude. The main reason to explain these discrepancies is that in this Equal masses scenario the cancellation among diagrams is even stronger than in the previous ones, since all of the SUSY mass parameters are identical. This strong cancellation makes that the non-decoupling dominant terms of all the diagrams completely cancel. The remnant terms in the form factors proportional to (m Hx /m SUSY ) 2 are not sufficient to reproduce the full one-loop results and then, to obtain a better convergence in this RR case, one should include the MIA subdominant decoupling contributions, proportional to (M W /m SUSY ) 2 . In order to check this expected better convergence, we have computed the most relevant diagrams providing the most important O(M 2 W /m 2 SUSY ) corrections in the MIA for this particular RR case in the Equal masses scenario. We include these analytic results at the end of Appendix B. Our numerical estimates of the LFVHD rates for this RR case after including these additional O(M 2 W /m 2 SUSY ) corrections are also displayed (in green) in figure 9, for comparison. We can clearly see that there is indeed a better convergence to the full result. However, as we have already said in all those cases where the disagreement MIA/full is clearly manifest, the predicted rates are very tiny and irrelevant for phenomenological purposes.
On the other hand, the different behaviors with tan β of the full LFVHD rates depending on  72), and knowing that, at large tan β, the total Higgs decay widths go as Γ tot (H, A) ∼ (tan β) 2 and Γ tot (h) is approximately constant with tan β. The partial widths of the h → τμ decay in the LL and RR cases, for generic SUSY masses, go as (tan β) 2 and the H, A → τμ decay widths are proportional to (tan β) 4 , therefore all the corresponding branching ratios grow as (tan β) 2 . By contrast, in the LR case, Γ(h → τμ) ∼ (tan β) −2 and Γ(H, A → τμ) are independent of tan β, thus BR(h, H, A → τμ) ∼ (tan β) −2 .
From figures 6-9 we learn that the only delta that can lead us to phenomenologically interesting LFVHD rates is δ LL 23 . In order to try to find the largest LFV Higgs branching ratios, we are going to investigate the quantities BR(h → τμ)/|δ LL 23 | 2 and BR(H, A → τμ)/|δ LL 23 | 2 that are deltaindependent when computed with the MIA. First, the contour lines of these two observables in the [m SUSY , tan β] plane are displayed in figure 10, within the Equal masses scenario with m A = 800 GeV. In both contour plots, the shaded red area is excluded by the current experimental upper limit for τ → µγ channel, BR(τ → µγ) < 4.4 × 10 −8 [54], and the shaded blue area represents the 95% C.L. excluded regions by the negative searches by ATLAS and CMS for neutral MSSM Higgs bosons decaying to a pair of τ leptons [55,56]. It is clear again the non-decoupling behavior with m SUSY of the LFVHD rates and their growth with tan β. unfortunately they are excluded by the τ → µγ upper limit and/or the ATLAS and CMS searches for MSSM Higgs bosons. The maximum values for these delta-independent rates, allowed by data, are BR(h → τμ)/|δ LL 23 | 2 ∼ 3 × 10 −8 and BR(H, A → τμ)/|δ LL 23 | 2 ∼ 5 × 10 −5 , very far away both from the current LHC sensitivity to these LFV processes [5,6].
Finally, we show in figure 11 the contour lines of BR(h → τμ)/|δ LL 23 | 2 (left panel) and BR(H, A → τμ)/|δ LL 23 | 2 (right panel) in the [m A , tan β] plane predicted in the MIA within the Equal masses scenario with m SUSY = 4 TeV, being the shaded blue area the 95% C.L. excluded regions by the negative searches by ATLAS and CMS for neutral MSSM Higgs bosons decaying to a pair of τ leptons [55,56]. The fact of fixing m SUSY = 4 TeV ensures us that the predictions are in agreement with the τ → µγ upper limit, as can be inferred from figure 10. In this case, the known decoupling behavior of BR(h → τμ) in the large m A limit is manifest on the left panel. The largest value for BR(h → τμ)/|δ LL 23 | 2 is 1 × 10 −5 , however it is again excluded by the ATLAS and CMS searches for neutral MSSM Higgs bosons. The largest h → τμ rates allowed by data are of O(10 −7 ), out of the reach of the present and next future LHC experiments. Fortunately, the prospects for H → τμ and A → τμ are much more promising, as we can see on the right panel of figure 11. The MIA predictions for BR(H, A → τμ)/|δ LL 23 | 2 are practically independent on m A and increase quadratically with tan β as expected. It reaches values, allowed by data, up to 3.5 × 10 −4 for large m A and tan β, not very far from the current LHC sensitivity. It is important to mention that our predictions of the LFVHD rates are identical for the τμ andτ µ final states, since we are assuming real δ LL 23 , and in order to compare our results with the ATLAS and CMS reported data, we have to multiply our rates by a factor of 2. Our maximum branching ratio is then of O(10 −3 ), only one order of magnitude lower than the current percent-level sensitivity achieved at the LHC [5,6].

Conclusions
In this work we have analyzed in full detail, both analytically and numerically, the decay rates of the neutral MSSM Higgs bosons into a lepton and an anti-lepton with different flavor: h, A, H → l klm (m = k). Our computation of the LFV partial widths Γ(h, A, H → l klm ) is a one-loop diagrammatic one, but different to previous analytic computations in the literature. Here it has been performed by the first time using the simple approximation provided by the MIA, which works with the electroweak interaction slepton and sneutrino eigenstates,l L,R i andν L i , with i = 1, 2, 3, and treats perturbatively the mass insertions changing lepton flavor, ∆ AB ij with AB = LL, LR, RL, RR and i = j. By using the MIA at the first order in the dimensionless parameters expansion δ AB ij , we have found compact analytic results for all the form factors involved in the LFVHD amplitudes in terms of the well known 3-and 4-point scalar one-loop integrals, and the relevant MSSM parameters, namely, the soft masses m L i , m R i , M 1 , and M 2 , the Higgs sector input mass m A , tan β, and the µ parameter. Then, by performing an expansion of the loop integrals in powers of the external momenta and keeping just the leading and next-to-leading terms, we have been able to find a set of simple analytic formulas, both for each contributing diagram and for the total sum, with all the relevant contributions explicit. These relevant contributions consist of two qualitative different parts that we have analyzed and presented separately: The leading non-decoupling contributions of O((m h,H,A /m SUSY ) 0 ) that tend to a constant value for asymptotically large m SUSY , and the next-to-leading decoupling contributions of O(m 2 h,H,A /m 2 SUSY ). At this point, we would like to emphasize that an alternative analytic computation to ours could be done by starting instead with the full analytic results of the form factors of [37], given in terms of the physical sparticle masses and rotation matrices, then performing a Taylor expansion in powers of ∆ AB mk and keeping the first order in this expansion. However, this is not an easy task since such a computation would involve a systematic Taylor expansion of all the physical slepton masses and rotation matrices elements, keeping all the relevant terms that will contribute to O(∆ AB mk ) in the form factors, and expressing them in terms of the EW basis parameters like the soft masses, etc. This kind of computation has not been completed yet, to our knowledge, for the LFV form factors of the three neutral Higgs bosons to a comparable level of our MIA computation, i.e. dealing with all the four slepton mixing cases LL, LR, RL, and RR, and keeping in the final results both the leading nondecoupling contributions of O((m h,H,A /m SUSY ) 0 ) and the next-to-leading decoupling contributions of O(m 2 h,H,A /m 2 SUSY ). We have also analyzed numerically the MIA results for the most interesting case of h, H, and A decays into τ and µ leptons. After an exhaustive comparison with the full one-loop results, we have concluded that the MIA provides indeed quite accurate predictions for the explored mixing parameters range, |δ AB 23 | < 1. We have detected only a few cases, for specific choices of the model parameters, in which there occur strong cancellations among contributing diagrams, mainly due to some degree of degeneracy in the mass parameters, where the MIA does not provide a good result as compared to the full one-loop computation. This happens for instance in the case of the Equal masses scenario with the non-vanishing flavor mixing input given by δ RR 23 . In this case, we have checked by an explicit computation that to achieve a better convergence of the MIA with the full results one must include in addition the next-to-leading decoupling contributions of O(M 2 W /m 2 SUSY ) which we have not taken into account generically in this work. Nevertheless, we wish to emphasize that this detected mismatch MIA/full is not important at all for phenomenological purposes since the predicted rates in those cases are very tiny and therefore irrelevant. Furthermore, it should be noticed that, for the heavy m A M W values considered here, it is only in the case of the lightest Higgs boson where generically the two types of next-to-leading corrections of O(M 2 W /m 2 SUSY ) and O(m 2 Hx /m 2 SUSY ) could be comparable in size and therefore, a priori, equally relevant. However, we have found that in the heavy SUSY masses scenario of our interest here, with m SUSY > 1 TeV, these corrections are below O(10 −13 ), and the maximum rates found for the lightest Higgs decays, allowed by data, are experimentally unreachable, being at most of O(10 −7 ). Hence, we have focused our interest here on the LFV heavy Higgs bosons decays.
In summary, we have presented in this work a set of simple analytic formulas for the form factors and the associated effective vertices, computed within the MIA, that we think may be very useful for future phenomenological studies of LFVHD and for their comparison with data. Finally, we have also concluded from our numerical results of the LFVHD rates, presented in contour plots in the [m A , tan β] and [m SUSY , tan β] planes, that for the most promising case of δ LL 23 mixing, one can obtain maximum allowed values (by τ → µγ experimental constraints and MSSM Higgs boson searches at the LHC) of up to BR(H, A → τ µ) ∼ 10 −3 (adding both final state τμ andτ µ rates), not far from the present experimental sensitivity accomplished at the LHC. In the case of the lightest MSSM Higgs boson h, the rates are much smaller and clearly not reachable at the LHC.

A Relevant Feynman Rules
The relevant Feynman rules for the present computation are collected in figures 12-16.
The notation and conventions are: S (x) Here and through the paper we use the short notation: s α = sin α, c α = cos α, s β = sin β, c β = cos β, t β = tan β, s α+β = sin(α + β), c α+β = cos(α + β), and t W = tan θ W . P L,R = (1 ∓ γ 5 )/2 are the usual L, R projectors. M W and M Z are the W ± and Z gauge boson masses, respectively. g and g are the gauge coupling constants of SU (2) L and U (1) Y , respectively. In the propagators, p denotes the flowing momentum and 1 denotes the identity in spinor space.

B Analytic expressions of the form factors
Here we present the analytic results of the form factors, F (x)AB L,R with AB = LL, LR, RL, RR, in eq. (32), from all the diagrams in figures 2, 3, 4, and 5. The contributions from each diagram are explicitly separated (with an obvious notation by a subscript referring to the corresponding diagram) and expressed in terms of the relevant one-loop functions, C 0 , C 2 , D 0 , andD 0 . These functions are given in Appendix C. Figure 15: Feynman rules for the relevant Higgs-ino-ino vertices.  , where we have detected that there are strong cancellations among diagrams and these contributions play a relevant role in obtaining a better convergence between the MIA and the full results. The main contributions at this order come from diagrams with two extra gaugino-Higgsino insertions in the internal fermion propagators of diagrams 5i, 5j, 6e, 8o, 8p, 8t; or one extra insertion gaugino-Higgsino and one extra of typẽ l L k −l R k in diagrams 5i, 5j; or only one extra insertion of typel L k −l R k in diagram 6e; or considering a new "type 6 like" diagram -pure Bino exchange-with vertex H x −l R k(m) −l R k(m) (no chirality flip). After this computation we have found that to include these new O(M 2 W /m 2 SUSY ) contributions into this RR form factor one should replace (F In the large tan β limit we obtain that this correction in eq. (71) grows linearly with tan β for h and quadratically for H and A. More specifically, we get for the heavy Higgs boson H (and similarly for A):F and C µ (q 1 , q 2 , m 1 , m 2 , m 3 ) = 2 i=1 q µ i C i (q 1 , q 2 , m 1 , m 2 , m 3 ) .
The particular values of the relevant loop functions for zero external momenta are the following: where a = m 2 1 , b = m 2 2 , and c = m 2 3 .
D 0 (0, 0, 0, m 1 , m 2 , m 3 , m 4 ) = 1 where a = m 2 1 , b = m 2 2 , c = m 2 3 , and d = m At non-zero external momenta all these integrals can be Taylor expanded for heavy internal particle masses as compared to the external momenta, m 2 i q 2 j , and expressed generically as their values at zero external momenta plus corrections given by functions with extra powers of the small O(q 2 j /m 2 i ) quantities. For instance, by keeping just the O(p 2 1 /m 2 i ) corrections in C 0 (p 2 , p 1 , m 1 , m 2 , m 3 ) we get: C 0 (p 2 , p 1 , m 1 , m 2 , m 3 ) = C 0 (0, 0, m 1 , m 2 , m 3 ) +b(a − c) 2 (−2ac + b(b + c))log(b) with a = m 2 1 , b = m 2 2 , c = m 2 3 , and d = p 2 1 . And similarly for other loop functions. For the present computation we have computed all the relevant Taylor expansions including the O(p 2 1 ) corrections with p 2 1 = m 2 Hx , for C 0,2 , D 0 , andD 0 , but we omit to show them here for shortness. Here we include instead just the simplest case, for illustrative purposes, that corresponds to taking all the involved SUSY masses to be equal, the so-called Equal masses scenario, keeping just the dominant and the leading subdominant contributions in the previously commented Taylor expansions. In this case, we get the following simple formulas: