Flavor techniques for LFV processes: Higgs decays in a general seesaw model

Lepton flavor violating processes are optimal observables to test new physics, since they are forbidden in the Standard Model while they may be generated in new theories. The usual approach to these processes is to perform the computations in the physical basis; nevertheless this may lose track of the dependence on some of the fundamental parameters, in particular on those at the origin of the flavor violation. Consequently, in order to obtain analytical expressions directly in terms of these parameters, flavor techniques are often preferred. In this work, we focus on the mass insertion approximation technique, which works with the interaction states instead of the physical ones, and provides diagrammatic expansions of the observables. After reviewing the basics of this technique with two simple examples, we apply it to the lepton flavor violating Higgs decays in the framework of a general type-I seesaw model with an arbitrary number of right-handed neutrinos. We derive an effective vertex valid to compute these observables when the right-handed neutrino masses are above the electroweak scale and show that we recover previous results obtained for low scale seesaws. Finally, we apply current constraints on the model to conclude on maximum Higgs decay rates, which unfortunately are far from current experimental sensitivities.


Contents 1 Introduction
Lepton flavor violating (LFV) processes are optimal observables to test new physics hypotheses. This is particularly true for LFV transitions in the charged sector, as they are forbidden in the Standard Model (SM), and extremely suppressed if one minimally introduces the observed light neutrino masses. Consequently, any observation of charged LFV transition would be a clear evidence of new physics beyond the SM (BSM). Moreover, since in several BSM theories this kind of processes are induced via quantum corrections with new particles running in the loops, exploring LFV transitions in the intensity frontier allows us to probe the existence of these new particles even if they are too heavy to be directly produced in any other experiment.
Many BSM theories are constructed following the same steps that succeed in the case of the SM. First, we write the most general Lagrangian with the chosen symmetries and field content, which includes some parameters that we may consider as the fundamental parameters of the theory. The easiest way of doing this is by choosing a field basis in which the conservation of the symmetries, in particular the gauge symmetries, is explicit, implying in most of the cases that gauge interaction are diagonal in this basis and that some particles are massless. Therefore, we refer to this basis as the interaction basis in general, and also as gauge or flavor basis for the particular cases of gauge or fermionic fields, respectively. As a second step, we assume that some of these gauge symmetries are spontaneously broken to a smaller gauge symmetry group, as in the case of the electroweak symmetry breaking (EWSB), providing masses to some of the fields. Nevertheless, these new mass terms do not need to be diagonal in the interaction basis, and some non-diagonal terms may appear, mixing different interaction fields. The basis in which the mass terms are diagonal is called the mass or physical basis, since it is in this basis where the parameters could be directly related to observables. The relation between the two basis is given by a series of unitary rotation matrices. In some simple cases, we may be able to find analytical expressions for these rotation matrices in terms of the original parameters in the interaction basis, however this is not always possible and we often need to use numerical methods.
The physical basis is the natural choice to compute any transition in quantum field theory (QFT), since we can properly define a loop expansion for any observable, such as the LFV processes LFV Obs. Present Upper Bounds (95% C.L.) BR(H → µe) -3.5 × 10 −4 CMS (2016) [3] BR(H → τ e) 4.7 × 10 −3 ATLAS (2019) [4] 6.1 × 10 −3 CMS (2018) [5] BR(H → τ µ) 2.8 × 10 −3 ATLAS (2019) [4] 2.5 × 10 −3 CMS (2018) [5]  we mentioned before. The reason is that, in this basis, particles have well-defined propagators, meaning that they will keep their identity unless they interact with other two -or more-particles. The resulting expressions will be given in terms of the physical parameters, i.e., physical masses and rotation matrices. Nevertheless, in many cases it is desirable to have expressions that, albeit approximate, are given directly in terms of the fundamental parameters of the interaction basis. Indeed, this is particularly interesting in the case of flavor transitions, since we can often track its origin back to few parameters in the original Lagrangian, e.g., the non-diagonal Yukawa couplings in many models.
One possibility in order to obtain expressions in terms of the fundamental parameters is to compute the amplitudes first in the physical basis, and then expand them using the relations between the two basis. An interesting possibility in this direction is given by the Flavor Expansion Theorem [1,2], which provides a recipe to translate analytically the expression in the mass basis to the flavor basis. Unfortunately, this technique may suffer from limitations when the external momenta and the masses within the loops are such that the involved loop integrals suffer from branch cuts [1]. As this will be the case for the observable we will be interested in, we will not consider this technique any further. We refer to the original references for more details.
Alternatively, we will use another technique, the so-called mass insertion approximation (MIA), which is very powerful when computing flavor transitions. The main idea behind this technique is to perform a new computation, independent to that in the physical basis, working directly in the interaction basis. In general, the mass matrices will not be diagonal in this basis and the non-diagonal terms will provide two-point interactions, the so-called mass insertions, which will allow a particle to transform into another one along its propagation. Then, the propagator of a particle is understood as successive insertions of this kind and, therefore, this technique provides a diagrammatic expansion of any observable. Although computing the full series of mass insertions would reproduce the complete result in the physical basis, in general this is not possible to do. Nevertheless, in the case where the mass insertions are small -smaller than the diagonal mass terms -, we can treat these insertions perturbatively and compute the diagrams in the interaction basis to a given order in the MIA expansion. For this reason, the MIA technique is very useful to compute LFV processes, since the strong experimental bounds suggest that any kind of parameter leading to LFV transitions should be small.
The particular model we will choose to perform the MIA computations will be a general type-I seesaw model [69][70][71][72][73], where an arbitrary number of right-handed (RH) neutrinos are added to the SM. This kind of models are very well motivated from the observation of neutrino oscillations [74][75][76][77], and it is well-known that they may lead to sizable LFV transitions, see for instance the recent review [78] and references therein. Furthermore, it has been shown that the MIA works very well in this context [13,79].
The paper is organized as follows. We start by describing the basics for a MIA computation in Section 2, providing two simple examples with the explicit computations in the case of small and large mass insertions. Then, in Section 3 we apply this technique to the case of the LFV Higgs decays in a general seesaw model. We use it to derive an effective vertex in the case of heavy RH neutrinos, which helps to compare the results with current experimental bounds from other observables and to conclude on maximum allowed LFVHD rates. Finally, we conclude in Section 4. Further details on the computation and heavy mass expansions are provided in the Appendices.

Basics for a Mass Insertion Approximation computation
In many models, spontaneous symmetry breaking is behind particle masses. It leads to quadratic terms in the Lagrangian, which implies a mass matrix in the interaction basis. In general, this matrix is non-diagonal and when it is diagonalized the physical basis is derived. This latter basis is in general the chosen basis to perform QFT computations, as it is possible to properly define a loop expansion for a given observable. Nevertheless, when doing this one usually loses track of the parameters in the interaction basis. On the other hand, having analytical expressions for a given observable in terms of the fundamental parameters of the theory, it is possible to extract information about them from the experimental measurements directly.
In order to work with the fundamental parameters in the computation of a given observable, the mass insertion approximation provides a powerful tool. This method is a diagrammatic diagonalization of the mass matrix in the interaction basis: in this approach the diagonal entries are considered as the mass parameters, while the non-diagonal ones are interpreted as two-point interactions (the so-called mass insertions) of the corresponding states. In this context, the propagator of a given state is constructed from the successive mass insertions connecting two different fields and the interaction states are dressed with these consecutive interactions. The exact diagonalization corresponds to a complete resummation of the infinite mass insertions that can occur in the propagation. In general, it is not possible to do this exact resummation, and therefore an approximation is used: as in a Taylor expansion, a dimensionless parameter is defined as the ratio of the non-diagonal mass insertion over the diagonal mass parameter, and its magnitude defines how many mass insertions must be taken into account to achieve a given precision in the expansion. In a general model, the hierarchy between the different mass scales defines different dimensionless parameters.
In this Section, we present two examples as an illustration of the application of this technique. The first one corresponds to a situation in which the non-diagonal terms are smaller that the diagonal ones. Then we show that the first two terms in the MIA expansion reproduce the computation in the physical mass basis to that order. The second example represents the opposite situation: the non-diagonal parameters are larger than the diagonal ones, and we need to perform a complete resummation of the infinite mass insertions.

First example: small mass insertions
As a first example of the MIA application, we consider a toy-model composed of three real scalar fields ρ, Φ 1 and Φ 2 in the interaction basis, with the following the Lagrangian: with an implicit sum over I, J is understood. For simplicity, we assume a real and symmetric squared mass matrix but not aligned in the interaction space: and a cubic interaction between the scalar ρ with the Φ fields that is diagonal in interaction space: The mass matrix M 2 can be diagonalized by an orthogonal matrix O: defining a physical basis φ ± , with physical masses given by Therefore, the Lagrangian in the physical basis is: Let us consider the one-loop contributions to the self-energy of the ρ scalar field: in the physical basis, they come from the Lagrangian of Eq. (6) and correspond to the sunset topology of Fig. 1, with φ ± running in the loop. They are expressed as function of the one-loop integrals defined in App. A: Notice that, in this example, the matrix O is not involved in the computation in the physical basis. This is due to the simplicity of our model, where we assumed λ = λ1 and therefore O is not present in Eq. (6). In a more general case, the expressions in the physical basis are given terms of the physical masses and the rotation matrices.
In order to illustrate how the MIA works, we assume that the non-diagonal entry in the mass matrix M 2 of Eq. (2) is much smaller than the diagonal one. Then the dimensionless parameter m 2 /M 2 1 is defined and it controls the diagrammatic expansion in the interaction basis. Now, considering the interaction fields Φ 1,2 of Eq. (1) running in the loop with an associated mass parameter M and the two-point interaction ∼ m 2 Φ 1 Φ 2 as the mass insertion, the systematic procedure is to add successive mass insertions up to a given order. In Fig. 2, the first two contributions in the MIA are shown: they correspond respectively to the leading order (LO) with no mass insertions, and to the next to leading order (NLO) with two insertions. In this example we cannot close the loop with an odd number of mass insertions in the sunset topology, since the cubic interactions are diagonal. Notice that all these MIA diagrams are at the same loop level, that of the corresponding diagrams in the physical basis, however they are of different order in the MIA expansion. Therefore, in this approach we have an expansion for the self-energy at the one-loop level as where the dots are contributions with more mass insertions and, thus, suppressed by higher powers of the dimensionless parameter m 2 /M 2 . The LO contribution in the MIA has the same type of diagrams than in the physical mass basis, but the interaction states Φ 1,2 are running in the loop now. Then, On the other hand, there are six diagrams contributing to the NLO order in the MIA. Each one has two mass insertion, so they are proportional to m 4 . Explicitly, the one-loop integral for the middle-left diagram in Fig. 2 is given by, Then, the NLO contribution in the MIA is: The analytical comparison between the MIA and the physical basis results is obtained when the physical masses (and the matrix rotations if any) are expressed in terms of the gauge parameters and the physical basis expressions are expanded up to a given order. In this example, since we have used the MIA up to two insertions, we need to expand the expression in the physical basis up to O(m 4 /M 4 ). From Eqs. (5) and (7), the physical basis computation of the scalar ρ self-energy in terms of the gauge interaction parameters is This two-point B 0 one-loop function, participating also in the LO contribution of the MIA in Eq. (9), is given by where γ E is the Euler-Mascheroni constant and µ is the usual scale for dimensional regularization.
We can now expand the expression obtained in the physical basis, Eq. (12), under the assumption of m 2 M 2 . At zero order, m = 0, this equation trivially leads to the LO MIA contribution in Eq. (9). The next terms in the expansion are of order m 4 /M 4 , and lead to the NLO MIA expression in Eq. (11). Similarly, one could check that higher order terms m 8 /M 8 , m 12 /M 12 , . . . will correspond to higher order terms in the MIA expansion with 4, 6, . . . insertions.
We remark again that the simplicity of the present toy-model allows us to compare explicitly the MIA results with the expansion of the physical basis results, due to the analytical diagonalization of the 2 × 2 mass matrix. In a more complex situation, this diagonalization is only performed numerically, and thus the dependence on the interaction parameters is missed. Moreover, the computational effort could be huge for higher dimension mass matrices. In that context, the MIA diagrammatic computation is a powerful tool in order to work with the interaction parameters explicitly. As we said, the MIA results correspond to a perturbative calculation in a dimensionless parameter.

Second example: large mass insertions
Now we analyze a situation in which a non-diagonal entry of the mass matrix is larger than a diagonal one, i.e., the corresponding dimensionless parameter results larger than 1. In that case, Figure 3: Diagrammatic interpretation of the dressed propagators (thick lines) for the same and opposite chiralities as an infinite series of successive mass insertions (black dots) between two massless propagators (thin lines).
an exact resummation of this large mass insertion is needed. In particular, we consider a Dirac spinor ψ = P L ψ + P R ψ = ψ L + ψ R with mass M and momentum p. The corresponding quadratic terms of the free Dirac Lagrangian are where we have a matrix of dimension 2 in chiral space (P L,R = (I ∓ γ 5 )/2). This approach is equivalent to having two massless fermions ψ L and ψ R interacting through the mass insertion M . The corresponding massless propagators are the inverse of the kinetic terms: As before, the dimensionless parameter of the MIA expansion should be the ratio between the mass insertion M and the mass of the fermions. However, the latter is zero in this example, implying an infinitely large expansion parameter. This fact can be solved by defining a dressed propagator that accounts for a resummation of all the insertions of this kind. In the chiral basis, there are four types of propagators depending on the chiralities of the connected fermions (LL, LR, RL and RR), as showed in Fig. 3. Here, the thin lines represent the massless propagators, the blacks dots are the mass insertions and the thick lines correspond to the dressed propagators (after the resummation of the successive two point interactions). The dressed propagators that connect two fermions with the same chirality contain an even number of mass insertions, while the ones connecting two opposite chiralities have an odd number of mass insertions. Explicitly, the propagator connecting two left-handed fermions (LL) corresponds to the geometric series: and the propagator connecting two right-handed fermions (RR) is obtained with the interchange P L ↔ P R : In the same way, the propagator connecting opposite chiralities (LR) comes from the geometric series: and the propagator RL results from the LR after interchanging P L ↔ P R : From the previous relations, we can interpret that the successive non-diagonal two-point interactions dress the propagator providing the corresponding masses. It is important to connect this approach with four types of propagators with the standard one of the Dirac propagator i( / p+M ) p 2 −M 2 : in a generic process with a Dirac propagator, the MIA approach produces four diagrams with the LL, LR, RL and RR propagators. Adding these contributions from Eqs. (16)(17)(18)(19), the complete Dirac propagator is restored. This procedure works in a generic context of two interacting states, as we will see in the next Section for the type-I seesaw model.

MIA in practice: LFV Higgs decays in a general seesaw model
In order to better illustrate the discussion in the previous Section, we apply next the MIA technique to the particular example of LFV Higgs decays in a general type-I seesaw model (GSS), where N right-handed neutrinos are added to the SM. The full computation in the neutrino mass basis 1 was done in [9] -see also [7] -, and the final expressions after correcting some typos can be found in [80]. The MIA technique was applied to the particular case of the inverse seesaw model with three pairs of degenerate sterile fermions [13]. Here, we generalize these results to a GSS and discuss how to apply them to the particular case of low scale seesaw models, recovering the results of [13] in the proper limit. Finally, we apply the constraints from the global fit analysis in [81] to conclude on the maximum allowed H → k m rates.

Model setup for the MIA
We consider a general type-I seesaw model where the SM is extended with N right-handed neutrinos. The corresponding Lagrangian is given by where L is the SM left-handed lepton doublet and Φ = iσ 2 Φ * with Φ the SM Higgs doublet. The fundamental parameters of the model are then the neutrino Yukawa coupling Y ν , which is a 3 × N complex matrix, and the Majorana mass matrix M which is a N × N symmetric matrix that violates lepton number in two units. The C-conjugate is defined as usual as ψ c = Cψ T , where we can choose C = iγ 2 γ 0 . After the EWSB, this Lagrangian leads to a neutrino mass matrix that, in the flavor basis (ν c L , ν R ), reads with the Dirac mass matrix defined as m D = vY ν and v 174 GeV. In the seesaw limit, the non-diagonal entries m D are smaller than the diagonal M and, therefore we can perform a MIA computation, which will be defined as a perturbative expansion in powers of m D /M . Moreover, and despite the fact we will be interested in expressing our results in terms of the EW parameters m D and M , we recall that in this seesaw limit the physical masses of the heavy neutrino will be approximately given by the Majorana mass matrix M , and that the active-sterile mixings in the physical basis will be of the order m D /M . Thus, our MIA computation will be in this sense an expansion in terms of active-sterile mixings.
As discussed in the previous Section, the first step in a MIA computation is choosing the proper basis. Despite the fact that the MIA works in the flavor basis, it is not mandatory to work with the full model in this basis, and it is actually convenient to choose the basis for each sector independently. For the present exercise of computing the LFV Higgs decays at one-loop, we will choose a hybrid basis: we will work with the flavor basis for the neutrino sector, while the rest of the fields will be taken in their mass basis. The latter will apply to the external fields, Higgs boson and charged leptons, as well as to the gauge and Goldstone bosons in the loops, which we will treat in the Feynman-'t Hooft gauge 2 . By doing this, we will obtain a useful expression for the LFVHD rates in terms of the new flavor parameters in Eq. (20) and those already known SM parameters in the physical basis.
Moreover, in the following we will choose the ν R basis such that the Majorana mass matrix is real and diagonal. Notice that we can do this without any loss of generality, and it will only imply that, for models where M is not diagonal at the beginning, we will need to diagonalize it first, as we will do explicitly when discussing low scale seesaw models in Section 3.3. If M is diagonal, then the only lepton flavor changing mass insertion will be the Dirac mass term and, consequently, our MIA computation will be defined as an expansion of successive mass insertions of m D . Nevertheless, in the computation of our observable there is another source of LFV coming from the Y ν of the cubic interactions between the ν R with the L doublet and the H and Goldstone bosons, see the first term of Eq. (20). Since the actual source of LFV is the same in both interaction and mass insertion, it will be convenient to consider the Yukawa coupling as the relevant LFV parameter for the expansion and organize our contributions in powers of the Y ν , as we will see later.
Once we have chosen to work in the flavor basis for the neutrino sector, the Majorana mass matrix M GSS can be understood as the collection of all the relevant mass insertions. From the one side, the already mentioned m D will mix the ν L and ν R fields and we will denote it with a cross, as in Fig. 4. Since every insertion of this kind will introduce a new ν R field, each new insertion is expected to be suppressed by inverse powers of the heavy mass M and, therefore, we can treat the m D mass insertions perturbatively. On the other hand, the mass term M can be understood as a -lepton number violating -mass insertion between the ν R and ν c R fields. This mass insertion, however, is not small and, thus, we need to resum all possible insertions of this kind. Following the discussion in the previous Section 2.2, we can define two kinds of dressed propagators, a lepton number conserving (LNC) one with any even number of insertions, and a lepton number violating Figure 4: Relevant (dressed) propagators and LFV mass insertion for a general seesaw. All the other Feynman rules needed for the computation of H → k¯ m at one-loop can be found in [13].
(LNV) one with an odd number of insertions: where M a is the corresponding element of the diagonal mass matrix M . Notice that we are interested in computing a lepton number conserving process; hence, it will be enough to consider the first of these dressed propagators.
With this setup, the computation is basically the same than that performed in [13], with the dressed propagator in Eq. (22) playing the role of the fat-propagator in [13]. All the other Feynman rules relevant for the computation of the LFVHD are the same, so we refer to [13] for further details and conventions.

The MIA computation and the heavy mass expansions
We are interested in computing the LFV process H(p 1 ) → k (−p 2 )¯ m (p 3 ), whose decay amplitude can be generically decomposed in terms of two form factors F L and F R [9], In order to further simplify our expressions, one could neglect the masses of the charged leptons with respect to the Higgs boson mass. Nevertheless, the form factor F L(R) is proportional to m k(m) , so we cannot fully neglect lepton masses and we need to keep the leading term. Using the fact that charged lepton masses are hierarchical, we work under the hypothesis m m m k . Then, it is enough to consider the F L form factor for the H → k¯ m decay, keeping the leading contribution in m k and neglecting any additional contribution from charged lepton masses 3 . Then, the partial decay width can be written as: Equivalently, in this case where m m m k , F R dominates in the CP-conjugated process H → m¯ k . In the Feynman-'t Hooft gauge and in the neutrino mass basis, these form factors receive contributions from the diagrams in Fig. 5, We remark that the definition of mass basis depends on the perturbation order considered, as loop corrections will generically modified the mass matrix adding non-diagonal terms, which need to be rotated away. We chose to work with the tree-level mass basis and, consequently, we need to include the self-energy corrections to the external legs, last line in Fig. 5. Alternatively, one can work with the one-loop level mass basis, which would rotate away those diagrams. Nevertheless, both techniques are equivalent at the one-loop level.
The full analytical expressions for these form factors can be found in [80]. These expressions are given in terms of the physical neutrino masses and the unitary matrix that relates the flavor and physical basis, so the analytical dependence on the initial parameters in Eq. (20) is lost. Morover, evaluating numerically these expression could be time consuming, in particular when the amount of right-handed neutrinos is large. Therefore, it would be useful to have expressions which are given directly in terms of the flavor parameters in Eq. (20), even if they are approximations.
As we already know the full result in the physical basis, we could in principle apply the Flavor Expansion Theorem proposed in [1]. Nevertheless, this technique requires that the external momenta are smaller than the masses running in the loops, and this is not the case in the decay process H → k¯ m due to the fact that the external momentum of the on-shell Higgs is larger than the ν L and W masses. Therefore, we do not consider this technique for this computation. We perform instead a MIA expansion, since it can be applied even when light masses are running in the loops.
In the MIA, this process can be computed as an expansion of the relevant LFV mass insertions, which in our case are the Dirac mass terms, or equivalently the Yukawa couplings. Thus, each of the diagrams in Fig. 5 can be computed in the MIA as, The fact of having only even powers of Y ν is related to the right-handed neutrinos, whose presence is needed to induce the LFV transition. Since they only interact via the Yukawa couplings, each (dressed) RH propagator will introduce two Y ν , one at each edge of the propagator. This means that the LO terms O(Y 2 ν ) will come from MIA diagrams with only one RH propagator, the NLO corrections O(Y 4 ν ) from diagrams with two, and so on. Being the Yukawa coupling perturbative, it should be enough to compute the first contributions to this expansion. Moreover, the addition of RH propagators will introduce inverse powers of M and ensure the convergence of Eq. (27) as a perturbative expansion in terms of O(m 2 D /M 2 ). We collect the MIA results in App. A, as well as the relevant Feynman diagrams entering the computation. We include the complete O(Y 2 ν ) terms, which give a good description of the observable when the Yukawa couplings are small. Nevertheless, for large Yukawa couplings, we need to compute also some of the O(Y 4 ν ) terms [13], which are not as suppressed as we may naively expect from the above discussion. These dominant O(Y 4 ν ) terms are also given in App. A. In order to better understand this point, it is useful to analyze our results when the Majorana scale is heavier than the electroweak scale. D /M 4 ) and will be, therefore, negligible. Nevertheless, it turns out that some diagrams lead to O(Y 4 ν ) terms that are suppressed only by two inverse powers of M , and thus they are important to take into account for a good MIA prediction. Indeed, for very large couplings O(Y ν ) ∼ 1, these terms may become the dominant ones. In the particular process we are considering, they come from diagrams of type (1), (8) and (10)  It is interesting to discuss a bit more the presence of this kind of unsuppressed O(Y 4 ν ) terms in different LFV observables induced from heavy neutrinos. Besides in the LFV Higgs decays, similar contributions were found in the context of LFV Z boson decays [79], which at the same time suggests that they are also present in LFV 3-body decays, such as k → 3 m , due to the strong correlation between these two latter observables [82,83]. On the contrary, these terms are not present in the case of LFV radiative decays k → m γ [10]. The difference could be tracked to the fact that neutrinos do not couple to the photon, but they do couple to the H and Z bosons, leading for example to diagram (1) in Fig. 5 for the LFVHD, which is not present in the k → m γ process 4 . The fact that the radiative decays are different for these other LFV processes is very interesting, as the former are usually the most constraining LFV processes, however this may not be true at very large Yukawa couplings due to these additional terms. Now, collecting all the relevant terms of O(M −2 ), we arrive to the following effective vertex for the LFV H → k¯ m decay, where we have defined: For the physical values of m H = 125 GeV and m W = 80.4 GeV we have r(m 2 W /m 2 H ) ∼ 0.31. We recall again that this expression is valid under the assumption of heavy Majorana masses M v and in the seesaw limit m D M , since we only kept the O(M −2 ) terms and we performed the computation at NLO in the MIA expansion. As we will see later, these terms are enough to reproduce the full results to a good accuracy in the parameter space that is still allowed by current constraints.
In Eq. (28), we have to sum over the indices a and b, which run over the RH neutrinos. In general, all of them will contribute and the indices will run from 1 to N . Nevertheless, in some interesting cases some of the RH neutrinos might be very heavy and they will completely decouple from the observable. Since the contribution to any very heavy neutrino to Eq. (28) is negligible, decoupling a RH neutrino is indeed equivalent to restricting the range of a and b to those (nondecoupled) right-handed neutrinos, which are still light enough to contribute.
Another interesting limit corresponds to the case with complete degenerate RH neutrinos, i.e., M 1 = ... = M N ≡ M . In this case, the effective vertex corresponds to Notice that, even if have focused on the H → k¯ m channel, the effective vertex V eff H m¯ k of the CP-conjugated process H → m¯ k can be easily obtained by conjugating Eq. (28): which is equivalent to interchanging the flavor index of the charged leptons k and m in the Yukawa couplings. Finally, the branching ratio for the process V eff H k¯ m can be computed by just plugging the corresponding effective vertex in Eq. (25), where Γ tot H is the total width of the Higgs boson. In order to illustrate the accuracy of the effective vertex, we show in Fig. 6 the predictions for H → τē computed using the full expressions in the mass basis (solid lines) as well as using the approximated expression in Eq. (32). Moreover, we quantify the agreement by means of the ratio R, defined as the approximated prediction over the full one. We choose here a particular realization of the seesaw model that we will introduce in Eq. (50), although we found similar results for other examples. The differences between the two panels are the heavy neutrino masses, chosen to be degenerated in the left and hierarchical in the right. The overall conclusion from this figure is that the effective vertex in Eq. (28) works very well in both cases, as long as the condition m D M is fulfilled.

Connection to low scale seesaw models
As a particular but interesting application of the effective vertex in Eq. (28), we apply it to the so-called low scale seesaw models (LSS), such as the inverse [85][86][87] or linear [88] seesaw models. These models are of great phenomenological interest since they can introduce relatively light heavy neutrinos with large Yukawa couplings, and still accommodate naturally the observed light neutrino masses. Moreover, this will also allow us to compare with existing results in the inverse seesaw model [13].
The common feature in all these models is the imposition of an approximate conservation of lepton number [89], which is only violated by some small parameter related to light neutrino masses. The difference between each particular low scale seesaw realization is precisely the nature or origin of this small LNV parameter. However all of them share the same lepton number conserving limit. In our case of study, the LFVHD are LNC processes, and thus the small LNV parameters will not play any important role and we can be neglected. This means that we can expect to have the same LFVHD rates for all these low scale seesaw models.
The LNC low scale seesaws could be realized by adding n ν pairs of new fermionic singlets with opposite lepton number, which we denote ν a R and ν a S respectively, with a = 1, . . . , n ν . For the purpose of this discussion, we are just interested in the neutrino mass matrix 5 , which in the (ν c L , ν R , ν S ) basis reads as where Y LSS and M R are respectively 3 × n ν and n ν × n ν matrices, and the zeros have the proper dimensions so the total M LSS matrix is a (3 + 2n ν ) symmetric matrix. In order to apply our results from the previous Section, the first step is to rotate the heavy neutrino sector to its diagonal form, 5 For more details on LSS, see for instance [90].
where now the dimensions of the new Dirac m D and Majorana mass M diag matrices are 3 × 2n ν and 2n ν × 2n ν , respectively, and U is a unitary matrix rotating the neutrino sector. Without any loss of generality, we can assume that M R is already diagonal and, therefore, the diagonalization becomes trivial, where the unitary matrix V just contains rotations of π/4 and i factors to make the entries of M diag possitive. In general, this unitary rotation is defined by four blocks of n ν × n ν Finally, the new Dirac mass matrix is given by We have now all the pieces needed to compute the H → k¯ m process, we just need to plug the Dirac matrix m D = vY ν of Eq. (37) and the M diag of Eq. (35), in the effective vertex of Eq. (28). For instance, we can consider the same setup as in [13], where all the entries of M R are degenerate. In that case, we can use the Eq. (30), with resulting in the effective vertex which is agreement with the result in [13], obtained for the particular case of the inverse seesaw model. Notice that, even that this equation seems to be the same as Eq. (30), it is now expressed in terms of the parameters of the low scale seesaw parameters in Eq. (33), whose physical interpretation is different from the parameters in Eq. (21).

Numerical analysis of the LFV Higgs decays
We conclude this Section by applying the derived effective vertex to study how large the LFVHD rates could be in a GSS model after having considered possible constraints from other observables. For that purpose, we will follow the global fit analysis done in [81], where two different scenarios were considered: a model with only 3 heavy RH neutrinos (3N-SS), and a general seesaw with an arbitrary number of them, as in Eq. (20). In both case, the Authors obtained upper bounds on the η matrix, a small Hermitian matrix encoding the deviations from unitarity in the light neutrino mixing. In our case, this matrix can be expressed as 6 , It is interesting to analyze first the 3N-SS case, as it is simpler. If we assume again the LNC limit, then it can be implemented by Notice that in this LNC limit light neutrinos are strictly massless, but they can be accommodated by introducing small LNV parameters in these matrices [81]. Nevertheless, since the LFVHD do not violate lepton number, these small LNV parameters will not be relevant for our observable and, therefore, we neglect them in the following. In this scenario, the effective vertex in Eq. (28) becomes, where we have used that Alternatively, it can be also written in terms of η as, (44) Notice that the observable vanishes in the Λ → ∞ limit, as it is manifest in Eq. (43). This decoupling behavior is hidden when we express it in terms of η, but this form is useful to conclude on maximum allowed LFVHD rates in this model. Due to the O(Y 4 ν ) terms, the maximum rates will be obtained at the largest value of Λ that allows to saturate Eq. (41) without spoiling the perturbativity of the Yukawa couplings. Assuming a perturbativity bound of |(Y ν Y † ν ) km | < 4π, a rough estimation points to Λ ≈ 10 TeV and consequently, where we have defined BR(H → k m ) ≡ BR(H → k¯ m )+BR(H → m¯ k ). The differences between the three channels have a double origin. From the one side, the fact that these decays are proportional to charged lepton masses suppressed the H → µe with respect to the other two by a factor of m 2 µ /m 2 τ . On the other hand, current bounds in Eq. (41) are a bit stronger for η τ µ , and much more stringent for η µe due to the strong constraints coming from µ → eγ [91]. This double suppression makes of H → µe an even more challenging observable for experiments.
The effective vertex in Eq. (28) has the advantage of being more general than previously computed ones. In particular, it allows to explore scenarios where the heavy neutrinos have different Majorana masses, as we show in Fig. 7. In order not to generate too large masses for light neutrinos, we consider an scenario such as the one in Eq. (33), but with only two pairs of Majorana neutrinos (n ν = 2) contributing to the Higgs decays, assuming that any possible additional neutrino has decoupled from this observable. Moreover, and in order to avoid the strong bounds in the µ-e sector, we have considered a simplified case for the Yukawa couplings where none of the heavy neutrinos couple to muons (left panel) or to electrons (right panel). The rest of the entries of the Yukawa coupling matrix are set to one, again for simplicity. This figure has been done using Eq. (28), although we have checked that the full computation leads to the same results in the relevant area allowed by the constraints of Eq. (41).
In this Fig. 7 we find again the decoupling behavior with each of the individual heavy masses, as we discussed before. We can also see that the predictions for H → τ e and H → τ µ are the same in both panels, as we have assumed a simplified scenario with equal size Yukawa couplings to the relevant flavors. Nevertheless, the bounds of Eq. (41) are stronger in the muon sector than in the electron one, which again implies that the allowed rates are larger for H → τ e than for H → τ µ. Finally, we notice that the results are symmetric with respect to M R 1 and M R 2 , although this is again a consequence of our simplified hypothesis of equal Yukawa couplings, and in general this will not be the case.
In a more general scenario, the η matrix will constrain a combination of all the heavy neutrino mixing to the active state, each of these contributions being of order ξ a = m a D /M a , where ξ = m D M −1 is the usual seesaw parameter. Generalizing the description of [13] to the case of . Shadowed areas are disallowed by global fits constraints (purple) and the non-perturbative regime for the Yukawa couplings (yellow), which we define as |(Y ν Y † ν ) km | > 4π for any element.
having N RH neutrinos, we can think of this parameter as made of three N -vectors: which leads to This implies that the diagonal entries of the η matrix constrain the modulus of these N -vectors, while the diagonal ones set upper bounds on the complex scalar products between them. Moreover, this geometrical picture is also useful to find solutions that accommodate the current bounds given in Eq. (41), as we only need to set the modulus and angles between these vectors, and to explore the implications on our LFV observable. Let us consider a final example to illustrate this latter point. In order to avoid issues generating too large contributions to light neutrino masses, an elegant solution is to consider a low scale seesaw realization of the model, as we introduced in Eq. (33). If we assume, for example, that there are two pairs of sterile fermions (n ν = 2) contributing to our observable and that they have the same mass, M R 1 = M R 2 ≡ M R , an interesting value for the Yukawa matrix is the following: This example is useful as it leads to a η matrix with a very similar pattern than that in Eq. (41). Then, depending on the value of the global strength factor Y 0 and the mass M R , we can define which part of the parameter space is allowed. We show in Fig. 8 the predictions for the three LFVHD channels in this particular example, although we expect similar results for other models with more RH neutrinos as long as they lead to the similar η matrices. The first thing we see from this figure is that the effective vertex (dashed lines) matches very well the full prediction (solid lines), as we already saw before. Nevertheless, the simple expression of Eq. (28) allows us to easily understand the dependence on the different parameters of the model. Finally, we can also use this figure to deduce how large LFVHD rates can we expect after having considered the bounds in Eq. (41) as well as perturbativity bounds. As we discussed above, the largest possible values are obtained for heavy masses at around 10 TeV, when the bounds from the global fit analysis and perturbativity cross. Unfortunately, the branching ratios in the allowed white area are too small and far from current experimental sensitivities.

Conclusions
In this work, we have discussed the importance of having expressions for lepton flavor violating transitions which are expressed directly in terms of the fundamental parameters of the model. These expressions are helpful to better understand the observable, as well as to compare it with experimental observations in order to constrain the interaction parameter space. In this context, the mass insertion approximation technique is a powerful tool, which we have reviewed with two simple examples.
We have then studied the LFV Higgs decays in the context of a general type-I seesaw model with an arbitrary number of right-handed neutrinos. We applied the MIA technique to this model, which allowed us to derived an effective vertex for H k m after integrating out the heavy right-handed neutrinos. This analytical expression is useful to understand the behavior of the observables with the fundamental parameters of the neutrino sector, i.e., the Yukawa couplings and the heavy Majorana masses. Moreover, it also provides an alternative way to quickly evaluate these observables to a very good approximation, without the need of long numerical evaluations of the full result in the physical basis.
Finally, we have made the connection to the phenomenologically interesting case of low scale seesaw models. After explicitly checking that we recover existing results for the inverse seesaw case, we have evaluated the LFVHD rates taking into account current bounds from global analysis, as well as perturbativity bounds for the Yukawa couplings. Unfortunately, the predicted rates for the LFVHD in the allowed area are far from current experimental sensitivities and, consequently, they do not provide a competitive way of probing the existence of these heavy Majorana neutrinos. project FPA2016-78645-P, and by the Spanish MINECOs "Centro de Excelencia Severo Ochoa" Programme under grant SEV-2016-0597.

A MIA Form factors
The analytical computation of the one-loop effective vertex in the MIA approach are written in terms of the one-loop functions for D dimensions in the Passarino-Veltman notation [92], following these definitions and conventions: As explained in the text, we consider the hierarchy between charged lepton masses m m m k and thus we only consider the left-handed form factor F L in Eq. (24). In order to simplify our expressions we neglect the contributions of O(m 3 ) in the form factors because they are very suppressed and we recovered the full predictions without them. We present the analytical results for the left-handed form factors computed within the MIA to one-loop order and considering the leading order contributions, O(Y 2 ν ), and the next to leading order corrections, O(Y 4 ν ), as explained in the text. We collect in Figs. 9-10, the relevant diagrams in the Feynman-'t Hooft gauge. We follow the notation for the diagrams given in [13]. Missing diagrams correspond to suppressed contributions proportional to higher lepton masses powers.
For shortness, we introduce the factor f k = 1 32π 2 m k m W , so that the results of the O(Y 2 ν ) dominant contributions diagram by diagram are: On the other hand, the results of the dominant O(Y 4 ν ) contributions are: In the above expressions, an implicit sum over i = 1 . . . 3 and a, b = 1 . . . N has to be understood, and the arguments of the one-loop integrals are: Although we have derived these expressions assuming the approximations in Eq. (57), we have checked that they are valid as long as p 2 2 = 0 and in the branch p 2 1 < 4m 2 W . We checked that the above expansions are in very good agreement with the numerical results from LoopTools [93]. We also checked that the one-loop functions with two different Majorana masses M a and M b reproduce the corresponding expressions in the degenerate limit M a = M b = M R given in [13].
An important comment is in order: an usual approximation to do those type of integrals is the zero external momentum (p 2 1 = p 2 2 = p 2 3 = 0). If we apply this approximation, the Higgs boson mass coming from Eq. (57) Although this equation corresponds to the third line of Eq. (58) in the limit p 2 1 → 0, it is not enough to reproduce the full result for the LFVHD. The reason is the fact that the contributions of the on-shell Higgs boson momentum are competitive with the W boson mass in the loop.
Finally, the LO and NLO contributions of the left-handed form factors (grouped in convenient pairs) in the large Majorana mass limit are: