Lepton flavour violation in the MSSM: exact diagonalization vs mass expansion

The forthcoming precision data on lepton flavour violating decays require precise and efficient calculations in New Physics models. In this article lepton flavour violating processes within the Minimal Supersymmetric Standard Model (MSSM) are calculated using the method based on the Flavour Expansion Theorem, a recently developed technique performing a purely algebraic mass-insertion expansion of the amplitudes. The expansion in both flavour-violating and flavour-conserving off-diagonal terms of sfermion and supersymmetric fermion mass matrices is considered. In this way the relevant processes are expressed directly in terms of the parameters of the MSSM Lagrangian. We also study the decoupling properties of the amplitudes. The results are compared to the corresponding calculations in the mass eigenbasis (i.e. using the exact diagonalization of the mass matrices). Using these methods, we consider the following processes: $\ell \to \ell' \gamma$, $\ell \to 3 \ell'$, $\ell \to 2\ell'\ell"$, $h \to \ell\ell'$ as well as $\mu \to e$ conversion in nuclei. In the numerical analysis we update the bounds on the flavour changing parameters of the MSSM and examine the sensitivity to the forthcoming experimental results. We find that flavour violating muon decays provide the most stringent bounds on supersymmetric effects and will continue to do so in the future. Radiative $\ell\to\ell^\prime\gamma$ decays and leptonic three-body decays $\ell\to 3\ell^\prime$ show an interesting complementarity in eliminating"blind spots"in the parameter space. In our analysis we also include the effects of non-holomorphic $A$-terms which are important for the study of LFV Higgs decays.


Introduction
So far, the LHC did not observe any particles beyond those of the Standard Model (SM).
Complementary to direct high energy searches at the LHC, there is a continuous effort in indirect searches for new physics (NP). In this respect, a promising approach is the search for processes which are absent -or extremely suppressed -in the SM such as lepton flavour violation (LFV) which is forbidden in the SM in the limit of vanishing neutrino masses. The experimental sensitivity for rare LFV processes such as → γ, µ → e conversion in nuclei and → µ + µ − or → e + e − will improve significantly in the near future, probing scales well beyond those accessible at foreseeable colliders. Furthermore, the discovery of the 125 GeV Higgs boson h [1,2] has triggered an enormous experimental effort in measuring its properties, including studies of its LFV decays. The most recent experimental limits on the LFV processes are given in table 2 in section 5.
Many studies of LFV processes within the MSSM (and possible extensions of it) exist (see e.g. refs.  and ref. [30] for a recent review). In this article we revisit this subject in the light of the new calculational methods which have been recently developed [31,32]. These methods allow for a systematic expansion of the amplitudes of the LFV processes in terms of mass insertions (MI), i.e. in terms of off-diagonal elements of the mass matrices. We show that a transparent qualitative behaviour of the amplitudes of the LFV processes is obtained by expanding them not only in the flavour-violating off-diagonal terms in the sfermion mass matrices but also in the flavour conserving but chirality violating entries related to the tri-linear A-terms as well as in the off-diagonal terms of the gaugino and higgsino mass matrices. This procedure is useful because in the MI approximation we work directly with the parameters of the Lagrangian and can therefore easily put experimental bounds on them. We compare the results of the calculations performed in the mass eigenbasis (i.e. using a numerical diagonalization of the slepton mass matrices) with those obtained at leading non-vanishing order of the MI approximation, in different regions of the supersymmetric parameter space and considering various decoupling limits. Of course, the MI approximation [33,34] has already been explored for many years as a very useful tool in flavour physics. However, a detailed comparison between the full calculation and the MI approximation is still lacking, partly because a fully systematic discussion of the MI approximation [31] to any order and the technical tools facilitating it [32] have not been available until recently.
Concerning the phenomenology, we summarise and update the bounds on the flavour violating SUSY parameters, show their complementarity and examine the impact of the anticipated increase in the experimental sensitivity. We investigate in detail the decay h → µτ showing the results in various decoupling limits and analyse the role of the socalled non-holomorphic A-terms [35][36][37][38][39][40][41][42], which are usually neglected in literature. We also avoid simplifying assumptions on the sparticle spectrum and assume neither degeneracies nor hierarchies among the supersymmetric particles.
This article is structured as follows: in section 2 we establish our conventions and present the results for the 2-point, 3-point, and 4-point functions related to flavour violating charged lepton interactions in the mass eigenbasis, i.e. expressed in terms of rotation JHEP06(2018)003 matrices and physical masses. Section 3 contains the formula for the decay rates of the processes under investigation. In section 4 we discuss the MI expansion and summarise important properties of the decoupling limits M SUSY → ∞ and M A → ∞. In section 5 we present the numerical bounds on LFV parameters obtained from current experimental measurements and discuss the dependence of the results on the SUSY spectrum. We also discuss the correlations between the radiative decays and the 3-body decays of charged lepton as well as the non-decoupling effects in LFV neutral Higgs decays. Finally we conclude in section 6. All required Feynman rules used in our calculations are collected in appendix A. The definitions of loop integrals can be found in appendix B. In appendix C we explain the notation for the "divided differences" of the loop functions used in the expanded form of the amplitudes. The expression for the 4-lepton box diagrams and for the MI-expanded expression of the amplitudes are given in the appendices D and E, respectively.

Effective LFV interactions
In this section we collect the analytical formula in the mass eigenbasis for flavour violating interactions generated at the one-loop level. 1 We use the notation and conventions for the MSSM as given in refs. [43,44]. 2 In our analysis, we include the so-called non-holomorphic trilinear soft SUSY breaking terms: which couple up(down)-sfermions to the down(up)-type Higgs doublets. Here, as throughout the rest of the paper, capital letters I, J = 1, 2, 3 denote flavour indices and the small letters i = 1, 2 are SU(2) L indices.

γ − − interactions
We define the effective Lagrangian for flavour violating couplings of leptons to on-shell photons as 2) The SM contribution to F JI γ is suppressed by powers of m 2 ν /M 2 W and thus completely negligible. In the mass eigenbasis the supersymmetric contributions to F JI γ come from the diagrams displayed in figure 1. Let us decompose F γ in the following way 3) 1 Note that these expressions are not valid in the flavour conserving case where additional terms should be included and renormalization is required. 2 The conventions of [43,44] are very similar to the later introduced and now widely accepted SLHA2 [45] notation, up to the minor differences summarised in the appendix A.

JHEP06(2018)003
N n (C n ) Here, V abbreviates the tree-level lepton-slepton-neutrino and lepton-sneutrino-chargino vertices, i.e. the subscripts of V stand for the interacting particles and the chirality of the lepton involved. The super-scripts refer to the lepton or slepton flavour as well as to the chargino and neutralino involved. The specific form of the chargino and neutralino vertices V L(R) is defined in appendix A and the 3-point loop functions C ij are given in appendix B. F γA (F γLB ) denotes the parts of the amplitude which is (not) proportional to the masses of fermions exchanged in the loop. F γRB can be obtained from F γLB by exchanging L ↔ R on the r.h.s. of eq. (2.4). Gauge invariance requires that LFV (axial) vectorial photon couplings vanish for onshell external particles. However, off-shell photon contributions are necessary to calculate three body decays of charged leptons. The vectorial part of the amplitude for the γ vertex can be written as where q = p I − p J and Γ JI γL is at the leading order in p 2 /M 2 SUSY momentum independent and reads Γ JI γL =  Γ JI γR can be obtained by replacing L ↔ R. Again, the loop functions C 01 , C 02 are defined in appendix B.

Z − − interactions
In order to calculate the three body decays of charged leptons as are considered in section 3.3 it is sufficient to calculate the effective Z − − interactions in the limit of vanishing external momenta. The Wilson coefficients of the effective Lagrangian for the Z coupling to charged leptons are generated at one-loop level by the diagrams shown in figure 2 and can be written as Here, Γ ZL(R) denote the contribution originating from the one-particle irreducible (1PI) vertex diagram and Σ V L(R) is the left-(right-)handed part of the lepton self-energy defined as Contrary to the left-and right-handed magnetic photon-lepton couplings, which change chirality, the Z¯ I J coupling is chirality conserving. Therefore, the Wilson coefficients of the left-handed and right-handed couplings are not related to each other but rather satisfy F IJ ZL(R) = F JI * ZL(R) . In the mass eigenbasis the vectorial part of the lepton self-energy and JHEP06(2018)003 the 1PI triangle diagrams are given by (see appendix A for definitions of vertices V ) at vanishing external momenta with obvious replacements L ↔ R for Σ JI V R , Γ JI ZR .

LFV Higgs interactions
To compactify the notation, we denote the CP-even Higgs boson decays by H K 0 →¯ I J , where, following again the notation of [43,44], H ≡ H 1 0 , h ≡ H 2 0 . As usual, we denote CP-odd neutral Higgs boson by A 0 .
In order to study h → decays precisely, we keep the terms depending on the external Higgs mass. Therefore, we assume the following effective action governing the LFV Higgslepton interaction: In addition, to calculate the µ → e conversion rate one needs to include the effective Higgsquark couplings. For this purpose, one can set all external momenta to zero and consider the effective Lagrangian However, in this article we consider only the lepton sector and therefore do not give the explicit forms of Higgs quark couplings. The relevant 1-loop expressions in the same notation as used in the current paper are given in ref. [53] and the formulae that take into account also non-decoupling chirally enhanced corrections and 2-loop QCD corrections in the general MSSM can be found in refs. [54][55][56]. 3 JHEP06(2018)003 Figure 3. Slepton-neutralino diagrams contributing to the H K 0 → I¯ J and A 0 → I¯ J decays in the MSSM (the mirror-reflected self-energy diagram is omitted). At the 1-loop level there are eight diagrams contributing to the effective lepton Yukawa couplings. The ones with slepton and neutralino exchange are displayed in figure 3, while diagrams with the chargino exchange can be obtained by the obvious replacements N → C, L →ν.
The expressions for F h and F A are obtained from 1PI triangle diagrams and the scalar part of lepton self-energies (see eq. (2.9)) while the chirality conserving parts of the selfenergies are absorbed by a field rotation required to go to the physical basis with a diagonal lepton mass matrix. Therefore, where the Z R denotes the CP-even Higgs mixing matrix (see appendix A) and the scalar self-energy contributions are evaluated at zero momentum transfer and given by: The neutralino-slepton contributions to the 1PI vertex diagrams can be written as (the symbols in square brackets denote common arguments of the 3-point functions) 4

Box contributions
4-fermion interactions are also generated by box diagrams. The corresponding conventions for incoming and outgoing particles are shown in figure 4. We calculate all box diagrams in the approximation of vanishing external momenta. The effective Lagrangian for the 4-lepton interactions involves the quadrilinear operators where X, Y stands for the chirality L or R. 5 The Wilson coefficients of these operators are calculated from the box diagrams in figure 4 and are denoted by B JIKL N XY with N = V ,S, or B JIKL T X . The operator basis in eq. (2.17) is redundant. First, we note that Second, there are Fierz relations among different operators: 5 Note that the upper L index in box formfactors denotes the sfermion flavor while the lower L subscript denotes its chirality, even if both symbols are identical. Also, recall that (¯ J σ µν PL I ) × (¯ K σµν PR L ) = 0.

JHEP06(2018)003
Eqs. (2.18) to (2.20) must be taken into account when deriving the effective Lagrangian. The case with both J = K and I = L covers the decays τ ∓ → µ ∓ e ∓ ± with = e or µ, but does not appear in µ ∓ decays. We can therefore specify to I = 3 for the effective Lagrangian. Furthermore, we can choose either (J, K) = (1, 2) or (J, K) = (2, 1) without the need to sum over both cases: the Fierz identities in eq. (2.19) permit to bring all operators into the form (e . . . τ ) × (µ . . . ) (corresponding to the case (J, K) = (1, 2)) or into an alternative form with e interchanged with µ. Thus we have with J = K and J, K, L ≤ 2, (2.21) as the four-lepton interaction in the Lagrangian. Note that the "+h.c." piece of L JK and an analogous expression for B JIKL T X .

Leptonic operators with J = K and I = L
The case J = K occurs for the decays µ ± → e ± e ± e ∓ and τ ± → ± ± ∓ with , = e, µ.
Thanks to the Fierz identities in eq. (2.19) we may restrict the operator basis to with X, Y = L, R and X = Y. (2.23)

JHEP06(2018)003
The four-lepton piece of the effective Lagrangian for the decay I∓ → J∓ J∓ L± reads: with L, J < I. (2.24) For the matching calculation it is useful to quote the tree-level matrix elements of the operators: The result is The Fierz identities further imply the equalities These operators do not appear in lepton decays, but trigger muonium-antimuonium transitions and describe muon or tau pair production in e − -e − collisions at energies far below M SUSY . Their Wilson coefficients are tiny in the MSSM.

Operators with two leptons and two quarks
The analogous Lagrangian for the 2-lepton-2-quark interactions reads Again, we consider only purely leptonic contributions here in detail and do not give explicit expressions for the 2-lepton-2-quark box diagrams. The relevant expressions in the mass eigenbasis can be found using formulae of appendix D and inserting proper quark vertices from refs. [43,44] into these.

Observables
In this section we collect the formulae for the LFV observables in terms of the effective interactions defined in section 2. All the processes listed here will be included in the future version of the SUSY FLAVOR numerical library calculating an extensive set of flavour and CP-violating observables both in the quark and leptonic sectors [66][67][68].

Radiative lepton decays: I → J γ
The branching ratios for the radiative lepton decays I → J γ are given by Here we used Γ( I → eνν) ≈ G 2 F m 5 I /(192π 3 ) for the tree-level leptonic decay width and the factors Br(µ → eνν) ≈ 1, Br(τ → eνν) = 0.1785 ± 0.0005 [69] are introduced to account for the hadronic decay modes of the τ lepton.
Even though in our numerical analyses we restrict ourselves to LFV processes, we remind the reader that the expressions for the anomalous magnetic moments and electric dipole moments of the charged leptons can be also calculated in term of the quantities defined in eq. (2.4) and read: Figure 5. Diagrams contributing to I → J K¯ L decay. I): 1PI irreducible box diagrams; II): penguin diagrams with V = Z, γ, h, H or A. For K = J crossed diagrams must be also included.

h(H) →¯ I J decays
The decay branching ratios for the CP-even and CP-odd Higgs bosons read: with F IJK h , F IJ A defined in eq. (2.14). Note that summing over lepton charges in the final state, I+ J− and J+ I− , would produce an additional factor of 2.

I → J K¯ L decays
The LFV decays of charged lepton into three lighter ones can be divided into 3 classes, depending on the flavours in the final state: (A) → : three leptons of the same flavour, i.e. µ ± → e ± e + e − , τ ± → e ± e + e − and τ ± → µ ± µ + µ − , with a pair of opposite charged leptons.
Class (C), representing a ∆L = 2 processes, is tiny within the MSSM: it could only be generated at 1-loop level by box diagrams suppressed by double flavour changes, or at the 2-loop level by double penguin diagrams involving two LFV vertices. Therefore, we will not consider these processes in our numerical analysis.
In order to calculate Br( I → J K¯ L ) we decompose the corresponding amplitude A as The relevant diagrams are displayed in figure 5. A 0 contains contributions from 4-lepton box diagrams and from penguin diagrams (including vector-like off-shell photon couplings, see eq. (2.5)) which in the limit of vanishing external momenta can be represented as the JHEP06(2018)003

4-fermion contact interactions.
A γ is the on-shell photon contribution originating from the magnetic operator (see eq. (2.2)) which has to be treated separately with more care as the photon propagator becomes singular in the limit of vanishing external momenta. We further decompose A 0 for the two cases (A) and (B) according to its Lorentz structure: with X, Y = L, R. Note that the amplitude A in general contains a second term which is obtained from the one given in eq. (3.6) by replacing (p J ↔ p J ). However, one can use Fierz identities to reduce it to the structure given in eq. (3.6). The basis of Dirac quadrilinears Γ Q is the same as the one used to decompose 4-lepton box diagrams in eq. (2.17): and Γ Q is obtained from Γ Q by lowering the Lorentz indices. The amplitudes originating from on-shell photon exchange are given by The full form of the coefficients C (A,B) N , C γ is displayed in table 1, where we compactified the expressions by using the following abbreviations for the Higgs penguin contributions: 6 Note that in eq. (3.7) and eq. (3.9) we do not explicitly display flavour indices, but they are specified in table 1.
Neglecting the lighter lepton masses whenever possible, the expression for the branching ratios can be written down as (for comparison see [23]): where N c = 1/2 if two of the final state leptons are identical (decays (A)), N c = 1 for decays (B) and X γ denotes the contribution to matrix element from the photon penguin A γ ,

JHEP06(2018)003
Decay (A) Decay (B) Table 1. Coefficients C N , C γ of eq. (3.7) and eq. (3.9) for decay types (A) and (B). B QXY ,B T X denote the irreducible box diagram contributions (see eq. (2.21)), the terms with F Z stem from the Z penguin Lagrangian (eq. (2.7)), V γ is the sum of the vector-like photon contributions (eq. (2.5)), Higgs contributions are defined in eq. (3.10) and the coefficients F γ of the magnetic operator are defined in eq. (2.2).
including also its interference with the A 0 part of the amplitude (m denotes the mass of the heaviest final state lepton)

µ → e conversion in nuclei
The full 1-loop expressions for the µ → e conversion in Nuclei depend on both the squark and slepton SUSY breaking terms. Thus, in principle the resulting upper bounds on the slepton mass insertions to some extent depend on the squark masses. Therefore, we do not include µ → e conversion in nuclei in our numerical analysis. 7 However, for completeness we collect here the complete set of formulae required to calculate the rate of this process.

JHEP06(2018)003
µ → e conversion in nuclei is produced by the dipole, the vector, and the scalar operators already at the tree level [71]. Following the discussion of ref. [72] we use the effective Lagrangian where N = V, S and X, Y = L, R with the operators defined as (3.14) Using the notation introduced in previous sections, the corresponding Wilson coefficients can be expressed as For this process, a Lagrangian involving only quark, lepton and photon fields is not sufficient. Instead, an effective Lagrangian at the nucleon level containing proton and neutron fields is required. It can be obtained in two steps. First, heavy quarks are integrated out. This results in a redefinition of the Wilson coefficient of the gluonic operator [73] with an analogous equation for C R gg . Second, the resulting Lagrangian is matched at the scale of µ n = 1 GeV to an effective Lagrangian at the nucleon level. Following [74] the transition rate Γ N µ→e = Γ(µ − N → e − N ) can then be written as

JHEP06(2018)003
where p and n denote the proton and the neutron, respectively. The effective couplings in eq. (3.17) can be expressed in terms of our Wilson coefficients as with analogous relations for L ↔ R. The Wilson coefficients in eqs. (3.18) and (3.19) are to be evaluated at the scale µ n . The nucleon form factors for vector operators are fixed by vector-current conservation, i.e. f (u) Sp/n can be borrowed from a lattice calculation [78]. 8 In summary, one has (3.20) The form factor for the gluonic operator can be obtained from a sum rule. In our normalisation The quantities D N , S (p/n) N , and V (p/n) N in eq. (3.17) are related to the overlap integrals [81] between the lepton wave functions and the nucleon densities. They depend on the nature of the target N . Their numerical values can be found in ref. [71]: for gold and aluminium, respectively.
Finally, the branching ratio is defined as the transition rate, (see eq. (3.17)), divided by the capture rate, the latter given in ref. [82]:

JHEP06(2018)003 4 Mass eigenstates vs. mass insertions calculations
For each process, we have given the exact one-loop expressions calculated in the mass eigenbasis (ME). These formulae are compact and well suited for numerical computations, however, do not allow for an easy understanding of the qualitative behaviour of the LFV amplitudes for various choices of the MSSM parameters. Therefore, in this section we expand the Wilson coefficients in terms of the "mass insertions", defined as the off-diagonal elements (both flavour violating and flavour conserving) of the mass matrices. Such an expansion allows us to: • Recover the direct analytical dependence of the results on the MSSM Lagrangian parameters.
• Prove analytically the expected decoupling features of the amplitudes in the limit of a heavy SUSY spectrum. In the case of Higgs boson decays, we also identify explicitly the terms decoupling only with the heavy CP-odd Higgs mass M A (which also determines the heavy CP even and the charged Higgs masses). The decoupling properties also serve as an important cross-check of the correctness of our calculations.
• Test the dependence of the results on the pattern of the MSSM spectrum and the size of the mass splitting between SUSY particles.
• Better understand the possible cancellations between various types of contributions and correlations between different LFV processes.
The mass insertion expansion in flavour off-diagonal terms has been used for a long time in numerous articles on the subject. However, often various simplifying assumptions have been made, i.e. some terms have been neglected or a simplified pattern of the slepton spectrum was considered. This is understandable as a consistent MI expansion of the amplitudes for the LFV processes in the MSSM, mediated by the virtual chargino and neutralino exchanges, is technically challenging. The standard approach used in literature is to calculate diagrammatically the LFV amplitudes with the "mass insertions" treated as the new interaction vertices. We follow the common practice and normalise such slepton mass insertions to dimensionless "∆-parameters": 9 where M 2 LL , M 2 RR , A l , A l are the slepton soft mass matrices and trilinear terms. As lepton flavour violation is already strongly constrained experimentally, it is sufficient to expand the amplitudes up to the first order in flavour-violating ∆'s. For instance, the JHEP06(2018)003 effective vertices listed in section 3 take the schematic form: The MSSM contributions to F LL , . . . , F BLR can be classified according to their decoupling behaviour, distinguishing the following types (M denotes the average SUSY mass scale): 1. Effects related to the diagonal trilinear slepton soft terms or to the off-diagonal elements of supersymmetric fermion mass matrices, decoupling as v 2 /M 2 .
2. Effects related to the external momenta of the (on-shell) Higgs or Z 0 bosons, decoupling as M 2 h /M 2 or M 2 Z /M 2 (we did not include the M Z dependence as it is not necessary for the considered processes).
3. Non-decoupling effects related to the 2HDM structure of the MSSM. Such contributions are constant in the limit of a heavy SUSY scale M but, in case of the SM-like Higgs boson h, decouple with the CP-odd Higgs mass like v 2 /M 2 A (the effective couplings of heavier H, A bosons do not exhibit such a suppression). They are proportional either to the lepton Yukawa couplings or to the non-holomorphic A l terms.
The structure of the box diagrams is more complicated as they carry 4 flavour indices. Their MI expansion is given in appendix E.5. All box diagram contributions decouple at least as v 2 /M 2 .
Calculating consistently the quantities F LL , . . . , F BLR to the order v 2 /M 2 is not trivial for chargino and neutralino contributions. If the MI expansion is used only for the sfermion mass matrices but the calculations for the supersymmetric fermions are done in the mass eigenbasis, the direct dependence on the Lagrangian parameters is hidden and the decoupling properties of the amplitude cannot be seen directly. However, one can also treat the off-diagonal entries of the chargino and neutralino mass matrices as "mass insertions". With such an approach, the final result is expressed explicitly in terms of Lagrangian parameters, but the computations can get very complicated. At the order v 2 /M 2 one needs to include diagrams with all combinations of two fermionic mass insertions (each providing one power of v/M 1 , v/M 2 or v/µ) or flavour diagonal slepton terms originating from trilinear A-terms (providing powers of vA l /M 2 , vA l /M 2 ). Thus, to obtain an expansion of the F 's in eq. (4.2), one needs to formally go to the 3rd order of MI expansion, adding all diagrams with up to two flavour conserving and one flavour violating mass insertion. Therefore, the number of diagrams grows quickly with the order of the expansion and such a method is tedious and prone to calculational mistakes.
In our paper, we employ a recently developed technique using a purely algebraic MI expansion of the ME amplitudes listed in section 3, without the need for direct diagrammatic MI calculations ("FET theorem") [31], automatised in the specialised MassToMI Mathematica package [32,83]. The use of this package and full automation of the calculations JHEP06(2018)003 allows us to perform the required 3rd order MI expansion for a completely general SUSY mass spectrum, without making any simplifying assumptions. Such a result would be very difficult to obtain diagrammatically, as in the intermediate steps of the calculations (before accounting for the cancellations and simplifications between various contributions) the expressions may contain up to tens of thousand terms, even if the final results collected in appendix E are again relatively compact. In detail: • We perform the expansion always up to the lowest non-vanishing order in the slepton LFV terms, taking into account the possible cancellations. Compared to previous analyses, we consider the non-holomorphic trilinear soft terms as well.
• In the MI expanded expressions we include all terms decreasing with the SUSY mass scale as v 2 /M 2 SUSY (or slower), where M SUSY denotes any of the relevant mass parameters in the MSSM Lagrangian (apart from the soft Higgs mass terms): diagonal soft slepton masses, gaugino masses M 1 , M 2 or the µ parameter.
• We do not assume degeneracy or any specific hierarchy for the sleptons, sneutrinos or supersymmetric fermion masses.
• In calculating the LFV Higgs decays we keep the leading terms in the external Higgs boson mass (m 2 h /M 2 SUSY ).
The full set of the expanded expressions in the MI approximation for the photon, Z 0 and CP-even Higgs leptonic penguins and for the 4-lepton box diagrams is collected in appendix E. We illustrate the accuracy of the derived MI formulae in figure 6. The plots show the ratio of the MI expanded couplings over the ones obtained in the mass eigenbasis with exact diagonalization. For this purpose, we start from the following setup where all mass parameters are given in GeV: Next, to see the decoupling effects we scale this spectrum uniformly up to slepton masses of 2 TeV. For each of the six penguin Wilson coefficients describing the transition between 2nd and 3rd generation, as a function of the average slepton mass. The accuracy of left-handed (right-handed) Wilson coefficients is illustrated with red(blue) lines. As can be seen from figure 6, the accuracy of MI expanded amplitudes is very good even for light SUSY particles and for M SUSY > 500 GeV always better than 95%.  Many analyses published to date for simplicity did not include the complete set of the contributions scaling like v/M order and/or assumed a partially or fully degenerate SUSY spectrum. This procedure is inconsistent with non-zero off-diagonal elements of mass matrices, because the latter enforce unequal eigenvalues of the corresponding mass matrix. To illustrate the numerical effects arising from the incorrect neglection of SUSY mass splitting we plot the ratio of our expressions in the MI approximation for penguin Wilson coefficients calculated for degenerate slepton masses (equal to 300 GeV rescaled by a common factor; other parameters as in eq. (4.3)) and the exact mass eigenbasis formulae (calculated with non-degenerate sfermion spectrum of eq. (4.3)) in figure 6. The accuracy of left-handed (right-handed) MI expanded Wilson coefficients with degenerate slepton spectrum is shown in green(brown). In this case discrepancy is much larger, of the order of 10%-40%, and does not disappear when increasing the total SUSY scale.
Some papers on the LFV in the MSSM, like e.g. refs. [7,84], deal with general SUSY spectra. In order to compare the accuracy of the MI approximation derived in our analysis with previous works, we plotted in figure 7 the ratios of Br(τ → µγ) and Br(h → τ µ) calculated using the exact (ME) and MI-expanded formulae scanning over randomly chosen MSSM mass spectra. In particular, in figure 7 we assume tan β = 5, α−β = −π/2−π/100, m h = 125 GeV, diagonal A terms which are proportional to the lepton Yukawa couplings Br(h->τμ), MI/ME µ, M 2 , mτ L , mμ L , mτ R , mτ R ∈ (500, 1000) . (4.5) As can be seen from upper left panel of figure 7, even for ∆ 32 LL = 0.5 the accuracy of our MI expansion is better than about 15%. This can be compared with the corresponding right panel of figure 8 in ref. [7] -there the difference between MI and ME calculation for the same value of ∆ 32 LL = 0.5 is 20%-70%, also the spread of points around the parabolic shape arising from neglected (∆ 32 LL ) 2 terms is much larger, 50% against a maximum of 10% in our approach. It is worth noting that the agreement in the lower left plot of figure 7 is almost perfect everywhere, which could be attributed to the fact that terms of higher order in ∆ LR are suppressed by additional v/M powers (when A terms are scaled linearly with M , like we choose) and thus small even for large ∆ LR values.
Accuracy of our MI expansion for Br(h → τ µ) shown in the right upper panel of figure 7 is worse, up to 20%, because we include only non-decoupling LL terms in our formulae, which is not a fully satisfying approximation for a SUSY scale in the range of 500-1000 GeV. Concerning the expansion in ∆ 32 LR (where only decoupling terms contribute), we include them consistently and the accuracy is much better. This can be compared with ref. [  considering the same process. The agreement for ∆ 32 LL = 0.5 in upper left panel of figure 6 ("general scenario") in [84] is better than ours, as they consistently include all LL terms, not just non-decoupling ones. However, for ∆ 32 LR (lower left panel of figure 6 in [84]) numerical accuracy of our formulae seems to be similar or even better. In general, no significant deviations should be expected here, as for this process our approach and the analysis of ref. [84] are equivalent up to the chosen calculational technique (FET vs. diagrammatic MI calculation) and, eventually, the selection of the included or neglected contributions.

Generic bounds on LFV parameters
As outlined in the introduction, flavour violation in the charged lepton sector is strongly constrained experimentally. In table 2 we collect the current and expected future experimental bounds on the processes discussed so far.
Assuming the absence of fine-tuned cancellations between different flavour violating parameters, the order of magnitude of the bounds on a given flavour violating entry ∆ can be obtained by assuming that it is the only source of flavour violation. At the lowest order  in the MI expansion, any LFV observable X scales like ∆ 2 : where f is a known (non-negative) function of diagonal mass parameters -for any given process it can be extracted from the expanded expressions listed in appendix E. Thus, the experimental bound on ∆ from a given measurement can be written as: where by X exp we denote one of the current experimental bounds listed in section 5.1 and X future is the expected future sensitivity.

Dependence on the mass splitting
The formulae derived in the previous sections allow to analyse how the bounds on LFV mass insertions depend on the splitting between different SUSY masses. However, any process involving transition between the generations I and J depends in general, even at lowest order in the flavour violating MI's, on many mass parameters: µ, gaugino masses M 1 , M 2 , left and right diagonal slepton soft masses mẽ LI , mẽ LJ , mẽ RI , mẽ RJ , and for the Higgs decays also on M A or on α angle. To simplify the discussion, we only take into account the bounds from → γ decays, which are currently most constraining.
In figure 8 we illustrate the dependence of the upper bounds on the ∆ parameters originating from µ → eγ on the mass splitting between left and right-handed sleptons for We have chosen here an average SUSY mass scale of M = 800 GeV, higher than M = 400 GeV used in tables 3-5, to avoid the experimental bounds on slepton masses even in the case of a large splitting between the left and right-handed masses. The features of plots in figure 8 can be understood using the expanded expressions for effective photon couplings collected in appendix E.1. As an example, let us consider the interesting cancellation between different contributions in the case of ∆ 12 RR (right upper panel of figure 8). For our parameter setup, the coefficient X eµ γN 2 multiplying the RR parameter (see eq. (E.4)) can be reduced to the form where f (x L , x R ) is a known, although complicated, dimensionless, rational and logarithmic function of mass ratios whose analytical form can be obtained using eq. (E.4), the loop integrals collected in appendix B and the definitions of divided differences from appendix C. The properties of this function can be examined analytically and numerically. One finds • For x R in the wide range 0.1−4 the function f vanishes for x L ∼ 0.45 (the exact value depends only weakly on x R ). As a result, the bounds on ∆ 12 RR disappear completely for m L ∼ 0.45M .

JHEP06(2018)003
• For large values of x R > ∼ 5 the position where the function f becomes zero shifts towards bigger values of x L . In addition, in this limit f is suppressed by an overall factor 1/x R , thus the bounds on ∆ 12 RR become weaker for a larger values of x R .
• For large values of x L the function f depends on x R only. Therefore, the contour lines become horizontal.
• For small values of x L the function f behaves like 1/x L . Thus, the bounds on ∆ 12 RR become stronger.
A similar analysis can be done for the bounds on ∆ 12 LL . However, the coefficient multiplying ∆ 12 LL contains contributions from both chargino and neutralino loops and does not vanish for any mass pattern. Therefore, there is no cancellation area in the upper left panel of figure 8. In this case, the bound on ∆ 12 LL is strongest for m L ∼ M and m R < ∼ M . For the case in which the left slepton masses are much lighter or much heavier than the masses of the SUSY fermion, the bounds become weaker.
Bounds on LR parameters, both holomorphic and non-holomorphic, are typically 1-2 orders of magnitude stronger than for LL and RR ones. In this case, the coefficient X eµ γN 1 multiplying the LR terms has a much simpler functional form. Therefore, it never vanishes and in addition is explicitly symmetric (as follows from the properties of divided differences) under the exchange of slepton mass arguments, as visible in both lower panels of figure 8. Furthermore, one can see that bounds on LR parameters are strongest for m L , m R < ∼ M and become weaker when the slepton masses are much heavier than the chargino and neutralino masses. More quantitatively, X eµ γN 1 is proportional to the divided difference of the function C 12 , which for x ≡ x L = x R (corresponding to the diagonal of lower plots in figure 8) has the simple asymptotic behaviour From the form of eq. (5.7) it is immediately visible that the bounds become constant for small x and fall like 1/x 2 for large x, as illustrated in the plots. Using the formulae collected in appendix E, a similar discussion can be, if necessary, performed to explain the features or cancellation areas of other plots presented in this section. However, as the general analytical formulae in the MI approximation are rather complicated, we illustrate here other scenarios with numerical plots only. Figure 9 shows similar bounds assuming identical left-and right-handed slepton masses which however differ among the generations, so that we choose bounds on ∆ 12 LL are strongest for small splitting between slepton and SUSY fermion masses, while the bounds on ∆ 12 RR are, apart from the cancellation region, stronger for mμ < ∼ M . It is obvious from the form of X eµ γN 1 in eq. (E.4) that the bounds on the LR parameters, both holomorphic and non-holomorphic, have an identical behaviour as in the case of the m L − m R splitting plotted in figure 8, with the replacements x L ↔ x e , x R ↔ x µ .
Finally in figure 10 we assume an identical mass of m = 400 GeV for all sleptons but vary M 1 = M 2 and µ. The results are displayed as a function of x 2 = M 2 /M , x µ = µ/M (we do not plot small values of |µ| < 100 GeV which are excluded by the direct searches for charginos and neutralinos). The structure of cancellation areas is more complicated, but again the "blind spots", where the bounds on MI's disappear, exist only for ∆ 12 RR . As expected from the form of X eµ γN 1 in eq. (E.4), the bounds on ∆ 12 LR , ∆ 12 LR are at leading order independent of the µ parameter. They are also correlated with the bounds displayed in lower plots of figure 8, as for a fixed slepton mass and varied M 2 the coefficient X eµ γN 1 is now proportional to so that again the bounds saturate for small x 2 and fall like 1/x 2 2 in the opposite limit. Similar plots constraining 13

Correlations between LFV processes
The correlations between various leptonic decays, in particular radiative and 3-body charged lepton decays, are important for designing new experiments searching for the LFV phenomena. In the photon penguin domination scenario the ratio of decay rates for both JHEP06(2018)003 processes is given by the simple formula: In this case the decision which measurement is more promising depends purely on experimental accuracy achievable for each of them. However, other types of contributions, like Z-penguin and box diagrams, can modify the ratio (5.9). Such contributions may be particularly important for a "blind spot" scenario, like the weakened limit on ∆ RR for some ratios of slepton and gaugino masses.
In figure 11 we plot the quantity R defined as as a function of the SUSY mass splittings, in the same scenarios as described in figure 8 and figure 9. We assume non-vanishing ∆ 12 LL and ∆ 12 RR terms. For LR terms, both holomorphic and non-holomorphic, a photon penguin dominated scenario is always realised and R is very close to 1.
As one can see from figure 11, radiative and 3-body decays are almost always closely correlated, with R differing from 1 by a few % at most. Exceptions are only possible for parameter combinations for which Br( → γ) becomes small due to cancellations or some other type of suppression, like in scenarios with large mass splitting (compare figures 8 and 9). Simultaneously, Br( → 3 ) is given by the more complicated expression (3.11), which in the limit of small photon penguin contribution becomes the sum of positive terms and cannot vanish. Thus, although both decays are usually strongly correlated and only relative experimental sensitivities decide which of them has better chances to discover generic LFV effects mediated by the slepton sector, for some particular ranges of MSSM parameter searches for 3-body charged lepton decays are a safer choice, allowing to avoid blind spots appearing for such setups due to the suppression of → γ decay rates.

Non-decoupling effects in LFV Higgs decays
LFV Higgs decays in the SM are absent at the tree level and strongly suppressed also at the loop level. Examining LFV Higgs boson decays within the MSSM is very interesting because, contrary to other processes discussed in this paper, some contributions to the Higgs decay amplitudes proportional to the lepton Yukawa couplings or to the non-holomorphic trilinear slepton soft terms do not decouple in the limit of heavy SUSY masses and can be potentially large.
As can be seen from tables 5, for an average SUSY mass scale of M = 400 GeV and the parameter setup of eq. (5.3) the upper bounds on the flavour violating parameters from Higgs decays are much weaker than from the other processes. However, the bounds from Higgs decays on the ∆ IJ LL , ∆ IJ RR and on the non-holomorphic LR terms ∆ IJ LR do not scale like 1/M 2 . Thus, comparing the limits on ∆ 13 LR and ∆ 23 LR entries from h → and → γ decays one can check that e.g. for tan β = 20 and The bounds on ∆ IJ LL , ∆ IJ RR are obtained assuming that the flavour diagonal A l terms vanish, so that all non-decoupling LL and RR contributions are proportional to the Yukawa couplings (see eq. (E.23)). In this case the Higgs decays become most constraining for slightly higher SUSY scales, again for tan β = 20 and α angle of eq. The Higgs decays in supersymmetric extensions of the SM have already been studied e.g. in [84,[102][103][104][105][106][107][108][109]. In this section we analyse within the general MSSM the decays of the lighter CP-even Higgs boson h. The mass eigenstates formulae for the MSSM contributions to the effective leptonic Yukawa couplings of h are given in eqs. ∆ IJ LR are constrained to O(1) by the vacuum stability conditions and the requirement of the absence of charge and colour breaking (CCB) minima of the scalar potential (see e.g. discussion in [110]).
The Higgs mixing angle α is subject to strong radiative corrections from the squark sector and thus from the point of view of pure leptonic sector can be treated as a free parameter. However, the allowed values of the Higgs mixing angles α, β are limited by the existing experimental constraints (see e.g. figure 6 in appendix B of ref. [107]), thus also the overall pre-factor cos 2 (α−β) cos 2 β in eq. (5.16) can be at most O(1). Summarising, the maximal Br(h → I¯ J ) which can be generated with the non-holomorphic trilinear terms is O (10 −4 ), not much below the current experimental sensitivities collected in table 2 (including decoupling contributions does not change this conclusion even for a light SUSY spectrum [106,107]). Further searches may therefore find the effects of non-holomorphic trilinear terms or provide stricter bounds on them.
Similar analysis could be done for non-decoupling contributions proportional to ∆ IJ LL and ∆ IJ RR parameters. However, in this case non-decoupling terms are proportional also either to the diagonal A l soft terms or to lepton Yukawa couplings, so the formulae become complicated and a more involved numerical analysis is required. Terms proportional to ∆ IJ LL and ∆ IJ RR multiplied by diagonal A l terms can generate similar LFV Higgs decay rates as the flavour off-diagonal non-holomorphic A l -terms. However, assuming that all non-holomorphic terms vanish, and including only the Yukawa suppressed contributions one has a much stricter bound Br(h → I¯ J ) < ∼ 10 −4 (Y I l ) 2 in the MSSM. For a complete phenomenological analysis of LFV Higgs decays in the MSSM one would need to go beyond the one-loop analysis of this article. First, one would need to perform the matching of the MSSM on the 2HDM with generic Yukawa couplings including the resummation of the higher order chirally enhanced effects (see for example [54][55][56]). Then, one has to calculate the loop effects for flavour observables within this generic 2HDM [111].

Conclusions
New precision data in the lepton flavour sector are expected to come in the foreseeable future. The search for beyond the SM effects will require precision and efficient calculations in various BSM models. In this article lepton flavour violating processes within MSSM have been calculated using the Flavour Expansion Theorem, a recently developed new technique of a purely algebraic mass-insertion expansion of the amplitudes [31]. Both flavour-violating off-diagonal terms and flavour-conserving mass-insertions are considered. The expansion in the flavour conserving off-diagonal mass terms leads to a transparent qualitative understanding of the coefficients in front of the flavour violating mass insertions (see eq. (4.2)) in various decoupling limits. Most flavour violating one-loop amplitudes decouple as v 2 /M 2 where M is one of the soft SUSY breaking mass parameters. The exception are the Higgs flavour violating decays where the amplitudes decouple as v 2 /M 2 A . We find that our full MI approximation, both in flavour violating and flavour conserving off-diagonal mass terms is an excellent approximation to the calculations in the mass eigenstates basis for a very broad pattern of supersymmetric spectra, in particular for highly non-degenerate spectra. This is JHEP06(2018)003 useful because in the MI approximation we work directly with the Lagrangian parameters and can constrain them with experimental limits.
On the physics side, the considered processes are: → γ, → 3 , → 2 , h → as well as µ → e conversion in nuclei. The bounds on the flavour changing parameters of the MSSM have been updated and their sensitivity to the forthcoming experimental results in different channels has been discussed. We have emphasised that, given the foreseen experimental progress, precision measurements of different processes have very different potential for the discovery of supersymmetric effects. The radiative and leptonic muon decays are likely to remain the most important source of information on supersymmetric LFV. The leptonic decays play a complementary role to the radiative ones in eliminating some "blind spots" of weakly constrained by the latter LFV mass insertions. This is illustrated in sections 5.2 and 5.3. Our complete analytical MI expansion facilitates the investigation of the LFV processes when the SUSY spectra are non-degenerate and finding such "blind spots" with suppressed branching ratios and regions of correlations between various processes. This is illustrated in sections 5.2 and 5.3. The LFV Higgs decays are discussed in some detail. For the supersymmetric spectrum of order of 1 TeV, the current experimental limits on the LFV Higgs decays give several orders of magnitude weaker bounds on lepton violating MI than the radiative lepton decays. However, for the superpartner masses of several TeV Higgs decays provide stronger bounds than the latter because the bounds from Higgs decays do not scale with superpartner masses. We have also analysed the role of the so-called non-holomorphic A-terms in the flavour-violating Higgs boson decays, which can give branching ratios not much below the present experimental sensitivity.

A MSSM Lagrangian and vertices
Throughout this article we use the notation of refs. [43,44] which is very similar to SLHA2 conventions [45], up to minor differences listed in table 6.
For completeness, we collect here the definitions of the mass and mixing matrices for the supersymmetric particles and the relevant MSSM Feynman rules. The slepton and sneutrino mass and mixing matrices are defined as:
where, as usual, we use tan The neutralino and chargino mass and mixing matrices can be written down as: We also use the following abbreviation for the matrix Z R parametrizing the mixing in the CP-even Higgs sector: Below we list the vertices used in calculations of the LFV processes expressed in terms of the mixing matrices defined above.

B Loop integrals
We define the following loop integrals for 2-point and 3-point functions with non-vanishing external momenta p and q: .
In our expanded results we need only the integrals above, their derivatives and higher point 1-loop integrals calculated at vanishing external momenta. Let us define . (B.2) In common notation L 2n 3 = C 2n , L 2n 4 = D 2n , L 2n 5 = E 2n etc. For i ≥ 3 one has: (with the exception of L 2 3 ≡ C 2 having also an infinite part, which however is always cancelled out in flavour violating processes and is thus not given here explicitly).
To simplify our formulae, we use the relation In addition, we define the following integrals:

C Divided differences
The expansion of the amplitudes given in the mass eigenbasis in terms of mass insertions can be naturally expressed [31] by the so-called divided differences of the loop functions.
In case of a function of a single argument, f (x), divided differences are defined recursively as: As can be easily checked, a divided difference of order n is symmetric under permutation of any subset of its arguments. It also has a smooth limit for degenerate arguments:

JHEP06(2018)003
To compactify the formulae for functions of many arguments, we use the notation where the order of the divided difference is defined by the number of arguments inside curly brackets. Then, for example a divided difference of the 1st order in the 1st argument and of the 3rd order in the 2nd argument for the function of 3 variables, g(x, y, z), can be written down as: with squared masses in the denominator. The FET expansion works for any transition amplitude, also in the case of nonvanishing external momenta or for multi-loop calculations. However, it is particularly effective for 1-loop functions with vanishing external momenta, due to the fact that the notion of the divided differences is naturally encoded in the structure of such functions: a divided difference of a n-point scalar 1-loop function is a (n+1)-point function (see eq. 3.13 in ref. [31] for generalisation to the case of non-vanishing external momenta). Thus, for example one has . . .
We use such relations extensively to find cancellations between various terms and to identify the lowest non-vanishing order of mass insertion expansion for a given process.

D Box diagrams in the mass eigenstates basis
There are four types of box diagrams with four external leptons involving slepton (sneutrinos) and neutralinos (charginos) in the loop, displayed in figure 13. Both chargino-sneutrino and neutralino-slepton pairs contribute to diagrams A) and B), while only neutralinos (Majorana fermions) can be exchanged in the "crossed" diagrams C) and D). Using whenever necessary Fierz identities, the amplitudes describing each of the diagrams N = A, B, C, D can be brought into the form with Γ V = γ µ , Γ S = 1 and Γ T = σ µν . Note that for Γ T only the case X = Y is non vanishing. Assuming that the generic couplings for an incoming lepton I -an incoming JHEP06(2018)003 Figure 13. MSSM box diagrams with 4 external charged leptons. scalar particle S k and an outgoing fermion f i takes the form the contribution from diagram A) in figure 13) to the Wilson coefficients B QXY can be written down as: where D 0 , D 2 above are the abbreviations for 4-point loop functions with respective mass arguments,

JHEP06(2018)003
Using the same notation, the contributions from diagram B), C), D) are: The contributions to 2-quark 2-lepton operators can be obtained from diagrams A) and C) by replacing K and L with q K and q L as defined in eq. (2.28). Therefore, the expressions for B q QXY can be obtained replacing vertices of leptons K and L by the relevant quark-squark vertices. Such vertices are not listed in appendix A but can be found in refs. [43,44]. The explicit form of dd box amplitudes can be also found in appendix A.3 of ref. [112].

E Effective lepton couplings in the leading MI order
We list below the MI expanded expressions for the leptonic penguin and box diagram amplitudes. For penguins we follow the decomposition of eq. (4.2), with F XY denoting functions of flavour diagonal SUSY parameters multiplying the respective slepton mass insertions: To compactify the notation, we also introduce the abbreviation where X, Y = L or R. add to or cancel terms generated from F γLB , F γRB . Thus, in the expressions below we give the sum of both types of contributions. The chargino contributions contain only terms proportional to LL slepton mass insertions (see appendix C for the notation of divided differences and curly brackets around the function arguments) The non-vanishing neutralino contributions are:  The leading v 2 /M 2 SUSY terms in the effective Z¯ I J vertex defined in eq. (2.7), expanded to the 1st order in LFV mass insertions, depend on divided differences of scalar C 0 and C 2 3point functions. They can be expressed as higher point 1-loop functions (see appendices B and C). We give here the expressions using explicitly scalar 4-, 5-and 6-point functions D, E and F .

E.3 CP-even Higgs-lepton vertex
The dominant MI terms in the effective CP-even Higgs -lepton couplings (see eq. (2.13)) can be split into four classes, defined as (below we give the sum of neutralino and chargino contributions, the latter appearing only as single term depending on sneutrino masses in eq. (E.23) and follow notation of eq. (4.2)): 1. Contributions proportional to non-holomorphic A l trilinear terms, 10 non-decoupling for M SUSY v: IJ RR M 1 D 0 (|M 1 | , mẽ LI , mẽ RI , mẽ RJ )A II L 10 For comparison with commonly used notation of the Higgs mixing angles, note that (v1Z 2K R − v2Z 1K R )/v1 = sin(α − β)/ cos β for K = 1 cos(α − β)/ cos β for K = 2 .

E.4 CP-odd Higgs-lepton vertex
For the processes considered in this article, the contribution from the LFV CP-odd Higgslepton vertex can become important only in the case of the three body charged lepton decays and only in the limit of M SUSY v, when photon, Z 0 and box contributions decouple. Thus, we give here only the dominant non-decoupling terms for this vertex.
As for CP-odd Higgs vertices, we give the sum of the neutralino and chargino contributions, the latter appearing only as single term depending on sneutrino masses in eq. (E.28): 1. Contributions proportional to non-holomorphic A l terms:

JHEP06(2018)003
Expressions listed below are valid only for ∆L = 1 processes, i.e. excluding combinations of indices I = J, K = L or I = K, J = L -for these one would also take into account flavour conserving diagrams. As mentioned in section 3.3, we do not consider MI expanded expressions for exotic ∆L = 2 processes.
The chargino diagrams contribute significantly only to the B V LL , all other contributions are at least double Yukawa suppressed and very small. The B V LL term is: Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.