Majorana Neutrino Masses in the RGEs for Lepton Flavour Violation

We suppose that the observed neutrino masses can be parametrised by a lepton number violating dimension-five operator, and calculate the mixing of double insertions of this operator into lepton flavour changing dimension-six operators of the standard model effective theory. This allows to predict the log-enhanced, but $m_\nu^2$-suppressed lepton flavour violation that is generic to high-scale Majorana neutrino mass models. We also consider the Two Higgs Doublet Model, where the second Higgs allows the construction of three additional dimension-five operators, and evaluate the corresponding anomalous dimensions. The sensitivity of current searches for lepton flavour violation to these additional Wilson coefficients is then examined.


Introduction
Neutrinos are elusive and enigmatic particles: uncoloured, uncharged, and very light. Nonetheless, their observed masses and mixing angles [1] imply that Lepton Flavour Violation (LFV) must occur, where we define LFV as flavour-changing contact interactions of charged leptons (for a review, see e.g. [2]). Since these do not occur in the Standard Model (SM), LFV is considered to be "New Physics", and searched for in a wide variety of experiments [1,3,4,5,6,7,8,9,10,11,12]. Neutrinos could also induce another kind of New Physics: if their small masses are "Majorana", they are Lepton Number Violating (LNV), and could for instance mediate neutrinoless double-β-decay [13]. Below the weak scale, such masses appear as renormalisable terms in the Lagrangian, but in the full SU(2) gauge invariant Standard Model, they arise as a non-renormalisable, dimension-five operator.
In this paper, we will assume that neutrino masses are Majorana, and that the scale Λ of New Physics in the lepton sector is large. We focus on the theory at scales above m W but below Λ, where it can be described in the framework of the Standard Model Effective Field Theory 1 (SMEFT). The neutrino masses can be parametrised by operators of dimension five, and LFV is parametrised by operators of dimension-six. Our aim is to obtain the log-enhanced loop contributions of two LNV operators to LFV processes. These can be calculated via renormalisation group equations (RGEs), and in particular we aim to calculate the anomalous dimensions that mix two dimension-five operators into a dimension-six operator. The renormalisation group running of the dimension-five operators has been extensively studied in the literature [16,17,18], and the mixing of the dimension-six operators among themselves have been evaluated at one-loop [19] in the "Warsaw"basis [20] of SMEFT operators. The mixing of two dimension-five operators into dimension-six operators was calculated in [21], using the Buchmuller-Wyler [22] basis at dimension-six. We perform this calculation using the "Warsaw"-basis, and our results appear to disagree with [21].
The mixing of neutrino masses into LFV amplitudes is O(m ν /m W ) 2 ln(Λ/m W ), so negligibly small, but completes the anomalous dimensions required to perform a one-loop renormalisation-group anlysis of the SMEFT at dimension-six. In addition, we explore an extention of the SMEFT with two Higgs doublets [23], where the second Higgs doublet lives at a scale m 22 between m W and significantly below the lepton number/flavour-changing scale Λ, and we impose that LFV at the weak scale is still described by the dimensionsix operators of the SMEFT. In this scenario, there are four LNV dimension-five operators above m 22 , but only one combination of coefficients contributes to neutrino masses. We calculate the mixing of these LNV operators into the LFV operators of the SMEFT, and estimate the sensitivity of current LFV experiments to their coefficients.
The paper is organised as follows. In Section 2 we introduce the notation of our standard model and two-Higgs doublet model calculation. The main results are presented in Section 3, where we discuss the general structure of our calculation and give the relevant counterterms, anomalous dimensions and renormalisation group equations. Section 4 discusses the phenomenological implications of both results before we conclude.
We provide the relevant Feynman rules, further details of the calculation (including a careful treatment of the flavour structures), and the renormalisation group in the Appendices A -C. The LFV operators of the SMEFT are recalled in Appendix D, and Appendix E gives the current experimental constraints on some LFV coefficients of the SMEFT at the weak scale. Appendix F provides a comparison with the previous calculation of [21] and Appendix G presents the lepton conserving contributions to the anomalous dimensions.

Notation and Review
The SM Lagrangian for leptons can be written as where Greek letters represent lepton generation indices in the charged-lepton mass eigenstate basis, [Y e ] is the diagonal charged-lepton Yukawa matrix, ℓ is a doublet of left-handed leptons, and e is a right-handed charged-lepton singlet. The explicit form of the lepton and Higgs doublets is which have hypercharge y ℓ = −1/2 and y H = 1/2 respectively. The covariant derivative for a lepton doublet is where τ a are the Pauli matrices. This sign convention for the covariant derivative agrees with [19]. Heavy New Physics can be parameterised by adding non-renormalisable operators to the SM Lagrangian that respect the SM gauge symmetries [22]. There is only a single operator at dimension-five in the SM, which is the Lepton Number Violating "Weinberg" operator [24] which is responsible for Majorana masses of left-handed neutrinos. The resulting effective Lagrangian at dimension-five is (ℓ α εH * )(ℓ c β εH * ) + C αβ * 5 2Λ (ℓ c β εH)(ℓ α εH) , (2.4) where ε is the totally antisymmetric rank-2 Levi-Civita symbol with ε 12 = +1, all implicit SU(2) indices inside brackets are contracted, and the charge conjugation acts on the SU(2) component ℓ i of the lepton doublet as ℓ i c = Cℓ i T . The charge conjugation matrix C fulfils the properties of the charge-conjugation matrix used in [25] 2 . The coefficient C αβ 5 is symmetric under the interchange of the generation indices α, β, the New Physics scale Λ is assumed ≫ m W , and the second term is the hermitian conjugate of the first.
In the broken theory, with H 0 = v + (h/ √ 2), v ≃ m t , this gives a Majorana neutrino mass matrix In the charged lepton mass eigenstate basis, this mass matrix is diagonalised by the PMNS matrix [m ν ] αβ = U αi m νi U βi . At dimension-six, we will be interested in SM-gauge invariant operators that violate lepton flavour; a complete list is given in Appendix D. Following the conventions of [20,19], they are added to the Lagrangian as: where X is an operator label and ζ represents all required generation indices which are summed over all generations. Of particular interest are the operators that can be generated at one-loop with two insertions of dimension-five operators, as illustrated in figure 1. With SM particle content, these operators involve two Higgs doublets and two lepton doublets, four lepton doublets, or three Higgs doublets and leptons of both chiralities. In the "Warsaw" basis [20], the possibilities at dimension-six are: where we normalise the "Hermitian" operators with a factor of 1/2 (see appendix D for a discussion) in order to agree with [20,19], and The choice of operator basis implies a choice of operators that vanish by the Equations of Motion (EOMs). For example i D / ℓ α − [Y e ] ασ He σ = 0 implies that the following operators In practice, if the coefficients C βα HDℓ (1) and C βα HDℓ(3) of these structures are present, they are equivalent to C βσ eH = C βα HDℓ (1) [Y e ] ασ + C βα HDℓ (3) [Y e ] ασ (and the hermitian conjugate relation).

In the case of the 2HDM
In this section, we consider the addition of a second Higgs doublet H 2 to the SM, of the same hypercharge as the SM Higgs (which we relabel H 1 ). The LFV induced by double-insertions of dimension-five operators could be more significant in this model, because there are several dimension-five operators, so neutrino masses cannot constrain them all. However, a complete analysis of LFV in the 2HDM would require extending the operator basis at dimension-six and calculating the additional terms in the RGEs, which is beyond the scope of this work. So for simplicity, we make three restrictions: 1. First, we consider only the dimension-six LFV operators of the SMEFT. This is the appropriate set of dimension-six operators just above m W , provided that H 2 has no vev, and that the mass m 22 of the additional Higgses is sufficiently high: m 2 W ≪ m 2 22 ≪ Λ 2 . In our phenomenological analysis we extend this range to the scenario m 2 W m 2 22 ≪ Λ 2 , by considering a Higgs potential where the additional Higgses are not directly observable at the LHC, and where the Yukawa couplings of H 2 are vanishing. Such a scenario would for example be realised in the inert two Higgs doublet model [26,27,28,29] and setting the scale m 22 close to the electroweak scale will not require the consideration of additional renormalisation group effects in the SMEFT.
2. Second, we suppose that at the high scale Λ no dimension-six LFV operators are generated. This is unrealistic, but allows us to focus on the LFV generated by double-insertions of the dimension-five operators.
3. Third, we suppose there is no LFV in the renormalisable couplings of the 2HDM (in particular, in the lepton Yukawas), so that when matching the 2HDM + dimension-five operators onto the SMEFT at the intermediate scale m 22 , no additional LFV operators are generated.
Consider first the renormalisable Lagrangian. The Yukawa couplings can be written [30]: where the flavour indices are implicit, and the basis in (H 1 , H 2 ) space is taken to be the "Higgs basis" where H 2 = 0. We suppose that [Y (1) ] and [Y (2) ] are simultaneously diagonalisable on their lepton flavour indices. The second Yukawa coupling changes the Equations of Motion for the leptons, so the 2HDM version of the equation-of-motion vanishing operators (given in eqn (2.9) for the single Higgs model) should be modified. As a result, the operators O HDℓ (1) and O HDℓ (3) should not be replaced only by the SMEFT operator O eH , as given in eqns (2.10), but also by an operator with an external H 2 leg. However, since we neglect dimension-six operators with external H 2 , we use the relations (2.9) and (2.10) also in the 2HDM case.
In this "Higgs" basis, the most general Higgs potential is In order to decouple the additional Higgses, we can, for instance, set m 2 12 = 0 and assume m 2 22 ≫ m 2 W , or leave m 2 22 free, and impose m 2 12 = λ 6 = λ 7 = [Y (2) ] = 0. At dimension-five in the 2HDM, there are four operators [16]: , but both terms are retained here because they are convenient in our Feynman rule conventions 3 .
Tree-level LFV is often avoided in the 2HDM by imposing a Z 2 symmetry on the renormalisable Lagrangian: if under the Z 2 transformation, H 1 → H 1 and H 2 → −H 2 , then [Y 2 ], λ 6 and λ 7 are forbidden. We will later discuss this case, but do not impose the Z 2 symmetry from the beginning, because it also forbids the C 21 , C 12 and C A coefficients at dimension-five.

Diagrams and Divergences
Diagrams with two insertions of the dimension-five operators are illustrated in figures 1 and 2. We focus on the lepton flavour violating diagrams of figure 1, and discuss the four-Higgs operators generated by figure 2 in Appendix G, because four-Higgs interactions are flavour conserving and arise in the SM. 3 The operator O 21 can also be written as 2(ℓ β ǫH * 1 )(ℓ c α ǫH * 2 ) +(ℓ β ǫℓ c α )(H * 1 ǫH * 2 ) using the identity (A.9), as done in the first reference of [16]. The Feynman rules arising from the (tree-level) Lagrangian of equations (2.1, 2.4, 2.6) are given in Appendix A. We use them to evaluate, using dimensional regularisation in 4 − 2ǫ dimensions in MS, the coefficient of the 1/ǫ divergence of each diagram of figure 1. These coefficients can be expressed as a sum of numerical factors multiplying the Feynman rules for the dimension-six operators of equations (2.7) and (2.10) (these Feynman rules are given in Appendix A), and then the EOMs are used to transform the operators of eqn (2.10) to O eH and O † eH . The required counterterm ∆C O for each of the dimension-six operators given in eqn (2.7) can be identified as (−1)× the numerical factor that multiplies its Feynman rule. This counterterm is added in the Lagrangian to the operator coefficient C O , resulting in a "bare" coefficient C O,bare = µ 2ǫ (C O + ∆C O ) that should be independent of the MS renormalisation scale µ. Note that the factor µ 2ǫ is chosen such that bare Lagrangian remains d-dimensional.
A more complete and rigorous presentation will be required in the next section, in order to derive the RGEs, so let us replace counterterms by Z factors in order to minimise notation and introduce the necessary factors of µ 2ǫ to obtain the correct dimensions. More details of the formalism and calculations are given in Appendix C.
We allow for multiple operators at both dimension-six and -five, and align the dimension-six coefficients in a row vectorC, and the dimension-five coefficients in a row vector C. Then the bare coefficients can be written where matrices wearing a hat act on the space of dimension-six coefficients, and matrices in square brackets act in the dimension-five space, soẐ represents the renormalisation of dimension-six coefficients amongst themselves, and [Z] represents the renormalisation of dimension-five coefficients. The quantity [Z] renormalises insertions of two dimension-five operators; [Z] ij k is a vector in the dimension-six space with index k, and a matrix in the dimension-five space with indices i, j. In the single Higgs model, i, j correspond to the flavour indices of the Weinberg operator, e.g. i = αβ, j = ρσ. The index k corresponds to the operator labels and flavour indices of dimension-six operators. The counterterms that renormalise the diagrams of figure 1 are then components of the vector C[Z] C † . All terms in the above expressions assume an implicit sum over flavour indices; the explicit flavour dependence is presented in Appendix B.
The first diagram of figure 1 has two Higgs and two doublet-lepton legs and so must be renormalised by the operators O Hℓ (3) and O Hℓ (1) , and the structures S HDℓ(1) and S HDℓ (3) . Since these all involve a derivative, the diagram is calculated for finite external momenta. The counterterms that we obtain from this diagram differ from those given in [21]; as discussed in Appendix F, it appears that the authors of [21] dropped one of the terms multiplying the 1/ǫ divergence. We check our result by attaching an external B µ or W a µ boson, in all possible ways, to the first diagram of figure 1, and verify that our counterterms also cancel the divergences of the 2-Higgs-2-lepton-gauge boson vertices generated by two insertations of the Weinberg operator (this is outlined in Appendix B.3). This diagram can be renormalised using the following counterterms: 3) where the last two counterterms represent divergences proportional to the structures S HDℓ(1) and S HDℓ (3) , which contribute to the renormalisation of C eH through the linear combination given in eqn (2.10). The middle diagram of figure 1 contributes to O βα eH , and the divergence it induces can be removed by the counterterm (16π 2 ǫ) −1 [C 5 C * 5 Y ] βα (where the flavour index order is doublet-singlet). Including also the counterterms for S βα HDℓ(1) and S βα HDℓ(3) (eqns (3.4,3.5)) gives Since the structures S HDℓ(3) and S HDℓ(1) are hermitian, they contribute to the renormalisation of both O eH and O † eH (see eqn(2.10)). Only the contribution to O eH is included in (3.6), because the hermitian conjugate in (2.6) generates a counterterm proportional to O † eH that absorbs the divergence of the the "conjugate" process of figure 1.
The third diagram of figure 1 contributes to the four-lepton operator O ραβσ ℓℓ , and the divergence it induces can be removed by the counterterm (3.7)

The 2HDM
In the 2HDM, we consider diagrams analogous to figure 1, but with insertions of any of the dimension-five operators given in eqn (2.13). The external Higgs lines are required to be H 1 , but the internal Higgs lines can be either doublet. The counterterms required to cancel double-insertions of the O 5 operator, discussed in the previous section, also arise in the 2HDM. In this section, we only list the additional contributions to the counterterms. We start again with the first diagram of figure 1, with O 21 or O A at the vertices. Since by construction, the Feynman rule for O 21 is identical to the rule for O 5 , double-insertions of O 21 require the same counterterms as given in eqns (3.2) to (3.5), but with C 5 , C * 5 replaced by C 21 , C * 21 . Double insertions of the antisymmtric operator O A require the counterterms: Finally, O A at one vertex and O 21 at the other require the contributions to the counterterms: It is straightforward to check, using respectively the antisymmetry and symmetry of C A and C 21 on flavour indices, that the combination Including also the additional counterterms for O βα HDℓ (1) and O βα HDℓ (3) in the 2HDM gives (3.13)

The Renormalisation Group Equations
The contribution of dimension-five operators to the Renormalisation Group Equations of dimension-six operators, due to double insertions, can be obtained following the discussion of Herrlich and Nierste [31]. The derivation is presented in Appendix C. Here we schematically outline the result. The bare Lagrangian coefficients are defined at one loop as in eqn (3.1), where the counterterm for one operator can depend on the coefficients of other operators. Recall that the bare coefficients are independent of the dimensionful parameter µ, and that the renormalised Cs are dimensionless.
where [γ] denotes the 4-dimensional anomalous dimension matrix, and we (unconventionally) 4 factor the 16π 2 out of the anomalous dimension matrices. While the −2ǫ term does not contribute in d = 4 dimensions to the mixing of the dimension-five operators, it plays an essential role in the renormalisation group equations of the dimension-six operators.
For the dimension-six coefficients, it is straightforward to obtain from eqn (3.1): whereγ is the one-loop anomalous dimension matrix for dimension-six operators [19] and [γ] = 2(16π 2 )ǫ[Z] is the anomalous dimension tensor. We give below the anomalous dimensions describing the one-loop mixing of double-insertions of dimensionfive operators into LFV dimension-six operators, in the 2HDM. The single Higgs model can be easily retrieved by setting C 21 = C A = C 22 = 0 in the equations below. The anomalous dimension tensor mixing a pair of dimension-five operators into a dimension-six operator is neccessarily a three-index object; below we sum over the two dimension-five indices, and give these summed components of the tensor as elements of a vector in the dimension-six operator space. These anomalous dimensions parametrise the mixing of figure 1 in the 2HDM (recall that a factor 1/16π 2 is scaled out of our anomalous dimensions): where the operator label and flavour indices on the left-hand-side refer to the dimension-six operator (the dimension-five indices are summed).
In the next section, we will need the RGEs for dimension-five operators. Recall that in the single Higgs model, [γ] is in principle a 9×9 matrix (or 6×6, if one uses the symmetry of C αβ 5 ), mixing the elements of C 5 among themselves. However, in the basis where the charged leptons are diagonal, [γ] is diagonal, and the anomalous dimension for the coefficient C αβ 5 of the Weinberg operator is [16]: where the Higgs self-interaction in the SM Lagrangian is λ 4 (H † H) 2 , and [Y f ] are the fermion Yukawa matrices.

Phenomenology
In order to solve the RGEs, it is convenient to define t = 1 16π 2 ln µ mW , in which case the one-loop RGEs for dimension-five and -six operator coefficients can be written as These are among the most familiar of differential equations, whose solutions have the form where 16π 2 t f = ln Λ mW . In these solutions, the anomalous dimension matrices were approximated as constant; this is not a good approximation, because the anomalous dimensions depend on running coupling constants, in particular the Yukawa couplings can evolve significantly above m W .
A simple solution to eqn (4.3) can be obtained by expanding the exponentials under the integral, as in eqn (4.2):C

The single Higgs model
In the SM case where there is only one Higgs doublet, there is only the Weinberg operator at dimension-five: a symmetric 3 × 3 matrix, whose entries are determined by neutrino masses and mixing angles (in the mass basis of charged leptons). We now want to estimate the contribution of double-insertions of this dimension-five operator to lepton-flavour violating processes. We neglect the "Majorana phases", suppose that the lightest neutrino mass is negligible, and and neglect the lepton Yukawas in the RGEs. Then the RG running of C αβ 5 between m W and Λ can be approximated as a rescaling, with γ ≈ λ − 3g 2 + 6y 2 t ≈ 3.5: For Λ ≤ 10 16 GeV, the log is ≤ 32.
We can now estimate the contribution of the neutrino mass operator to lepton flavour violating processes from eqn (4.4). We neglectC(Λ) and find that the contribution is 1 16π 2 ln Λ mW × the coefficients of eqns (3.17) to (3.20), that is, of orderC As expected, this is negligibly small, because C 2

The two Higgs doublet model
Experimental Neutrino data constrain the dimension-five operator in the one Higgs doublet model, so the lepton flavour violating effects estimated in eqn (4.6) are suppressed by the smallness of the neutrino masses. The situation changes in an extended Higgs sector, where more than one dimension-five operator is present. The operator O A cannot contribute to neutrino masses as it is anti-symmetric in flavour space and is hence unconstrained. In addition, the neutrino mass contribution of operators O 21 and O 22 is suppressed if the vacuum expectation value of the second Higgs doublet is small. Renormalisation group effects [16,17,18] will in general mix all operators, which could lift these suppression mechanisms at loop level. However the mixing factorises in the limit where λ 6 , λ 7 and Y (2) tend to zero: then the operators O 21 and O A will not mix into O 5 and O 22 and are hence not constrained by the observed neutrino masses. Furthermore, the mixing of O 22 into O 5 vanishes in the limit where in addition λ 5 tends to zero (see [32] for a symmetry argument). In the following we will study the sensitivity of lepton-flavour violating decays to these additional operators. We assume that the Wilson coefficients of the dimension-five operators are generated at Λ = 10TeV, while all other dimension-six Wilson coefficients are zero at this scale. To avoid constraints from the observed neutrino masses we consider the scenario where the second Higgs doublet has a negligible vacuum expectation value and a mass at the weak scale. The Higgs sector could be assumed to be close to that of an inert two-Higgs doublet model [26,27,28,29] and the dangerous couplings λ 6 , λ 7 and Y (2) are not generated radiatively. Renormalisation group running will then generate non-zero Wilson coefficients of several dimensionsix operators at µ ∼ v. Only those dimension-six operators that involve standard model particles are of interest to us, since the vanishing vacuum expectation value of the second Higgs doublet will suppress the contribution of the other operators after spontaneous symmetry breaking. Applying the constraints of Table 1 of the Wilson coefficients evaluated using eqn (4.4) neglecting the small log ln(m 22 /m W ), we find the following: the µ → 3e decays provide the greatest sensitivity to the additional dimension-five Wilson coefficients. In particular the left-handed contribution implies where we neglected the mixing of the dimension-five operators amongst themselves, as this would contribute at two-loop order to the lepton flavour violating processes. For the right-handed contribution we find which exhibits a weaker sensitivity. It is interesting to note that current experimental data is already sensitive to this parameter space of Wilson coefficients. The contribution to the µ → eγ is further suppressed by the smallness of the Yukawa couplings which puts these beyond current experimental sensitivity. We also checked that the current experimental situation for τ decays does not lead to significant constraints.

Summary
Motivated by neutrino masses and the expected progress in searches for lepton flavour violation, we calculated the leading one-loop contribution of a pair of lepton number violating dimension-five operators to the coefficients of lepton flavour violating dimension-six operators. The diagrams are given in figure 1. The dimension-five operators that we considered are the Weinberg operator, constructed out of SM fields and given in eqn (2.4), and three additional dimension-five operators that can be constructed in the Two Higgs Doublet Model, given in eqn (2.13). The dimension-six, lepton flavour violating operators of the SMEFT are listed in Appendix D, in the "Warsaw" basis, and the subset of these operators relevant for our calculation is given in eqn (2.7). A selection of constraints on their coefficients, evaluated at the weak scale, is given in Appendix E.
In section 3, we obtain the anomalous dimensions mixing two dimension-five operators into the lepton flavour violating operators of eqn (2.7). The required counterterms are given in eqns (3.2-3.7) for the Standard Model with a single Higgs, and in eqns (3.8-3.13) for the case of the Two Higgs doublet model. Then in section 3.3, we outline the derivation of the renormalisation group equations (Appendices B and C present our calculation and the flavour dependence of our result in more detail), and the resulting anomalous dimensions are listed in eqns (3.17-3.20). The mixing of two dimension-five operators into the lepton flavour conserving four-Higgs operator, via the diagram of figure 2 is given in Appendix G; however, we do not consider mixing into dimension six operators constructed with the second Higgs of the 2 Higgs Doublet Model. This completes the one loop renormalisation group equations of the standard model effective theory, up to operators of dimension six. It is amusing that the insignificant effect we calculate does not involve Standard Model couplings, so, in an expansion in terms of SM couplings, our result is the "leading" contribution to the one-loop RGEs of the dimension-six SMEFT 5 .
In the effective field theory constructed with Standard Model fields, the coefficient of the Weinberg operator is proportional to the neutrino mass matrix. So the lepton flavour changing amplitudes induced by double insertions of the Weinberg operator are ∝ (m ν /m W ) 2 ln Λ/m W , and far below current sensitivities. This is outlined in section 4.1. However, the situation is different in the 2 Higgs doublet model, as discussed in section 4.2: there are four operators at dimension five, and the neutrino mass matrix only constrains one combination. We evaluated the mixing of the four operators into lepton flavour violating operators of the standard model effective theory, and for a lepton number violating scale of 10 TeV we found that the current experimental value of µ → 3e is sensitive to the Wilson coefficients of these additional operators.
as [33] ψ so the amplitude M f i is where the SU(2) lepton indices are lower case, Higgs indices are upper case, ℓ αj and ℓ c αj represent a final state lepton and an initial state anti-lepton respectively. The factor i is the usual factor for Feynman rules and the factor (−i) is due to the calculation of M f i . This expression agrees with Feynman rule of Reference [21].
A Feynman-rule to attach a W -boson to the ℓ c line also will be needed. With the following identities [25] ℓ c = Cℓ one obtains (where the (-1) is for interchanging fermions) and recall that τ = τ † , so τ * = τ T .  Note that we have chosen a convention for our Feynman rules to eliminate any dependence on the momentum of the incoming lepton, p i , since all momenta are not independent.

B.1 Flavour dependence
We allow for multiple operators at both dimension-five and -six, and denote a particular Wilson coefficient by C ζ X , where X and ζ are the operator and flavour labels respectively. Then the bare Wilson coefficients of the dimension-six standard model effective theory Lagrangian can be written as where ζ, η and θ represent generation indices of an operator, and the renormalisation constants Z ζθ XY encode the mixing of dimension-six Wilson coefficients amongst themselves, which can be extracted from the anomalous dimensions of reference [19]. In the standard model, the mixing of two dimension-five Wilson coefficients and write the generation summation in the case of an operator involving four fermions explicitly as: The sum over generation indices reduces trivially for operators that involve less fermions. The corresponding renormalisation equation ensures that the pole of the one-loop off-shell matrix element of an insertion of two dimension-five operators is cancelled by its counterterm. Factoring out the common overall factor C αβ 5 C δγ *

5
we write: 1/ǫ denotes the 1/ǫ pole of a one-loop diagram and f | and |i are arbitrary off-shell final and initial states.
In calculations of the loop diagrams the following generation structures arose:

B.2 Four-lepton Green's function
In the following we will explicitly present the renormalisation of a Green's function involving four lepton doublets. When we consider double-insertions of dimension-five operators one additional operator that vanishes in the limit d → 4, a so-called evanescent operator, appears in our calculation. The exact definition of the evanescent operator in d dimensions is not important, but will induce a scheme dependence beyond one-loop. We use O αβγδ where the first term has a left-right chirality structure and i, j, k, l are SU(2) indices.
Denoting the flavour and SU(2) component of the final state f | = ℓ k,φ ℓ l,χ | and the initial state |i = |ℓ i,ψ ℓ j,ω by φ, χ, ψ, ω, and i, j, k, l respectively, we find for the third diagram of figure 1 f |Q αβ 5 (Q γδ 5 ) † |i | 1/ǫ = (ū ψi P L v ωj )(v φk P R u χl ) 64π 2 (δ ψδ δ ωγ + δ ωδ δ ψγ ) (δ χα δ φβ + δ φα δ χβ ) (δ il δ jk + δ ik δ jl ) , (B where we have used the hermiticity condition of the renormalisation constants 6  and O βα v (3) . The internal Higgs and lepton lines of the loop diagram may couple to B µ or W a µ bosons of the U(1) Y and SU(2) L groups respectively. Since the group structure of U(1) Y is trivial, we concentrate here on the calculation resulting from emission of a W a µ boson. The results for emission of B µ emission may be retrieved from these results by replacing the SU(2) L generators everywhere by U(1) Y generators, 1 2 τ a ij → y H,ℓ δ ij at the beginning of the calculation.
The renormalisation equation for the process where the tree-level matrix elements are replaced by their respective amplitudes and the SU(2) algebra has been simplified. Two diagrams must be evaluated for the double insertion of dimension-five operators with associated emission of a W a µ boson, which can couple to either the internal Higgs or internal lepton. These diagrams are denoted by D 1 and D 2 , and are shown in figure 7. Figure 7: Double insertions of dimension-five operators with associated emission of W a µ that mix into dimensionsix operators.
Calculating the diagrams and isolating the 1/ǫ poles gives , (B.10) we find the total amplitude of the double-insertion of dimension-five operators: In this form it is simple to set up simultaneous equations for the renormalisation condition by comparing the loop and tree amplitudes, (B.14) This underconstrained set of equations may be constrained by substituting in solutions for Z γδηκ,βα 55,v (1) and

C Renormalisation Group Equations
The bare Wilson coefficients of dimension-five operators can be written as is the renormalisation matrix, and µ is the renormalisation scale. The µ 2ǫ introduces an additional term proportional to ǫ into the d-dimensional renormalisation group equation This reduces to the renormalisation group equation in d = 4 dimensions where the 4-dimensional anomalous dimension matrix is independent of the choice of the overall factor µ 2ǫ . Therefore the µ 2ǫ term can be neglected when only considering mixing amongst operators of equal dimensions. In the case of mixing between operators of different dimensions a more careful treatment is required. At loop level, operators of different dimensions can mix via multiple operator insertions [31]. Consider the specific case of loop diagrams involving two dimension-five operators mixing into diagrams with a single dimension-six operator insertion. We denote dimension-six quantities with a tilde, quantities that mix dimension-five and -six with a hat, and dimension-five quantities without a tilde or hat. The bare dimensionsix Wilson coefficient isC whereC bare is µ-independent. Therefore the renormalisation group equation is whereγ θη Y X is defined analogously to equation (C.2), and where the explicit form in terms of generation indices is [γ αβ γδ AB ] † = [γ βα δγ AB ] * and δ αβ γδ AB = δ AB δ αγ δ βδ . The terms in the second line of the above equation only contribute beyond one-loop. Furthermore, the contribution to the renormalisation tensor Z ζθ,υ AB,Y is µ independent at one-loop and only the term proportional to 2ǫ contributes in our calculation. A comment regarding the sign of the 2ǫ contribution is in order. The factor in µ 2 ǫ in (C.5) generates a term proportional to −2ǫ, while the derivative of the dimension-five Wilson coefficients generates a contribution proportional to 2 × 2ǫ from (C.2). Hence the one-loop anomalous dimension matrix reads

D Operators
This Appendix lists dimension-six, SM-gauge invariant operators that change lepton flavour.The operators are in the Buchmuller-Wyler basis, as pruned in Grzadkowski et.al. [20], commonly refered to as the "Warsaw" basis. All operators are added to the Lagrangian +h.c., as given in eqn (2.6): where the flavour indices are represented by ζ, and are all summed over all generations. In the conventions of [20] and [19], the hermitian conjugate is not added for "self-conjugate" operators, for which ζ C ζ . So we define such operators with a factor 1/2 to avoid this double-counting.
The four-fermion operators involving β ↔ α flavour change and two quarks are: where ℓ, q are doublets and e, u are singlets, n, m are possibly equal quark family indices, and A, B are SU (2) indices. The operator names are as in [20] with ϕ → H; the flavour indices are in superscript.
In the case of four-lepton operators, the flavour change can be by one or two units. Notice that in the case of O ee and O ℓℓ , which are symmetric under interchange of the two bilinears (eγ µ µ)(τ γ µ τ ) = (τ γ µ τ )(eγ µ µ), there will be two equal coefficients that contribute to the Feynman rule: Then there are the operators allowing interactions with gauge bosons and Higgses. This includes the dipoles, which are normalised with the muon Yukawa coupling so as to match onto the normalisation of Kuno-Okada [2]: where Y β denotes the Yukawa coupling of a charged lepton e β in the mass basis, the double derivatives are defined in eqn (2.8), and we include factors of 1/2 for hermitian operators as discussed above eqn (D.1).

E Experimental bounds on coefficients
The aim of this appendix is to obtain experimental constraints on the coefficients of the LFV operators of eqn (2.7), evaluated at the weak scale m W . We are interested in this subset of operators because they are generated at one loop by double-insertions of dimension-five, lepton number changing (LNV) operators. Such constraints will allow an estimation of the sensitivity of LFV processes to the coefficients of LNV operators.
Recall that constraints and sensitivities are different. A constraint is an exclusion, which tells the range of values a coefficient cannot have. For instance, the dipole coefficient (evaluated at the muon mass scale) C eµ D,R (m µ ), cannot be larger than 1.05 × 10 −8 because the branching ratio searched for by the MEG experiment [3] is , and the current experimental search imposes this constraint. Sensitivity is often discussed when an observable depends on many coefficients, and gives the range of values where a coefficient could have been seen. For instance, among the many loop processes that contribute to µ → eγ, there are two-loop diagrams involving flavour-changing Higgs coupling C eµ eH (m W ). Calculating these diagrams and imposing that they saturate the current experimental bound gives eα e y t 8π 3 y µ C eµ eH (m W ) = 1.05 × 10 −8 .
Smaller values of C eµ eH are allowed(the experiment could not have seen them), but larger values are not excluded by MEG, because many other operator coefficients could contribute to the rate, with possibly cancellations.
The difference between an exclusion and a sensitivity is illustrated in figure 8, where the allowed region is the diagonal ellipse. The horizontal variable x is excluded outside the projection of the ellipse onto the x-axis (where the axis is thickened). But the experiment is only insensitive to x inside the intersection of the axis with the ellipse (dashed red line). Values of x between these two regions are allowed, provided that y has the appropriately correlated value.
Three ways to relate low-energy experimental bounds to the coefficients of operators at a higher scale are: 1. to calculate the sensitivity of an experimental process to a particular operator coefficient. This is usually simple.

2.
To express an experimental rate as a function of high-scale coefficients. This is slightly more difficult, because more coefficients are involved: each coefficient that contributes at the experimental scale will become a linear combination of high scale coefficients due the renormalisation group mixing.
3. To obtain constraints on coefficients at the high scale. This is more involved, because a sufficient number of experimental constraints must be combined, in order to obtain a finite allowed region in coefficient space (no "flat directions"). Then the allowed region must be projected onto the various axes, in order to obtain constraints.
The third option is the most useful, but beyond the scope of this work. Instead here, we partially follow the second option, as a contribution to the third: we consider experimental bounds on the dimension-six operators which are generated in RGE evolution by double-insertions of dimension-five operators that change lepton number. We aim to quote these bounds at m W . The processes in question are LFV Higgs and Z decays (which occur at the weak scale), and flavour-changing lepton decays at low energy (these bounds must be translated to the weak scale via the RGEs of QED and QCD). So we will not succeed in our aim of setting constraints on coefficients at m W , because the low-energy experimental bounds depend on many coefficients at the weak scale, and we do not include enough experimental bounds. In the following sections, we outline the calculations of the various rates, and summarise the experimental constraints on coefficients at m W in table 1.   (3) , O ℓℓ , O ℓe ) also apply to the conjugate coefficient. All the bounds apply to running coefficients evaluated at m W , and are for Λ = v ≃ m t . The combination of coefficients C penguin is defined in eqn(E.12) and before eqn (E.23), δ is defined after eqn (E.23), and g e R = 2s 2 W , g e L = −1 + 2s 2 W .

E.1 Rates and calculations
When the Higgs gets a vev, the "penguin" operators O Hℓ (1) and O Hℓ(3) generate a vertex involving the Z and two charged leptons. If the flavour-changing Z-fermion vertex is written in a SM-like form : −l α Z µ g 2cW γ µ (g V − g A γ 5 )l β , then The branching ratio can be written where 2.5 GeV is the Z width in the SM. Since O Hℓ (1) and O Hℓ(3) are hermitian, the conjugate process Z → l β l α neccessarily occurs at the same rate, so the BR to the experimental final state is and the bounds we obtain on the operator coefficients, evaluated at ∼ m W , are given in table 1.
The flavour-changing Higgs decays occur via the non-hermitian operator O eH . When the Higgs has a vev, it induces the Feynman rules for a flavour-changing Higgs vertex with two fermions: We calculate the flavour-changing branching ratio by comparing to BR(h → bb) = 0.575 ± 0.32 (from the Appendix of the Higgs Working Group Report [35], for m h = 125.1 GeV), assuming the Feynman rule for hbb We use a one-loop approximation [15] for the running b mass GeV. The operator O eH is not hermitian, but is always included in the Lagrangian +h.c.. So C eµ eH O eµ eH + h.c. will induce both h → e L µ R and h → µ R e L at the same rate: where downstairs there is a 3 for quark colour sums, and a 2 from the chiral projectors in the lepton decay. The experimental search sums the e L µ R and µ R e L final states, so we obtain and the resulting contraints are given in table 1.

E.1.3 Including the low energy decays
The flavour-changing τ and µ decays listed in table 1 occur at energies ∼ m µ , m τ , so the decay rates are usually written in terms of the coefficients of dimension-six operators from the QCD×QED invariant basis appropriate at low energies. These "low energy" coefficients, which we denote with a tilde C, can be expressed in terms of SMEFT coefficients at m W by running them up to m W , then matching the QCD×QED-invariant operator basis onto the SMEFT. This was performed in [34] for µ → eγ, so we use the results of [34] for the radiative decays of section E.1.5. Reference [36] studied the Renormalisation Group evolution, below the weak scale, of the coefficients who mediate µ → eēe (as well those for as µ → eγ and µ → e conversion); we use these results, combined with the weak-scale matching conditions of [34], for the discussion in section E.1.4 of three body leptonic decays of τ s and µs. The minor differences between µ and τ decays are discussed in section E.1.4.
In the EFT below m W , we use the basis of lepton-flavour-changing four-fermion operators introduced in [2,34] for µ ↔ e flavour change 7 . The operators and coefficients have as subscript their Lorentz structure (V, S, T ) and the chiral projection operators of the two fermion bilinears, and the flavour indices of the four fermions as superscript. They wear tildes to distinguish them from the coefficients of SMEFT operators. We restrict to the dipole and vector operators, and neglect the scalars and tensors, which will turn out to be irrelevant for our study of LFV operators generated by double-insertions of LNV operators. So the four-fermion operator basis below m W is where αβ ∈ {eµ, µτ, eτ }, f ∈ {e, µ, τ, u, d, s, c, b}, and αβσρ ∈ {eτ eµ, µτ µe}. In addition, below m W we consider the photon dipole operators because the SMEFT operators O Hl(1) , O Hl (3) and O eH match onto the dipole at m W . The current bounds on µ → eγ, τ → eγ and τ → µγ will give the best sensitivity to the coefficients C Hℓ(1) , C Hℓ(3) and C eH .
The first step is to translate the experimental bounds into constraints on operator coefficients at the experimental scale. For the three-body leptonic decays of the τ , it is convenient to define ). Then BR(τ → 3l) can be directly compared to the branching ratio for µ → 3e [2]: where 2 √ 2G F = 1/v 2 and the generalisation to τ decays is straightforward, after accounting for 2s as we now discuss.
We calculate the decay rates in the approximation that all final state fermions are massless. Factors of 2 can arise when there are two identically-flavoured fermions in the final state: there will be 2 diagrams, and a factor of 1/2 in the final-state phase space. Then there are two cases: a) if the identical fermions have the same chirality, there is constructive interference between the two diagrams (despite the fact that they have relative minus signs due to Fermi statistics), which doubles the rate. (This is consistent with µ → 3e rate of Kuno and Okada [2] given above.) b) if the fermions have different chirality, the interference is suppressed by final state masses (which are neglected), so the two for two diagrams cancels the 1/2 from phase space.
We set the dipole coefficients to zero, because they are better constrained by the radiative decays discussed in the next subsection (see table 1). Then it is clear that each upper bounds on a three-body leptonic decay of the τ or µ, implies six independent constraints on operator coefficients (evaluated at the experimental scale), those of interest to us are given in table 2.
C eτ ee V,LL < 2.8 × 10 −4 , C eτ ee V,LR < 4 × 10 −4 τ → eµµ 1.6 × 10 −7 C eτ µµ V,LR , C eτ µµ V,LL < 4 × 10 −4 τ → µee 1.0 × 10 −7 C µτ ee V,LR , C µτ ee V,LL < 3.2 × 10 −4 τ → µµµ 1.2 × 10 −7 C eτ µµ V,LL < 2.5 × 10 −4 , C eτ µµ V,LR < 3.5 × 10 −4 τ → eeµ 8.6 × 10 −8 The operator coefficients C X (m τ ) given in table 2 can be expressed in terms of coefficients at m W using the one-loop RGEs [34,36]: where [γ e ] is the one-loop anomalous dimension matrix of QED, ln mW mτ = 3.85, ln mW mµ = 6.64 and the approximate solution neglects the running of α e . The one-loop QED corrections involve photon exchange between two legs of the operator, which does not change the flavour or chiral indices, and also "penguin" diagrams, where two legs of the operator are closed in a loop, and a photon is attached, which turns into two external leg fermions. The "penguins" can change the chirality and flavour, and allow 2-lepton-2-quark operators to mix with the four-lepton operators. We therefore need a recipe for dealing with the quark-sector threshholds m b , m c and Λ QCD . We make the simplest approximation, which is to have a single low-energy threshhold at m τ , and run from m W → m τ with five flavours of quark, and we use this low-energy scale also for the decays of the µ. In this approximation, it is convenient to define the combination of operator coefficients where l ∈ {e, µ, τ }, q ∈ {u, d, s, c, b}, and Q q is the electric charge of the quark. Then the coefficients constrained in table 2 can be written Finally, the combinations of coefficients that are constrained by data can be matched at m W onto coeffi-cients of SMEFT operators [34] 8 : where l ∈ {e, µ} in the above equations, and g e R = 2 sin 2 θ W , g e L = −1 + 2 sin 2 θ W . In order to match the "penguin" coefficient of eqn (E.12) onto coefficients of the SMEFT, matching conditions for operators with a quark bilinear are also required: where all the coefficients are evaluated at m W , and δ = αe 4π ln mW mτ ∼ 1/400. This and other constraints from 3-body τ decays are given in table 1, where for compactness, [1 ± 12δ] is approximated as 1.
The radiative decays l β → l α γ provide some of the most restrictive bounds on lepton flavour violation. The branching ratio at m β can be written where the low energy dipole operators are added to the Lagrangian as in eqn (E.9). 8 These equations differ from [34] due to different conventions for operator normalisation and signs, and also due to some errors in [34]. The SMEFT basis used here is normalised according to [20], where there are "redundant" flavour-changing four-fermion operators, which are absent from the basis used below m W in [34], compare e.g. the left and right hand sides of eqn. (E.21). Then, the sign convention used here for the g f L,R and the Z-vertex Feynman rule agrees with the PDG but is opposite to that of [34]. Finally, in [34], there is an incorrect factor of 2 mutiplying the penguin coefficients which generate s and t channel diagrams; this 2 should not appear, because the four-fermion operator generates the same s and t channel diagrams.
The dipole coefficients evaluated at the experimental scale can be expressed in terms of SMEFT coefficients at the weak scale as [34] C αβ D,L (m τ ) = C βα * eγ (m W ) + eα e y t 8π 3 y µ C βα * eH (m W ) + where the contributions of scalar and tensor four-fermion operators were neglected, g e R and g e L are defined after eqn(E.21), and C αβ eγ = c W C αβ eB − s W C αβ eW . (E.27)

F Comparison with the Literature
The standard model calculation has been performed in Reference [21] in a different operator basis. We disagree with their final results even after transforming our results to their basis. To do this we specify our basis T , (F.1) and the one used in Reference [21] where the additional operators are defined as only the first two rows ofR have entries that are not proportional to an identity transformation. These two rows are determined by the following linear transformation 9 : The Wilson coefficients and renormalisation constants will consequently fulfil our hermiticity conditions in our basis, but not necessarily in the basis of Reference [21]. The counterterms of the Wilson coefficients transform in the same way as the respective Wilson coefficients under our change of basis, i.e. as δc =R T δC , where δC = (16π 2 )ǫ CZ C † represent the counterterms multiplied with (16π 2 )ǫ, while δc correspond to the analogous expression in theQ basis. Using the counterterms presented in Equations (3.2), (3.3) and (3.6), and the results of (F.8) we obtain δc = − 5 2 C 5 C * 5 , − 1 2 C 5 C * 5 , which fulfil the hermiticity condition of the overall Lagrangian, even though this is not immediately apparent due to the choice of basis. These results are in disagreement with the final results quoted in Reference [21]. Yet using the results quoted in the individual diagrams in Appendix B of Reference [21] we find agreement with the expression of eqn (F.7) apart from a global minus sign, which suggests that a different projection was performed. We explicitly checked that the diagrams given in Appendices B.1, B.3 and B.4 of their calculation have the opposite sign compared with our calculation, while we agree with their lepton conserving contribution presented in Appendix B.2. Following the explanations of the calculation it appears that part (the δδ part) of the diagram evaluated in Appendix B.1 of Reference [21] is projected onto an operator basis where the operators Q we find Again, this result does not agree with Reference [21]. Finally, note that projecting the results quoted for the individual diagrams in Appendix B of Reference [21], except the εε part, would give while projecting only the εε part on the non-hermitian basis yields δc εε = +2 C 5 C * 5 , 0, 0, 0 T . Summing these two terms would reproduce the results of Reference [21].

G Flavour Conserving Contribution
Even though the diagram in figure 2 cannot induce lepton flavour violation it contributes to the renormalisation of O eH and the corresponding operators that involve quarks. We also explicitly checked that the diagrams that involve six external Higgses vanish after summing over them. Denoting the trace over the product of the two dimension-five Wilson coefficients by Tr[C 5 C * 5 ] we find where Y u and Y d are defined as Γ u and Γ d of reference [20]. In addition we also generate the following mixing into operators that only comprise Higgs and gauge fields and write where the additional operators are defined as: