Gearing up for the next generation of LFV experiments, via on-shell methods

Lepton Flavor Violating (LFV) observables such as $\mu\to e\gamma$, $\mu\to 3e$ and $\mu N \to eN$ are among the best probes for new physics at the TeV scale. In the near future the bounds on these observables will improve by many orders of magnitude. In this work we use the SM EFT to understand the impact of these measurements. The precision reach is such that the interpretation of the bounds requires an analysis of the dimension-six operator mixing up to the two-loop level. Using on-shell amplitude techniques, which make transparent many selection rules, we classify and calculate the different operator mixing chains. At the leading order, on-shell techniques allow to calculate anomalous dimensions of SM EFT operators from the product of tree-level amplitudes, even for two-loop renormalization group mixings. We illustrate the importance of our EFT approach in models with extra vector-like fermions.


Introduction
In the SM, lepton numbers L e,µ,τ arise from accidental symmetries that guarantee the absence of neutrino masses and lepton flavor violating (LFV) processes. These symmetries are, however, not respected by higher-dimensional operators beyond the SM suppressed by some new physics scale Λ. In fact, already at the first order in the 1/Λ expansion (dimension-five operators), we find that lepton number can be broken in two units, ∆L i = 2, generating then neutrino masses. The smallness of the neutrino masses however forces the scale suppressing these operators Λ to be extremely large (or the corresponding Wilson coefficients extremely small), making the physical effects of these operators only visible in neutrino physics.
We are here interested in LFV processes where the relative lepton numbers are violated, but the total lepton number is preserved. These processes are among the best indirect probes for new physics at the TeV [1]. This is due to the fact that L e,µ,τ are not easy to arise accidentally in scenarios beyond the SM (BSM), where the matter content and interactions are much more extended than the simple leptonic sector of the SM. As a consequence, sizable BSM effects are expected to arise in LFV processes. They can be characterized by dimension-six operators (∼ 1/Λ 2 ) that, as they preserve the total lepton number, do not need to be as highly suppressed as dimensionfive operators. For this reason they can have an important impact in future LFV experiments involving the charged lepton sector.
We will consider in particular LFV processes with ∆L e = ∆L µ = 1. At present, the most competitive experimental measurements come from the processes µ → eγ, µ → eee and µN → eN (see Table 1), and these will improve in the near future. Especially, the sensitivities of µ → eee and µN → eN are expected to improve by four orders of magnitude by the mid-2020s [2]. These LFV processes are very clean observables from the theoretical point of view, since the contributions of the SM (even when implemented with neutrino masses) are much smaller than the present and BR(µ → eγ) BR(µ → eee) R(µN → eN ) BR(h → µe) Current 4.2 · 10 −13 [33] 1 · 10 −12 [34] 7 · 10 −13 [35] 6.1 · 10 −5 [36] Future 6.0 · 10 −14 [37] 1 · 10 −16 [38] 8 · 10 −17 [39]  future experimental sensitivities. We will also consider h → µe, although, as we will show, it does not place any important constraint. There are other LFV processes such as Z → µe and J/ψ → µe, but they are not competitive with those in Table 1. We also have processes involving flavor violations in the quark sector, e.g. K → µe, but these quark flavor transitions are quite constrained from other non LFV processes, so we will not consider them here.
We will use the Effective Field Theory (EFT) approach to characterize the BSM contributions to LFV processes as they are captured in the Wilson coefficients of the dimension-six operators. Our purpose is to understand at what order in the loop expansion the Wilson coefficients can enter into the different observables of Table 1. Depending on the experimental sensitivity, we will see that certain Wilson coefficients can be better bounded by µ → eγ even when they enter at the two-loop level.
This program will therefore require the calculation of anomalous dimensions at the two-loop level. We will use on-shell methods to perform these calculations. It has already been shown the efficiency of these methods to calculate anomalous dimensions at the loop level [3][4][5][6][7][8][9], where these can be reduced to products of tree-level amplitudes integrated over a phase space. Onshell amplitude methods are also able to show in a transparent way many selection rules hidden in the ordinary Feynman approach [10][11][12][13][14][15][16][17], mainly due to helicity [11] or angular momentum conservation [14,15]. This simplifies substantially the loop calculations. 1 Our calculation of the anomalous dimensions at the two-loop level will in particular show how µ → eγ can constrain the LFV couplings W µν e and hµe and some four-fermion interactions at an unprecedented level. Previous (partial) analysis can be found in Refs. [22][23][24][25][26][27][28]. 2 In Section 2 we review the LFV dimension-six operators of the SM EFT and the tree-level bounds derived from the processes in Table 1. In Section 3 we analyze which operators mix into the dipole operators responsible for µ → eγ . We find that there are various important two-loop anomalous dimension matrix elements that are unknown. Section 4 is then dedicated to compute these missing pieces of the two-loop renormalization group (RG) equation, via on-shell methods. In Sections 5 and 6 we analyze the RG mixing of dimension-six operators into the µ → eee and µN → eN processes, respectively. A global perspective of the new loop effects we have found from mixings into the LFV observables is summarized in Section 7. For illustration, in Section 8 we show the impact of our analysis on two simple UV completions of the SM EFT. Finally, we conclude in Section 9.

LFV dimension-six operator basis
The relevant dimension-six operators for our analysis of processes with ∆L e = ∆L µ = 1 are given by L τ a σ µν e (1) R HW a µν + C µe DB y µ g Λ 2L (2) R HB µν + (µ ↔ e) (2.1) is the SM muon Yukawa coupling and g, g are the SU(2) L and U(1) Y couplings as defined in Appendix B. We also have F L = Q L , L L and f R = u R , d R , e R the SM left-handed SU(2) L doublets and right-handed SU(2) L singlets respectively. We only specify the flavor indices for the µe transitions where the superindices 1, 2 corresponds to e, µ. We also have factored out a y µ in operators involvingL R where the µe transitions have a change of chirality. Therefore, in the interchange µ ↔ e in L 6 we keep fixed y µ .
The four-fermion operators have been separated into two groups: Eq. (2.4) are operators of type ψ 2ψ2 in Weyl notation, while Eq. (2.5) are operators of type ψ 4 (andψ 4 ). They have respectively total helicity h = 0 and h = ±2. This is important as we will see later for understanding the renormalization mixing. As compared to the Warsaw basis [40], we have made the replacement 3 and we have written differently the operator (via Fierzing) to make it clear that this is an operator of the type ψ 2ψ2 with total helicity h = 0.
It will be important for later to understand the impact of the operators (2.2) to write them in 3 Using Fierz identities, we have the relation O the unitary gauge: which shows that these operators induce a LFV coupling for the Z, W and hW, hZ. Custodial symmetry together with a L ↔ R parity can protect some of these couplings. For example, as explained in Ref. [41], these symmetries can impose C L + C L3 = 0 such that the LVF Zµe and hZµe couplings are zero. Therefore, for BSM sectors respecting these symmetries, the only LFV couplings will involve neutrinos which are difficult to detect, implying that no strong bounds can be derived from direct measurements such as W → µν e , eν µ or Z → ν µνe . As we will see, loop effects will be important to get better bounds from other LFV observables.

LFV experimental constraints at the tree level
The Wilson coefficients of the operators (2.2)-(2.5) are constrained from the LFV observables of Table 1. In this section we want to review the bounds obtained from a tree-level analysis. The result is presented in Table 2 in black color (see Refs. [22][23][24][25] for previous analysis). The aim of this analysis is to understand which loop effects mixing the different Wilson coefficients of (2.2)-(2.5) are important to consider in order to obtain better bounds. These loop calculations and the discussion of how much the bounds are improved will be presented in the next sections.
The current constraint on the dipoles C eµ,µe DW,DB of Eq. (2.1) is dominated by the experimental bound on BR(µ → eγ) at which they enter at tree level (see Eq. (3.3)). This leads to 1 and will improve by almost an order of magnitude in the future. Eq. (2.9) provides a very strong constraint that shows that renormalization effects to (C DW − C DB ) arising from other Wilson coefficients C i could also place strong constraints on these latter. These effects can be important even when they arise at the two-loop level, e.g.
where we estimate to obtain C i /Λ 2 1/(6 TeV) 2 . Therefore an analysis of anomalous dimension mixings into C DW,DB is needed up to the two-loop level. The Wilson coefficient (C DW − C DB ) also enters at tree level in the observables µ → 3e and µN → eN (see Eq. (5.1) and Eq. (6.2) respectively), but the present constraints are not so competitive, as shown in Table 2. Nevertheless, the spectacular prospects of improvement by four orders of magnitude in each of them will lead to an improved bound on (C DW − C DB ) that could overcome that from µ → eγ.
Let us consider now the four-fermion interactions (2.4) and (2.5). In particular, the operators of type µēēe enter at tree level to µ → 3e and provide the limits  They are clearly not as strong as that in (2.9), but in the near future this bound is expected to improve by two orders of magnitude, being then comparable to Eq. (2.9). Therefore loop effects mixing other Wilson coefficients into C µeee LL,RR,LR,RL should also be considered. As we will see, however, all the relevant LFV Wilson coefficients enter into the anomalous dimensions of C µeee LL,RR,LR,RL at the one-loop level, making two-loop effects not necessary. Similar conclusions can be reached for four-fermion interactions of type µēūu and µēdd that are constrained by µN → eN at a similar level as in Eq. (2.10).
The Wilson coefficients of Eq. (2.2) give rise to the LFV Zµe coupling (see Eq. (2.8)), and therefore can enter at tree level into the observables µ → 3e and µN → eN (see Eq. (5.1) and Eq. (6.2) respectively), providing a bound similar to (2.10). There are also direct bounds from Z → µe but these are not so competitive (one gets constraints of order Λ 5 TeV). It is important to remark again that the induced Zµe coupling is proportional to the combination (C µe L +C µe L3 ), and therefore tree-level constraints from µ → 3e and µN → eN are only on this specific combination. As we already discussed, the orthogonal combination (C µe L −C µe L3 ) only enters at tree level into LFV gauge boson couplings involving neutrinos, implying that no relevant bounds on this combination can be placed at this order. Loop effects however can be important as we will discuss in the next section. The parametrization of the bounds using these two particular combinations has also a theoretical motivation, since, as we already mentioned, BSM models can give sizable contributions to (C µe L − C µe L3 ) but not to (C µe L + C µe L3 ) as this latter is protected by a custodial symmetry [41]. We provide an example in Section 8.1.
The last operator to consider is (2.3) that can only enter at tree level in Higgs physics, specifically in h → µe. 4 The bound (see Table 2) is very weak, therefore it is expected that bounds from other observables, where C µe,eµ y can enter via mixing at the loop level, can be potentially stronger. We will indeed see that this operator enters at the two-loop level into µ → eγ and gets a much more severe constraint.
We conclude that there are some Wilson coefficients whose loop effects can be potentially relevant for the observables µ → eγ, µ → 3e and µN → eN . Below we will calculate the anomalous dimensions at the two-loop level for these C i that will allow us to obtain new competitive bounds.

µ → eγ in the SM EFT
The µ → eγ process arises from the Lagrangian term where d µe,eµ must be evaluated at the muon mass scale. The branching ratio is given by [42] BR (µ → eγ) = 384π 2 |d µe | 2 + |d eµ | 2 , where the large factor is because the LFV decay is a two-body process, while the dominating muon decay channel is a three-body one. This decay process is tightly constrained and further sensitivity is expected in the near future, see Table 1.
At tree level the only Wilson coefficients of the SM EFT that enter into d µe are the dipoles C µe DW,DB : where e = g sin θ W , and similarly for d eµ by exchanging µ ↔ e.
We are interested in operator mixing into the dipoles (3.1). There are various effects to be considered: (i) finite matching contributions from a new physics scale Λ, (ii) RG mixing from Λ to the electroweak scale m W , (iii) finite threshold corrections arising from integrating out the W , Z, Higgs and the heavy SM fermions, and (iv) RG mixing from the electroweak scale m W to m µ . At the leading order (one-loop level) the UV matching threshold (i) and SM IR matching (iii) lead to finite contributions to the dipoles. These type of corrections are model dependent, and can (partially) cancel against each other as they involve rational coefficients, apart from some overall couplings. We will show examples in which this occurs in Section 8.1. Here we are instead interested in contributions from RG running, which correct the tree-level dipoles by logarithms ln(Λ/m W ), raised to a suitable power, which are hardly possible to cancel against matching contributions of either type (i) and (iii). In particular we are interested in the leading RG mixings from any operator, other than the dipoles, into the dipoles. Our aim is to understand to which physics a precise measurement of (3.2) is sensitive. Therefore, we will only keep track of the leading RG mixing contribution from any dimension-six operator into (3.1). At the one-loop level, selection rules, mainly based on spurious SUSY [10] or helicity arguments [11], dictate the anomalous dimension mixing terms O i → O j . Of particular importance is the selection rule ∆n |∆h| [11], where ∆n = n f − n i and ∆h = h f − h i , being n (h) the number of fields (helicity) of the operator. There is only one exception to this rule,ψ 2 ψ 2 ↔ ψ 4 only when the loop involves the Yukawa product y u y d or y u y e . In the particular case of dipole operators that have n = 4, |h| = 2, the only operators that can mix into them are those with helicity |h| = 2. This reduces to the four-fermion operators in (2.5), apart from the orthogonal combination (C µe DW +C µe DB ). The one-loop RG mixing from O µeqq LeQu however is trivially absent becauseL (2) e (1) is external to the loop calculation and does not have the dipole structure (and similarly for O eµqq LeQu ). This argument leaves out only the following possibility for one-loop leading-log contributions: where we suppressed flavor indices. The calculation of the one-loop RG mixing in (3.4) was done in [43] and recently revisited in [5] with on-shell methods. The result is presented in Section 4.1.
At the two-loop level, there are other various dimension-six operators that can mix with the dipoles. We can classify them according to a logarithmic expansion. At leading order we have two-loop double-log contributions proportional to  Figure 1: Structure of the RG mixings from dimension-six operators into LFV dipoles. One-loop RG mixing can only proceed from left to right, generating mixing between operators belonging to different helicity boxes (shown in grey), and within operators belonging to the same helicity box. One-loop mixing from the category (n, h) = (4, 0) into (4, 2) is forbidden by spurious SUSY [10] or helicity [11] selection rules, except for the non-holomorphic contribution proportional to y u y d or y u y e , which mixes the ψ 4 andψ 2 ψ 2 four-fermion operators at one loop. In blue and red we indicate the relevant RG mixings for the LFV dipoles at one-and two-loop precision, respectively. product of Yukawa couplings y u y d or y u y e . There are double logs of the type C 2 i ln 2 (Λ/m W )/(16π 2 ) 2 as well. However these are sub-leading corrections to already existing one-loop effects. The second type of corrections are two-loop single-log contributions There are only three of such two-loop mixings: ψ 2 φ 3 , ψ 2ψ2 → ψφF , which were computed in [44,4], andψγψH † DH −→ F µνψ Hψ, which will be calculated here for the first time.
We summarize the structure of the one-and two-loop mixing pattern into dipoles in Figure 1. In the next section we present the calculation of the relevant anomalous dimensions of dipole Wilson coefficients up to the two-loop level. Our task is greatly simplified by using on-shell methods as we explain next.

The on-shell way
Let us start presenting the basic concepts necessary for the computation of anomalous dimensions from on-shell amplitudes. Our calculations are based on the formula [3] which allows to compute the leading correction to the anomalous dimension matrix γ (1) ji with i = j. On the left-hand side (l.h.s.) of (4.1) we have a form factor computed in the free theory, i.e. at zero order in the coupling, denoted as minimal form factor. For instance, the minimal form factors of the dipoles are given by where and τ are tensors of SU(2) L . 5 The form factors n|O j |0 (0) are polynomials of the kinematical variables {|i , |i]}, i.e. the spinor-helicity variables. 6 On the right-hand side (r.h.s.), the symbol '⊗' denotes a dLIPS phase-space integration over the intermediate states m | m m| connecting the on-shell scattering amplitude M( n ← m) and the on-shell form factor m|O i |0 . The integral over the phase space is often referred to as unitarity cut, or cut for short. When evaluating the r.h.s. with Feynman diagrams, a cut corresponds to putting on shell the internal propagators crossing the cut, which has the net effect of cutting the Feynman diagram into lowerpoint diagrams sewed together through the phase-space integral. The upper script (0) means that the r.h.s. of (4.1) is computed to leading non-trivial order so that the unitarity cut n|M ⊗ O i |0 (0) leads to a polynomial in the kinematical variables, matching the left-hand side of (4.1). Further details and examples can be found in [3,4,7].

One-loop mixing into LFV dipoles
When only two particles cross the cut in the r.h.s. of (4.1), and when only four particles are involved in either the form factor or the amplitude, we are left with the following phase-space integral [3] O where the amplitude M describes the x + y → 1 + 2 scattering process, and as usual we have extracted a total delta function 12|M|xy = (2π) 4 δ (4) (p 1 + p 2 − p x − p y )M (12; xy). The particles of the form factor F O i (12 · · · n) ≡ 12 · · · n|O i |0 are all outgoing. The integral can be easily performed if we write the spinors in the integrand in a basis spanned by two of the external with the measure given by dΩ 2 ≡ dφ 2π 2 cos θ sin θdθ. The rotation of |x], |y] into |1], |2] is obtained by complex-conjugation of (4.5). When identical particles cross the cut in (4.4), one should include an extra combinatorial symmetry factor of 1/2! in the phase space.
As we already explained in the previous section, at the one-loop level we only have Eq. (3.4). This computation was done in [5] using unitarity cuts, agreeing with [43,44]. As a warm-up, we review next such computation. For the case at hand, Eq. (4.4) is given by where on the l.h.s. the gauge boson must be attached in all possible ways to the Higgs (dashed line) or the fermion lines (solid). This results in a two-to-two amplitude given by where T a k = g δ a k is the SU(2) L tensor arising from the contraction of left-handed doublets. On the right-hand side we have the form factor given by It is now a straightforward matter to plug (4.7) and (4.8) into (4.4), perform the spinor rotations (4.5) and a few elementary integrals, leading to where N c = 3 is the number of colors and y t the top Yukawa coupling. We recognize the minimal form factor of the dipole (4.2) and therefore the anomalous dimension is . In the case of mixing into O DW , we set hypercharges Y t R = 0 and Y H = 1 and change the SU(2) L tensor to be T a k = g(τ α /2) a k on the amplitude side. Then we get the following result: From the last expression we recognize the dipole (4.3) and the corresponding anomalous dimension.

Two-loop mixing into dipoles
We want to calculate here the two-loop mixingψγψH † DH −→ F µνψ Hψ, which is the only one relevant for µ → eγ not yet calculated. The two-loop leading-log contributions to the r.h.s. can in principle involve three-particle cuts or two-particle cuts: (4.11) The first diagram involves a tree-level amplitude and a tree-level form factor, so that the threeparticle cut accounts for the two-loop factor. The second/third diagram involves a tree-level/oneloop amplitude and a one-loop/tree-level form factor which, together with the two-particle cut, make it to two-loop order. Bellow we will show that the second and third diagrams do not contribute to theψγψH † DH −→ F µνψ Hψ mixings because of simple helicity selection rules. Thus, all our non-trivial calculations will only involve three-particle cuts. For the transition ψγψH † DH → F µνψ Hψ, in (4.11) we only need to consider two external particles to the scattering amplitude and form factor.
The phase-space integral involving the three-particle cuts can be nicely simplified into the following angular integration [3] where the amplitude describes the x + y + z → 1 + 2 scattering process at tree level. The spinors in the integrand can be rotated in terms of a basis spanned by the two external spinors: and the measure is dΩ 3 = 4 cos θ 1 sin 3 θ 1 dθ 1 2 cos θ 2 sin θ 2 dθ 2 2 cos θ 3 sin θ 3 dθ 3 dδ 2π dφ 2π . In the case that n identical particles cross the cut, one should account for an extra symmetry factor of 1/n! in the phase-space integral, as we will see in some examples in the next sections. Note also that (4.12) includes the −1/π factor in the r.h.s. of (4.1).

Top Yukawa y 2 t contributions
We expect this type of contributions to be the dominant ones because they are proportional to N c y 2 t . We first explain in detail the mixing of O eµ L into O eµ DB through a top loop. The three-particle cut is given by (4.14) The disconnected gauge boson notation means that the gauge boson must be attached anywhere in the left-hand side of the cut, i.e.

= + + +
Summing over all such possible attachments of the gauge boson leads to the following tree-level scattering amplitude A ba , (4.16) where A ba = g ba is a SU(2) L tensor arising from the contraction of the left-handed doublets. We have computed this amplitude using BCFW recursion relation [45,46] with a Risager shift [47]. 7 In (4.16) we have eliminated the hypercharge of the left handed fermions in favour of the Higgs hypercharge and right handed fermions hypercharge. The hypercharges are given by (Y µ R , Y t R , Y H ) = (−1, 2/3, 1/2). 8 The tree-level form factor on the r.h.s. of the cut is given by where the SU (2) tensor is B ba kl = −δ b k δ a l . Next we plug (4.16) and (4.17) into (4.12), perform the rotations (4.13) 9 , and after some simple algebra and elementary integrations we are led to We recognize the minimal form factor of the dipole (4.2).
Regarding the two-particle cuts, we will next argue that the second and third diagrams of (4.11) vanish. For the sake of this argument, we will now consider all the particles of the amplitude in the r.h.s. of the cut as outgoing. Eq. (4.4) involves a physical {in} → {out} amplitude, but as customary in the scattering amplitudes literature, we will cross all particles of the amplitude to be outgoing. Consider first the second diagram of (4.11): the only potential contribution to such cut involves a tree-level amplitude with all negative helicity (outgoing) particles, which vanishes in the SM because such maximal helicity violating amplitude does not exist. Next we consider the third diagram in (4.11), which for our case is To obtain the l.h.s. scattering amplitude, one must sum over all the diagrams with the gauge boson attached anywhere on the particles to the left of the cut. If the gauge boson is attached to the lepton fermion line, the internal Higgs line goes on shell and the diagram factorizes into a one-loop dressing of the Higgs line and an all negative helicity tree-level diagram that vanishes. If 7 We emphasise that in (4.16) the particles are not outgoing, instead (4.16) describes the actual physical process x a y b z → 32. For crossing fermions we used the same rules as in [3], which can be derived by gluing factorized amplitudes that exchange a fermion, see also appendix of [5]. 8 The hypercharge of the left-handed fermions has been expressed in terms of those of the right-handed fermions and the Higgs Y ψ R + Y ψ L + Y H = 0. 9 Note that for the choice of momenta assignments in (4.14) the 12 [12] factor in (4.12) should now be 32 [32].
the gauge boson is attached elsewhere (i.e. Higgs lines or top loop), the sum of the diagrams is equal to 3x f (2, y), with f a function that depends only on the 2 and y spinors. This amplitude vanishes because it is not possible to build an invariant with 2y and [2y] which has zero helicity for the scalar (of momenta p y ) and helicity h = −1 for the gauge boson (of momenta p 2 ).
Therefore the contribution of O eµ L to O eµ DB proportional to N c y 2 t comes only from the threeparticle cut (4.14) and is given by −N c y 2 t (Y µ R + 2Y H ) /(16π 2 ) 2 . This particular contribution is actually zero for the SM hypercharges Y µ R + 2Y H = 0. This is a surprising accidental zero in the anomalous dimension matrix of the SM EFT operators. However, below we will show that other mixings of the typeψγψHDH 2-loop −→ dipoles are not zero. The computation of the mixings from O eµ R to the dipoles is quite analogous: only three-particle cuts matter (for exactly the same reasons as before). The only cut is given by Again, one must sum over all the diagrams where the gauge boson is attached anywhere on the l.h.s., so that we are led to the following five-point amplitude while the form factor is It is again straightforward to perform the rotation (4.13) and integrations in (4.12). 10 The SU(2) tensors A and B for the C eµ R contributions are given in the table of (4.24) below. We get It is now straightforward to generalize our computation for the remaining mixings from the operators O L , O L3 and O R to the dipoles. There are a few minimal changes that we now specify. First, we need to determine the tensors A and B. We find: 10 In this case, given the momenta assignments in (4.20), we have to consider a factor 12 [12].
We remark that the A and B tensors for O L and O L3 should be used in (4.16) and (4.17), while those of O R are used in (4.21) and (4.22). Next, we note that for the O DW computation, apart from using the correct A tensors given in (4.24), we also set the hypercharges Y e R = Y µ R = Y t R = 0, Y H = 1 in the amplitudes (4. 16,4.21), and the SU(2) generators are already included in A.
Finally, putting all the contributions together we get the following anomalous dimension matrix where in the last step we have neglected the terms of order O(y e /y µ ) ≈ 0. The contributions from C µe L , C µe L3 and C µe R into the dipoles can be easily obtained by exchanging µ ↔ e in the amplitudes and form factors. We find (γ C µe DB , γ C µe We note that the contribution of O L and O L3 to O DB is proportional to 2Y e R + Y H = 0, which seems to be an accident.

Higgs quartic λ contributions
Starting with the mixing from O L to the dipoles, we find two three-particle cuts: The amplitude on the first cut (left diagram) is given by where the tensors are given by C ab kc = g · 2δ a k δ b c . The form factor is where D c bla = (−δ c b δ l l ) l a . We are introducing these tensors in anticipation of the generalization to non-abelian SU(2) dipoles. It is now straightforward to compute the three-particle cut. We have to compute 42 [42] where s ij = (p i + p j ) 2 . Regarding the second cut (right diagram of (4.27)), we need the following amplitude where F b akc = 2(δ b k ac + δ b c ak ), and the tree-level form factor where G ac lb = g (−δ a l δ c b ). Once more, plugging this amplitude and form factor into (4.12) (there is no 1/2! symmetry factor to be included for this cut), performing the rotations and a few simple integrals, we get we find that the logs nicely cancel. Therefore the contribution of O eµ L to O eµ DB proportional to λ is given by 6λY H /(16π 2 ) 2 .
Next we provide the details of the O eµ R → O eµ DB mixing. We need to consider two three-particle cuts: The amplitude of the first cut is given by (4.28). The form factor is where D c bla = δ c b δ l l l a . For the second cut, the amplitude is where F b lkc = 2(δ b k lc + δ b c lk ). Note that there is a sign in (4.37) w.r.t. (4.31) due to crossing; this sign is crucial to get a local result, i.e. that the logs like (4.30) and (4.33) cancel. The form factor is given by where G c b = g δ c b . Finally we generalize our computations for the rest ofψγψHDH → dipole mixings. It is now a simple matter to do so because the various computations differ only on the C × D and F × G tensors, which are given in Appendix C. All in all we find that the quartic contribution to (γ C eµ DB , γ C eµ DW ) T = γ D3 · (C eµ L , C eµ L3 , C eµ R ) T is given by

Finite one-loop contributions at the electroweak scale
Following the EFT approach we have to integrate out the W , Z, h and top at the electroweak scale ∼ m W , and match with the Wilson coefficients of the EFT of photons and light fermions. In this process extra finite contributions to d µe,eµ may be generated, as can be found in [22]. We are interested in one-loop corrections arising from C L,L3,R that can compete with our two-loop calculation to the anomalous dimension of C DW,DB . These are one-loop diagrams involving a W or a Z (the Higgs contributions are suppressed by extra Yukawa couplings) in which C L,L3,R enters via the vertex corrections (2.8). From the W we get (neglecting terms proportional to m e ) ∆d eµ (m W ) = e 32π 2 while from the Z we get Notice that since s 2 θ W ≈ 1/4, these contributions roughly cancel for BSM theories that generate C eµ L = C eµ L3 = 0, as occurs in certain models that we will discuss later.

µ → eee in the SM EFT
The process µ → eee arises from the Lagrangian terms apart from the dipoles d µe and d eµ , that generate the branching ratio [42] BR (µ → eee) = 2 The Wilson coefficients entering into the g i at the tree level are only enter in the combination (C µe L + C µe L3 ), as they induce µ → 3e through the Zµe coupling of Eq. (2.8).
At the loop level, other Wilson coefficients can mix into the ones in Eq. (5.2). In particular, four-fermion h = 0 interactions involving other families or the combination (C µe L − C µe L3 ) can enter into the RGE of the Wilson coefficients of Eq. (5.2) at the one-loop level via gauge interactions, as they are also n = 4, h = 0 terms and the mixing is allowed by the ∆n |∆h| selection rule. The corresponding RGEs are given in the Appendix A and the bounds obtained will be discussed later in Section 7.
On the other hand, C LuQe , C LeQu , associated to four-fermion |h| = 2 interactions, and C µe y , related to an operator with n = 5 and |h| = 1, cannot mix at the one-loop level with the Wilson coefficients of Eq. (5.2), and can only enter into the µ → 3e observable by mixing into d µe,eµ as we explained for the µ → eγ case.

µN → eN in the SM EFT
The conversion µ → e in nuclei can arise from the four-fermion terms − 4G F √ 2 g u L,V (μ L γ µ e L )(ūγ µ u)+g u R,V (μ R γ µ e R )(ūγ µ u)+g u L,S (μ L e R )(ūu)+g u R,S (μ R e L )(ūu)+(u → d) +h.c. , (6.1) defined at the nuclei scale. Also the dipoles d µe and d eµ can enter into this observable via the photon splitting into quarks. The branching ratio is given by [48] L,S S (p) + g (n) L,S S (n) 2 R,S S (n) 2 , where ω capture is the nuclear capture rate of the muon and D, V (p,n) , S (p,n) are overlap integrals defined in [48]. We also define We have neglected the contribution from the s quark. At tree level the Wilson coefficients entering into the effective couplings (6.1) are and similarly for the down-sector with u → d, where g u Z = ( 1 2 − 4 3 s 2 θ W ) and g d Z = (− 1 2 + 2 3 s 2 θ W ). The loop mixing into the Wilson coefficients of Eq. (6.4) follows the same pattern as for the µ → 3e case, as the main difference is the replacement ee → uu, dd. The only new ingredient is the presence of the |h| = 2 operators with Wilson coefficients C µeuu LeQu that can receive one-loop corrections from C µeuu LuQe (and similarly for u → d). These coefficients however can also receive large corrections at the QCD scale [49] and will not be discussed any further here.

Constraints from anomalous dimension mixings
In Table 2 we present the bounds obtained when anomalous dimension mixing is considered. In blue we show the bounds coming from one-loop mixings into the Wilson coefficients of the observables, while in red are those in which the mixing is at the two-loop level. In purple we show the bounds for those Wilson coefficients entering into the observables by a two-step one-loop mixing, an effect of order (3.5). We remark that we are not including here any finite loop contributions, that could be larger but are also very model dependent as there can be cancellations depending on the details of the model (see next section for particular cases).
Let us start considering the bounds obtained from the two-loop mixings into the dipole transitions d eµ,µe that is the novel part of this article. The most interesting bound is on the combination (C µe L − C µe L3 ) which, as we explained above, has no serious constraint at the tree level. Our two-loop calculation of the mixing effect in Eq. (4.25) and Eq. (4.26) shows that one can get a bound from µ → eγ of order Λ/ C µe L − C µe L3 24 TeV (assuming a running from that scale). This is quite competitive as compared with bounds coming from a one-loop mixing into (C µe L + C µe L3 ) where this latter is strongly constrained from µ → 3e and µN → eN . The two-loop mixing effect leads to a bound only a factor ∼ 2 smaller than that coming from the one-loop mixings. This provides an interesting correlation between these 3 observables, in the sense that if one of them is measured in the near future, the other 2 should also be experimentally accessible. On the other hand, the bound derived on (C µe L + C µe L3 ) from the two-loop mixing into µ → eγ is much weaker (in part, because it is only generated from two loops involving the Higgs, Eq. (4.39) and Eq. (4.40)) than that from µ → 3e or µN → eN (see Table 2). This makes it unfeasible to see this particular BSM effect in µ → eγ.
The other interesting bound from a two-loop mixing into µ → eγ is for C µe y . One gets Λ/ C µe y 4 TeV that clearly overcomes the tree-level bound from h → µe. In fact, this result is already telling us that µ → eγ constrains this branching ratio to be BR(h → µe) 2 · 10 −8 , making it inaccessible at the LHC or even at future colliders. 11 Let us now move to one-loop mixing effects. Although this analysis has been previously done in the literature [22][23][24][25]1], we will provide here an understanding of the quality of the bounds from selection rules of operator mixings [11,10]. In particular, the only mixing at the one-loop level into the dipoles (n = 4, |h| = 2 operators) is C µett LuQe whose associated operator has n = 4, |h| = 2. This operator can be induced from leptoquarks. One gets a quite strong bound, Λ/ C µett LuQe 304 TeV. Due to the presence of the Yukawa coupling y t in the RGE (see Eq. (A.4)) these effects are much smaller for other families. By mixing into this Wilson coefficient C µett LuQe , other four-fermion operators can enter into µ → eγ by a two step one-loop mixing, and get a bound only slightly weaker (see purple bounds in Table 2). In this mixing the n = 4, h = 0 operators need to involve a y t Yukawa coupling (as dictated by the only exception to the selection rule ∆n |∆h| [11,10]), and as a consequence the mixing is only relevant for Wilson coefficients involving the top.
Finally, other bounds on µef f operators, where f can be any SM fermion of the 2nd and 3rd family, can arise by mixing to µeee, µeuu, µedd at the one-loop level. This can occur via gauge interactions as all these operators are n = 4, h = 0. As seen in Table 2 (for the 3rd family, but the same applies for the 2nd), the bounds range between 10 − 100 TeV, depending on the hypercharges of the states.

Impact on UV Models
As an application of our EFT analysis we would like to consider the impact of the discussed oneand two-loop anomalous dimensions for concrete BSM scenarios. In particular we will consider models with extra heavy fermions and BSM that violate lepton universality. 11 Bounds on C µlle,µqqe LR can also be obtained by a two-loop mixing into µ → eγ, but the contributions are proportional to y l,d , which leads to weak constraints.

Heavy vector-like fermions
Let us consider a heavy vector-like fermion with mass M that can either be a singlet (S), a hypercharged Y E = −1 state (E), or a SU(2) L doublet (D). We assume that they couple to the SM by mixing with the SM fermions: We would like to calculate their contributions to µ → eγ. Following the EFT approach, we must first integrate out these vector-like states at the scale Λ = M and match these contributions with the Wilson coefficients of (2.2)-(2.5). At tree level, we find and C µe L,R,L3 = (C eµ L,R,L3 ) * , as well as At the one-loop order, these heavy states also contribute to the dipole Wilson coefficients of Eq. (3.3). These can be extracted from the contributions to (g − 2) [50]. We find The coefficients C µe DW,DB for S, E, and C eµ DW,DB for D are both O(y e /y µ ) ≈ 0. Next, we have to evolve these Wilson coefficients from M to the electroweak scale. In this RG evolution the Wilson coefficients (8.2) and (8.3) mix at the two-loop level with C eµ,µe DW,DB , in the manner explained in the previous sections. At the electroweak scale, we must now match the theory to the EFT with only photons and light fermions. At the one-loop level, the dipoles can pick up finite terms ∆d eµ (m W ), which are given by using (4.41)-(4.42) and the tree-level matching coefficients in (8.2). For instance, we obtain ∆d eµ (m W ) = − We should still run these coefficients from m W to m µ , but we will not include these effects here as they can be found elsewhere.
Eq. (8.5) shows that the two-loop RG running can be sizable. For instance, for the singlet model S, setting the Yukawa couplings of the heavy fermions to one y (i) S = 1, the current bound on µ → eγ implies M 43 TeV. In this case the RG contribution accounts for approximately the 20% of the total magnitude of d eµ . For the doublet D model the RG contribution has a slightly larger impact. Setting again the heavy Yukawas y (i) D = 1, we find that M 54 TeV, with the RG contribution being a 25% of the total magnitude of d µe . In general, for exponentially larger values of M the RG contribution would dominate, but for relatively low values of M the importance of the RG contribution is model dependent.

BSM with lepton universality violations
BSM sectors that couple only to the muons of the SM have been recently proposed to explain some experimental discrepancies in the muon sector. A particular possibility is the operator which could arise from integrating a heavy vector boson that only couples to muons and to the i family quarks. In the presence of this lepton universality breaking from some BSM, lepton number is not anymore automatically preserved since the diagonalization of the SM Yukawa matrix y e leads, in the presence of Eq. (8.6), to muon number violations. In particular, the operator O µett LL3 is induced with where U L L ,Q L is the left-handed rotation that diagonalizes y e and y u . If y e,u are roughly symmetric, we can estimate U 21 L L ∼ m e /m µ and U Q L ∼ V CKM . Considering the constraint on the Wilson coefficient C µett LL3 from µ → eγ, we obtain the bound On the other hand, the operator (8.6) also induces a contribution to b → sµµ of order C µµbs According to Ref. [51], we need C µµbs LL3 /Λ 2 ∼ 1/(56 TeV) 2 in order to explain the experimental discrepancy in B → Kµµ. From Eq. (8.9) we see that µ → eγ allows for the possibility i = 2 but not i = 3 (unless of course some of the above assumptions are relaxed).

Conclusions
In this work we have analyzed the impact of Lepton Flavor Violation processes with ∆L µ = ∆L e = 1 on the SM EFT. The most stringent constraints arise from µ → eγ, µ → eee and the transition rate µN → eN , where a rich program of measurements is planned in the upcoming decade as summarized in Table 1. Given these spectacular prospects, our main goal here has been to understand at which loop order the different dimension-six operators of the SM EFT mix into these LFV observables. In particular, we have argued that the current and future precision reach of µ → eγ required the knowledge of operator mixings at the two-loop level.
We have shown that due to selection rules only one type of operators enter at the one-loop level into µ → eγ, and only two other types are doing it at the two-loop order. 12 This is sketched in Figure 1. The two operators mixing into µ → eγ at the two-loop order are |H| 2 Hψψ and H † D µ Hψγ µ ψ which at tree level induce LFV h, Z and W couplings. While the mixing from |H| 2 Hψψ was already calculated in [44,4], we have presented here the calculation of the H † D µ Hψγ µ ψ mixing. In particular, we have calculated the two-loop anomalous dimensions of C µe,eµ DW,DB arising from C µe L,L3,R . Our task was greatly simplified by using on-shell tools. In Section 4 we provided a lightning review of the on-shell methods that we used. Then, we explained in detail the calculation of the two-loop mixings, showing how this simply reduces to a product of tree-level amplitudes integrated over some phase space. We have also analyzed the operator mixing for the µ → eee and µN → eN observables, although in these cases we have shown that a one-loop analysis was enough.
An interesting application of our analysis has been to obtain a bound on (C µe L − C µe L3 ) and C µe,eµ y from µ → eγ that is competitive with bounds coming from other observables. In particular we have shown that the bound on C µe,eµ y constrains BR(h → µe) to be too small to be detected in future colliders. This interplay between the different bounds arising from the different LFV precision measurements was discussed in Section 7, and the actual bounds were shown in Table 2 where we indicated the loop order of the mixing of each operator into the process of interest.
Finally, we have shown a few illustrative examples of how to use our EFT analysis to understand the different BSM effects inducing µ → eγ. In special, we have considered models with heavy vector-like fermions and compared the different contributions coming from EFT matching at Λ and m W with those from running.
Our main message here has been to show that the next generation of LFV precision measurements will require the knowledge of renormalization effects at higher orders, where on-shell methods have been shown to be extremely suitable not only for performing calculations, but also for understanding the patterns behind them.

A One-loop anomalous dimensions relevant for LFV
In this Appendix we present the one-loop anomalous dimensions of the Wilson coefficients that enter at tree level into the observables µ → eγ, µ → eee and µN → eN . These results can be found, for example, at [43,52]. We are not interested in self-renormalization but only in mixing effects coming from other Wilson coefficients not present at that level. This gives us the leading order at which the coefficients enter into the LFV processes.

A.1 µ → eγ
The Wilson coefficient entering at tree level into µ → eγ is the combination (C DW −C DB ). This can only be renormalized at the one-loop level by operators with |h| = 2. In particular, the orthogonal combination (C DW + C DB ) can mix into (C DW − C DB ). Using we derive this mixing to be 3) The other |h| = 2 operators are Eq. (2.5). Among them, however, O LeQu leads to an amplitude A ∼ le qu antisymmetric under l ↔ e, so it cannot renormalize Eq. (4.2) and Eq. (4.3) that are symmetric under this exchange [5]. Therefore only O LuQe gives a nonzero contribution to the dipoles at the one-loop level (derived in detail in Section 4.1): There are several Wilson coefficients of Eq. (2.5) however that can enter into the renormalization of C LuQe . Knowing these effects allows us to understand which Wilson coefficients can renormalize the dipoles at the two-loop level with a double log. We have (A.6) The first terms of both equations correspond to a |h| = 2 operator, while the other terms correspond to h = 0 operators where the helicity selection rule ∆n |∆h| is violated due to the Higgs interchange that leads to a contribution ∝ y u y e [11,10].

A.2 µ → eee
The Wilson coefficients C µe L,L3 enter in Eq. (5.2) but only in the combination C µe L + C µe L3 . For this reason it is convenient to define C L± = C L ± C L3 , (A.7) as we are only interested in the mixing from C L− into C L+ . From we obtain the mixing d d ln µ C L+ = g 2 16π 2 where we neglected self-renormalization.
For the four-fermion µeee Wilson coefficients entering into Eq. (5.2), we have where we are interested in the case where q, u, d corresponds to the 2nd and 3rd family, since for the 1st family these Wilsons are already highly constrained at tree level by µN → eN . Also we are interested in the projection C L → C L− /2 and C L3 → −C L− /2. The renormalization from C µeτ τ LL,LL3,RR,LR,RL can be obtained from Eq. (A.10) by the replacement q, d → τ , N c → 1, and putting to zero the contribution from u.

A.3 µN → eN
The relevant anomalous dimensions of the Wilson coefficients C µeuu LL,RR,LR,RL of Eq. (6.4) can be obtained from those in Eq. (A.10) by the replacements Y L L → Y Q L and Y e R → Y u R .

B Conventions and minimal form factors
For the computations in Section 4, we work with 2-component Weyl spinors, and take all fermion fields to be right-handed. Then the SM dimension-four Lagrangian is where i, j = 1, 2 upper (lower) indices of the (anti-)fundamental representation of SU (2). All other indices like Lorentz, families etc. are contracted properly but not shown. We define σ µ = (1, σ) andσ µ = (1, − σ). The covariant derivatives are D µ = ∂ µ − ig(τ α /2)A α µ − ig Y B µ , where τ α are Pauli matrices and Y is the hypercharge.
Next we report the minimal form factors in both 4-component Dirac and 2-component Weyl fermions. Remember that they are Fourier transforms of position-space form factors evaluated at zero momentum, as defined in the text. Starting with the dipole operators: where Σ µν = i 2 [γ µ , γ ν ] such that (ψ L Σ µν ψ R ) † = ψ R Σ µν ψ L . Similarly in 2-component Weyl basis we define σ µν = i 2 (σ µ σ ν − σ ν σ µ ). The relation between conjugated fermion doublets in Dirac and Weyl is L i = C · ij · L * L,j where C is acting on Lorentz indices and = iτ 2 acting on SU(2) indices. Next, the current-current operators: = 2 13 [32](−δ k l δ n m ) , (B.5) where to distinguish the lepton flavor, we used letters E and e for electron doublet and singlet respectively, and M and µ for muon.

C Tensors
Here we give the SU(2) tensors used in the main text to compute the anomalous dimensions: while for the O R operator: We emphasize that for O L and O L3 these tensors should be used in (4.28,4.29) and (4.31,4.32), while for O R they should be used in (4.28,4.36) and (4.37,4.38).