Renormalisation-group improved analysis of $\mu\to e$ processes in a systematic effective-field-theory approach

In this article, a complete analysis of the three muonic lepton-flavour violating processes $\mu\to e \gamma$, $\mu\to 3e$ and coherent nuclear $\mu\to e$ conversion is performed in the framework of an effective theory with dimension six operators defined below the electroweak symmetry breaking scale $m_W$. The renormalisation-group evolution of the Wilson coefficients between $m_W$ and the experimental scale is fully taken into account at the leading order in QCD and QED, and explicit analytic and numerical evolution matrices are given. As a result, muonic decay and conversion rates are interpreted as functions of the Wilson coefficients at any scale up to $m_W$. Taking the experimental limits on these processes as input, the phenomenology of the mixing effects is investigated. It is found that a considerable set of Wilson coefficients unbounded in the simplistic tree-level approach are instead severely constrained. In addition, correlations among operators are studied both in the light of current data and future experimental prospects.


Introduction
Processes involving charged lepton flavour violation (LFV) are very promising places to search for physics beyond the Standard Model (BSM). In the SM with massless neutrinos these processes are forbidden, and even including neutrino masses and mixing these processes are suppressed by ratios m 4 ν /m 4 W . Thus, they are by far too small to be measured in any foreseeable experiment. Therefore, any observation of charged LFV would directly prove the existence of BSM physics. In fact, in many BSM scenarios, such processes can have sizeable decay rates which are naturally in the reach of forthcoming experiments.
Among the charged LFV processes, the category of muonic transitions is particularly interesting, providing three processes with complementary sensitivities to new physics (NP) and stringently constrained by experiment: the current limits from the MEG [1,2] and SINDRUM [3,4]  The future experimental prospects for µ → e transitions are also promising. With the upgrade to MEG II [5], the sensitivity of Br(µ → eγ) will reach ∼ 5 × 10 −14 . Furthermore, Mu3e will improve the sensitivity on µ → 3e by up to 4 orders of magnitude [6]. Concerning µ → e conversion in nuclei, the DeeMe experiment aims at an accuracy of 10 −14 [7], while Mu2e at FNAL and COMET at J-PARC [8][9][10] aim to improve the sensitivity to the conversion rate by four orders of magnitude compared to SINDRUM II. Moreover, the PRISM/PRIME project [11] aims to reach the remarkable limit of Br(µ → e) < ∼ 10 −18 .
LFV processes have been studied in detail in many specific extensions of the SM (for a recent review see for example [12] and references therein). However, in this article an effective-field-theory description of NP interactions is adopted. In this context, Kuno and Okada [13] reviewed µ → e flavour-changing processes and experiments, and the operator basis required to parameterise them.
Constraints on Wilson coefficients, usually at the scale of the experiments, have been compiled for LFV 2-lepton-2-quark operators [14,15] with a tau-lepton [16][17][18], and 4lepton operators [19]. In the quark flavour sector, a long-time effort allowed to establish the QCD running of Wilson coefficients (see [20] and references therein).
However, given that the electromagnetic coupling is much smaller than the strong one, α e α s , QED running is often neglected for the leptons. Czarnecki and Jankowski [21] proved that the self-renormalisation of the µ → eγ dipole reduces the coefficient at low energy. Also in the case of (g − 2) µ (lepton-flavour conserving effective interactions), the QED running and mixing among a basis of SM-induced operators has been included [22]. More recently, the one-loop running of an SU (2)-invariant operator basis [23,24] was presented in [25,26] and loop effects in the SU (2)-invariant theory were calculated [27,28] 1 . The results of [25,26,28] were used in [33], to translate the current µ → eγ bound from the experimental scale to the new-physics scale, using QED×QCD invariant operators below m W , and the SU (2)-invariant operators above 2 .
The purpose of this article is to give the renormalisation-group evolution (RGE) between the electroweak-symmetry-breaking (EWSB) scale and the scale at which the experiments are performed. We consider the Wilson coefficients that are relevant for µ → eγ, µ → 3e and coherent µ → e conversion and include the lowest non-vanishing order in QED and QCD (≡ leading order). Our analysis extends [33] by including the one-loop RGE of vector operators, and the two-loop mixing of vectors to the dipoles, as well as the inclusion of µ → e conversion and µ → 3e in the phenomenological analysis. In addition, we include the dimension-seven lepton-gluon operator that is relevant in µ → e conversion [35].
The outline of the paper is as follows: Section 2 introduces the QCD×QED invariant operators of our effective Lagrangian. The essential formulas for the transition rates of our three processes in terms of the Wilson coefficients of the operators are recalled in Section 3. Then, the RGE of these Wilson coefficients is discussed in Section 4 with the analytic formulas collected in Appendix A. These results are combined in Section 5 to obtain limits on various Wilson coefficients and discuss the complementarity of the three processes.

Low-energy Lagrangian
Following the Appelquist-Carazzone theorem [36], we consider an effective Lagrangian that is valid below some scale Λ with m W ≥ Λ m b . Therefore, it consists of all operators that are invariant under U (1) QED × SU (3) QCD and contain the fermion fields f ∈ {u, d, c, s, b, e, µ, τ }, as well as the QED and QCD gauge fields. In addition to the dimension-four QED and QCD Lagrangians, it contains higher-dimensional operators multiplied by dimensionless Wilson coefficients C. Having µ → e transitions in mind, we restrict ourselves to operators that induce such transitions and are flavour diagonal with respect to the other fields. Concretely, our Lagrangian reads with the explicit form of the operators given by Concerning the dipole operators O D L and O D R , the fields only add up to dimension 5. Thus, they cannot mix into the four-fermion operators (due to renormalisability arguments), although the latter can mix into the dipole operators. Since the dipole operators flip chirality we have defined them by including a factor of m µ , enhancing their dimensionality up to 6. Also, the prefactor e included in O D L and O D R is introduced for convenience (when considering the RGE).
Following [35], the only dimension 7 operators that we include in Eq. (2.8) are O L gg and O R gg . These operators are phenomenologically relevant because, as shown in [37], they encode the effects of scalar operators with heavy quarks (i.e. c, b) below the heavy-quark mass scale m Q , after the matching has been performed. In practical terms, these operators are suppressed by 1/(Λ 2 m Q ) rather than 1/(Λ 3 ). The normalisation of these operators has been chosen such that their Wilson coefficients do not run under QCD at one loop and the factor G F is included to resize the dimensionality down to 6.
In the scenario where BSM physics is realised at a scale Λ < m W , this NP directly gives rise to the higher-dimensional operators in L eff . If BSM physics is beyond the EWSB scale, as it is usually assumed, SU (2) invariant higher-dimensional operators are first generated in the Standard Model Effective Field Theory (SMEFT). Then, the higher-dimensional operators in L eff stem from the matching of the SMEFT to our theory, as performed at tree level in [33]. It is clear that L eff cannot be a valid description of nature above the EWSB scale, as it does not respect the SU (2) symmetry.
However, we stress the fact that if NP is realised not too far above the EWSB scale, then a reduced hierarchy between the NP and EWSB scales would not allow for potential large logarithms from the SMEFT RGE, while the hierarchy between the EWSB scale (our matching scale) and the muon/nuclear scale is always sufficiently large to give rise to important effects.

LFV muon decays
In this section, the expressions for the processes µ + → e + γ, µ + → e + e − e + , and coherent muon-to-electron conversion in muonic atoms µ − N → e − N in terms of the coefficients C i of the Lagrangian Eq. (2.1) will be presented. We restrict ourselves to the tree-level expressions. Loop effects are only included through the RGE of the Wilson coefficients. Thus, for the expressions that follow the Wilson coefficients are understood to be evaluated at the scale of the process.

µ → eγ
Currently, µ + → e + γ is the LFV muonic process with the most stringent experimental bound [1]. At the tree level, the Lagrangian in Eq. (2.1) results in the branching ratio where Γ µ is the width of the muon. The scale of the process is µ = m µ . Operators other than the dipole will enter this process only through the RGE.

µ → eee
The bounds on µ → 3e are already very precise, but Mu3e will significantly improve them in the coming years [6]. The branching ratio expressed in terms of Wilson coefficients where the interference term with the dipole operator is given by Notice that interference between the four-lepton operators can be neglected because it is suppressed by powers of m e . So in contrast to µ → eγ and coherent µ → e conversion, which only impose two constraints, the upper bound on Br(µ → 3e) sets independent constraints on several four-lepton operators. Through the RGEs, this process is also sensitive to operators involving quarks or other leptons.

µ → e conversion in nuclei
The coherent LFV muonic transition in nuclei is produced by the dipole, the vector, and the scalar operators already at the tree level [38]. Once relativistic and finite nuclear size effects are taken into account for heavy nuclei [38][39][40], the transition amplitudes for the three classes of operators exhibit different sensitivities to the atomic number Z. Hence, in principle different target atoms provide different limits on the coefficients of the involved class of operators [35]. Whether it is feasible to perform measurements for many different elements remains to be seen. The SINDRUM collaboration has presented limits for gold, titanium and lead [4,41,42], but the upcoming experiments mostly concentrate on aluminium.
For this process, the Lagrangian L eff as given in Eq. (2.1) is not directly suitable. Instead, a Lagrangian at the nucleon level containing proton and neutron fields is required. This Lagrangian is obtained in two steps. First, heavy quarks are integrated out. This results in a redefinition of the Wilson coefficient of the gluonic operator [37] with an analogous equation for C R gg . Second, the resulting Lagrangian is matched at a scale of µ n = 1 GeV to an effective Lagrangian at the nucleon level. Following [35] the transition rate Γ N µ→e = Γ(µ − N → e − N ) can then be written as where p and n denote the proton and the neutron, respectively. The effective couplings in Eq. (3.5) can be expressed in terms of our Wilson coefficients as with analogous relations for L ↔ R. The Wilson coefficients in Eqs. (3.6) and (3.7) are to be evaluated at the scale µ n .
The nucleon form factors for vector operators are known from the vector-current conservation, i.e. f (u) V n = 0. Hence, the sum in Eq. (3.6) is in fact only over q = u, d. The calculation of the scalar form factors is more involved. Following [43], the values of the up-and down-quark scalar couplings f (u/d) Sp/n from [44] (based on the two-flavour chiral perturbation theory framework of [45]) are used. The values of the s-quark scalar couplings f (s) Sp/n are borrowed from a lattice calculation [46] 3 . In summary, we use (3.8) The form factor for the gluonic operator can be obtained from a sum rule. In our normalisation we get Sp/n . for gold and aluminium, respectively. The branching ratio is defined as the transition rate, Eq. (3.5), divided by the capture rate. For the latter we use taken from [49].

Renormalisation-group evolution
The operators present in L eff , Eq. (2.1), will give rise to µ → e transitions. Thus, experimental constraints from these processes can be translated into bounds on the coefficients of various operators. However, in the first instance this procedure results in bounds on C(m µ ) or C(µ n ), the Wilson coefficients evaluated at the low scale. In order to gain more insight into possible NP scenarios, it is preferable to extract limits on C(Λ), i.e. on the Wilson coefficients at the high scale. If Λ < m W , this will be directly at the NP scale. If the BSM scale is higher than the EWSB scale, our theory will have to be matched to the SMEFT and, in a second step, the RGE within the SMEFT [25,26] will have to be carried out from m W to the scale, where the BSM theory is to be matched to the SMEFT.
Under the RGE, the various operators in L eff mix among each other. In our analysis, we take into account at least the one-loop anomalous dimensions for all the operators. Since the dipole operator plays a prominent role in all µ → e transitions, we also consider those two-loop effects of direct mixing into C D L and C D R that are "leading", i.e. those with vanishing one-loop contribution.
At one loop, the structure of the anomalous dimension matrix splits into two blocks, the vector operators and all the other operators. The dipole operator Q D L receives contributions from itself [21,28], Q S LL ee , Q S LL µµ and Q T LL hh [33,35], and the scalar and tensor operators also mix among themselves. In addition the dimension 7 gluon operators mixes into the quark scalar operators. On the other side, the vector operators mix among themselves [50], but are not connected to the remaining operators.
Many of the anomalous dimensions we need can be found in the literature. However, we have performed an independent computation of all one-loop anomalous dimensions in our basis. In order to perform such a calculation in an automated way, several openly available tools were used: the described model was implemented in FeynRules v2.3 [51] to obtain consistent Feynman rules 4 ; the FeynArts interface of FeynRules was exploited to produce a model file for the FeynArts v3.9 [52] and FormCalc v9.2 [53,54] packages; the combined packages FeynArts/FormCalc were employed to generate non-integrated amplitudes to be elaborated afterwards with the symbolic manipulation system Form v4.1 [55].
At two loops, also the vector operators mix into dipole operators. Since this mixing is not induced through one-loop diagrams, then the two-loop contributions O(α e ) are the leading order effect. In general, the corresponding anomalous dimensions are regularisationscheme dependent. This issue has received a lot of attention in the context of b → sγ transitions [56,57]. For a physical quantity, this scheme dependence is cancelled by the corresponding scheme dependence of the one-loop finite contributions. In fact, this contribution is of the same order in the coupling, namely O(α e ). Here, we present our two-loop anomalous dimensions for vector operators in the 't Hooft-Veltman (HV) scheme, where the finite one-loop contribution vanishes. They have been extracted from [57] and are given in Appendix A. This can also be done for those scalar operators (O S LR hh and O S RL hh ) that can be interchanged with vector operators by Fierz transformations. We remark that our approach, while phenomenologically useful, it is not completely self-consistent and should not be understood as a replacement of a two-loop precision calculation, but rather as a qualitative indication of a leading order effect.
Assembling the Wilson coefficients in a vector C(µ), the renormalisation-group running can be written in matrix form as where Γ i are the transpose of the anomalous dimension matrices. The analytic results for the Γ ST (scalar-tensor mixing) are shown in Table 1, and for the Γ V (vector mixing) and Γ V D (vector into dipole mixing) in the upper squared block and last row of Table 2, respectively. The complete expressions, keeping small mass ratios and charge factors are collected in Appendix A.
In the evolution, the operators involving bottom quarks, taus and charm quarks are integrated down to the various thresholds. Concerning the two-loop running of the dipole operator, we remark that the threshold effects result in vanishing one-loop contributions in the 't Hooft-Veltman (HV) scheme. Operators involving light quarks are not meaningful any longer for scales below µ n = 1 GeV. At about this scale, our theory should be matched to a nuclear effective theory (following the same approach that we adopted for the µ → e conversion in Section 3.3). Concerning the RGE evolution, we also integrate out the light quark operators at 1 GeV.
In order to illustrate the numerical importance of the various terms, we express the Wilson coefficients at a low scale in terms of the Wilson coefficients at µ = m W . To this with an analogous expression for the L ↔ R interchange. The structure of the solution to the RGE, Eq. (4.1), can then be written as The choice of the low scale in the various parts of C on the l.h.s. of Eq. (4.3) and the evolution from µ n to m µ deserves some comments 5 .
We start by noting that C D L , the first entry of C ST l and the first two entries of C V contribute to either µ → eγ or µ → 3e, Eqs. DV does depend on how the evolution from µ n to m µ is treated, these effects overall are so small that they do not noticeably affect the limits we will present in Section 5. Hence we can safely stop the evolution of the light-quark operators at µ n = 1 GeV, even those that are part of C V (m µ ).
Regarding the remaining two coefficients on the l.h.s. of Eq.  For the presentation of the evolution matrices in the remainder of this section we will set the light-quark masses to zero and stop the evolution of the light-quark operators at µ n . While this will affect some of the entries of Eq. Finally, the largest matrix is given in the form and the entries ofŪ V V (m µ , m W ) are given bȳ Of course, it is no problem to change the high scale from m W to another scale Λ, as long as m W ≥ Λ m b .

Phenomenological analysis
In this section we use the concepts presented previously for a phenomenological analysis 6 . For this purpose, we assume that the Wilson coefficients (which are generated by some underlying NP theory above the EWSB scale) are given at Λ = m W . In a first step, we use the RGEs presented in Section 4 to evolve the Wilson coefficients from the high scale m W to the low experimental scale (namely µ n and m µ , respectively). At the low scale, the predicted rates are then confronted with the experimental limits. This results in constraints on the various Wilson coefficients at the high scale. Due to mixing effects in the RGE, Wilson coefficients that are zero at the high scale can be non-zero at the low scale. Hence we are able to place bounds on coefficients which would be unconstrained if loop effects were disregarded.
Let us start by assuming that at the high scale m W only one Wilson coefficient at a time is non-zero. The corresponding bounds on the coefficients are shown in Table 3 both for the current and for the future experimental limits. Note that for µ → 3e we did not take into account efficiency corrections due to cuts on the transverse momentum applied in the experimental analysis. These corrections are in general small for 4-fermion operators but significantly reduce the sensitivity to the dipole operator. From Table 3 we can infer the following general structure of these limits: • Experimental bounds on the direct µ → eγ transition represent a powerful tool to test the Wilson coefficients of the dipole operator. Furthermore, the impact of mixing effects originating from some scalar and tensor operators can also be examined with high precision. However, future prospects for nuclear conversion are so good, that it could overtake the direct µ → eγ limits. The only (numerically accidental) exception is represented by C S µµ , that will be still better constrained by the next generation of µ → eγ experiments.
• A µ → 3e experiment is the most powerful tool to explore µ-e-e-e Wilson coefficients of four fermion operators, regardless of the Dirac structure of the operator. This is mainly due to the fact that such interactions produce the µ → 3e decay already at the tree level (see Eq. (3.2)) while it enters all other processes only via loop effects.
vector) currents instead: where X ∈ {L, R} and f ∈ {u, c, d, s, b, e, µ, τ }. In a simplistic tree-level approach, µ → e conversion is not sensitive to C XA f f and C XP f f . However, the remarkable outcome is that axial-vector operators mix into vector operators. This results in strong bounds from µ → e conversion in nuclei once the Wilson coefficients are evaluated at a scale higher than the experimental scale. Therefore, the common preconceptions that µ → e conversion is not sensitive to axial-vector currents is not true anymore once loop effects are taken into account. The corresponding results are presented in Table 4. In order to understand how the parity-selection rules work for the vectorial lepton-quark and lepton-tau operators, in Figure 1 we show the Feynman diagrams of their mixing at the one-loop level.
First, vectorial operators do not receive contribution from the diagrams in Figure 1b (they vanish when the wave-function renormalisation is included). After, we remark that the penguin diagram in Figure 1a generates operators in which the flavour conserving current has a vectorial structure, i.e. vector-vector and axial-vector operators. Then, from the contribution of the diagrams in Figure 1c one obtains a "maximally flipped" operatorial mixing: axial-axial into vector-vector, axial-vector into vector-axial, vectoraxial into axial-vector and vector-vector into axial-axial. By combining the contributions, the parity selection rules work as follows: vector-vector and axial-axial operators mix into vector-vector (with the contribution from the penguin diagram), axial-vector mixes into vector-axial, axial-vector and vector-axial mix into axial-vector (with the contribution from the penguin diagram), and vector-vector mixes into axial-axial. These results have been discussed also by previous literature [59,60].
In a next step let us compare the exploring power of current and future µ → eγ, µ → 3e and µ → e conversion experiments by directly relating the branching-ratio limits that are needed for the various processes to achieve a particular bound on a Wilson coefficient. For illustrative purposes, we single out two coefficients namely C S LL µµ and C V RR ee . In Figure 2, the current and future branching ratio for µ → eγ and µ → 3e experiments are compared to the future µN → eN prospects (where N is an aluminium nucleus).
Starting with the upper panel, the horizontal dashed red line for example indicates that a rather modest limit Br(µ → eγ) 10 −12 is as constraining on C S LL µµ as the future Mu3e limit Br(µ → 3e) < 5 × 10 −15 . In order for muon conversion to be more constraining a limit of Br(µN → eN ) < 10 −15 would be required. This is indicated by the vertical dashed red line. The future MEG II experiment will place the strongest limit on C S LL µµ unless the COMET or Mu2e experiments improve their expected limit to reach at least Br(µN → eN ) < 5 × 10 −17 . For this specific operator Mu3e will have less of an impact.  From the lower panel, we infer that the future Mu3e experiment will deliver the best bound on C V RR ee . In order to perform with similar standards, the COMET and Mu2e experiments will have to obtain a limit of Br(µN → eN ) < 2 × 10 −17 . Instead, in this specific context, the MEG II experiment will exhibit a generally weak sensitivity, orders of magnitude below the capability of the other experiments.
So far we have considered a scenario where only one coefficient at a time is nonvanishing at the high scale. Now, let us extend the previous results by assuming two Wilson coefficients are non-vanishing. For this purpose, we generated plots in which the parameter space is analysed in light of current and future experimental limits for all three processes. For a better understanding, we display them in a (pseudo-)logarithmic scale. In Figure 3 we show the allowed regions in the C S LL ee − C V RR ee plane (given at the scale m W ). Comparing current (solid lines) and future (dashed lines) limits on µ → eγ (in green) and µ → e conversion (in blue) to those of µ → 3e (in red) indicate that µ → 3e experiments are most sensitive to these Wilson coefficients. It is also interesting to note that the other experimental setups could be blind to specific regions in parameter space where cancellations occur, while this is never the case for µ → 3e. This is mostly due to the fact that these operators trigger µ → 3e already at the tree-level while they give rise to the other processes only via mixing effects. In this case, µ → e conversion experiments display a superior capability to probe vectorial four-fermion operators, as previously discussed. Even the new Mu3e experiment will be just a little better than the current µ → e conversion limit established by SINDRUM II more than a decade ago. However, the plot demonstrates an interesting complementarity among various experiments: assuming that the underlying theory produces a cancellation both in µ → e conversion and µ → 3e, then µ → eγ experiments will provide a complementary limit, finally closing the allowed region of the parameter space.
In Figure 5, we focus on the coefficients C D L and C V RR ee . Here, the current limit on the C D L coefficient comes from the MEG experiment, and the future one will be set by the , an interesting interplay between the observables implies that all of the future experimental limits are useful to ensure that no blind spots in parameter space exist.
With Figure 6 we conclude our review on correlations, by showing the allowed regions in the C D L − C S LR bb plane. The first information one obtains from the plot is that for this case µ → 3e experiments are less constraining than the other two experimental options. In the long term, µ → e conversion will set the best limits on each Wilson coefficient separately. However, there is a big portion of the parameter space where a cancellation in µ → e conversion occurs. Results from MEG II will play an important role to cover this region.
Of course, our choice of combinations for the free parameters is far from being exhaustive. However, the main message is that the interplay between the various experiments is crucial to cover all corners of the parameter space, also the ones in which cancellation can result in blind spots for one or even two specific experiments. are present at the large scale. A very efficient way to determine the impact of the experimental limits on a particular BSM model is to obtain the Wilson coefficients at the weak scale through matching and then use the RGE. A reasonable approximation for the RGE can be obtained by using the numerical evolution matrices given in Section 4. This determines the coefficients entering Eqs. (3.1), (3.2) and (3.5) and, hence, immediately indicates whether for the chosen parameters the model is still allowed or ruled out.

Conclusions and outlook
In this article, we have provided RGE improved predictions for the three µ → e processes µ → eγ, µ → 3e and µ → e conversion in nuclei. Working within the effective theory valid below the EW breaking scale, we have computed the complete one-loop anomalous dimensions for the contributing dim-6 operators taking into account QED and QCD effects. In addition, we have included the leading two-loop QED effects for the mixing of vector operators into the dipole operators and recalled the formula for the µ → eγ, µ → 3e and µ → e conversion rates.
Our O(α (s) ) RGE is renormalisation-scheme independent and can be used for the evolution of the Wilson coefficients from any matching scale Λ m W to the scale of the experiments, after a O(α 0 (s) ) matching has been performed. If the NP theory is realised well above the EWSB scale, this NP theory first has to be matched to the SMEFT and the RGE within the SMEFT [25,26] has to be applied for the evolution from Λ to m W . If a NP theory is realised below or not too far above the EW scale it is sufficient to match it directly to our Lagrangian and use the RGE discussed in this paper.
In our phenomenological analysis we have provided a numerical solution for the RGE. In a second step, we have summarised the resulting bounds on the Wilson coefficients (given at the scale m W ) under the assumption that only one Wilson coefficient at a time is non-zero (see Tables 3 and 4). Afterwards, we have shown the complementary of the three µ → e processes by pointing out the capability of covering regions of parameter space which would be blind spots for a single process.
The limits presented in this paper should be interpreted in light of the fact that they have been obtained under several simplifying assumptions. In particular, obtaining more accurate predictions for the rates as a function of the Wilson coefficients is not the main aim of including RGE contributions. More importantly, one obtains quantitatively new effects. For example 4-fermion vector operators with b, c or s quarks, which do not enter any of these processes directly, mix into contributing operators resulting in stringent constraints. Furthermore, operators with axial-vector currents, which do not enter µ → e conversion at tree-level, mix into contributing vector operators. Therefore, many more correlations among the µ → e processes are present once the RGE effects are taken into account.
The future prospects for observables involving µ → e transitions are intriguing. MEG II will improve the sensitivity on µ → eγ by nearly an order of magnitude, while the existing bounds on µ → 3e and µ → e conversion could even improve by four orders of magnitude. Interestingly, if µ → e conversion managed to improve further to Br ∼ 10 −18 it could be competitive with µ → 3e even for vector operators involving three electrons once loop effects are taken into account. Furthermore, the search for µ → e transitions in Kaon decays like K → µe or K → πµe (see [61] for a recent account) but also in LFV B (see for example [62,63]) and tau decays (see for example [18,[64][65][66]) can give complementary information. While, in our EFT defined below the EW scale, these processes are completely unrelated this situation would change once flavour symmetries are involved or if EW matching effects are considered.

A Anomalous dimensions
In this appendix, the running of the coefficients of the operators listed in Eqs. (2.2)-(2.8) are presented. We useĊ and C F = (N 2 c − 1)/(2N c ) with N c the number of colours. Q l = −1, Q u = 2/3, Q d = −1/3 are the charges associated to leptons, u-type and d-type quarks, respectively. The corresponding equations for the chirality-flipped operators can be obtained form the label interchange R ↔ L. In our computation, the covariant derivative is defined according to the convention of FeynRules v2.3: D µ φ = ∂ µ φ − ig s G a µ T a φ and D µ φ = ∂ µ φ − ieQ φ A µ φ, where φ is a generic field, G a µ and A µ are the gluon and photon gauge field respectively, T a is the colour matrix and Q φ is the electromagnetic charge associated to the field φ.
The running of the leptonic scalar and tensorial operators is summarised by the following equations:Ċ S LL = 12 α e Q 2 l C S LL for ∈ {e, µ}, Finally, the gluon operator is defined in Eq. (2.8) such that its Wilson coefficient does not run at one loop,Ċ L gg =Ċ R gg = 0 [67].