Completeness and complementarity for μ → eγ, μ → ee¯e\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ e\overline{e}e $$\end{document} and μA → eA

Lepton Flavour Violation (LFV) is New Physics that must occur, but is stringently constrained by experiments searching for μ ↔ e flavour change, such as μ → eγ, μ →ee¯e\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ e\overline{e}e $$\end{document} or μ → e conversion. However, in an Effective Field Theory(EFT) parametrisation, there are many more μ ↔ e operators than restrictive constraints, so determining operator coefficients from data is a remote dream. It is nonetheless interesting to learn about New Physics from data, so this manuscript introduces “observable-vectors” in the space of operator coefficients, which identify at any scale the combination of coefficients probed by the observable. These vectors have an overlap ≳ 10−3 with most of the coefficients, and are used to study whether μ → eγ, μ →ee¯e\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ e\overline{e}e $$\end{document} and μ → e conversion give complementary information about New Physics. The appendix gives updated sensitivities of these processes, (and a subset of τ → ℓ decays), to operator coefficients at the weak scale in the SMEFT and in the EFT below mW.


Introduction
The observed pattern of neutrino masses implies new particles and interactions (New Physics ≡ NP) in the lepton sector (beyond the Standard Model with left-handed neutrinos), so one of the best-motivated challenges in particle physics is to discover what it is. Knowing where to search would be a useful input to this enterprise. One possibility would to look everywhere where NP is not already excluded; another would be to explore regions suggested by motivated models. However, the guidance from models can be ambiguous, because they are legion, and for a given process, many models may predict rates just beyond the current experimental reach. So this manuscript takes the agnostic view that it is interesting to look for Lepton Flavour Violation (LFV) everywhere it is not excluded. The current aim is to identify where that is; once NP is detected, reconstructing it becomes interesting.
This manuscript focuses on NP that changes lepton flavour µ ↔ e (for a review, see e.g. [23]), for simplicity restricted to processes that are otherwise flavour diagonal, such as µ → eēe or µ → eγγ, but not K → µ ± e ∓ . In order to look for LFV "everywhere that it is not excluded" among such processes, it would be interesting to list observables such that: JHEP02(2021)172 2. the observables are complementary, in the sense that they give independent information about the NP, and cannot be predicted one from the other. In particular, we want to avoid searching for a branching ratio that is already excluded by the upper bound on another observable.
Instead of trying to construct such a list, we ask a more pragmatic question: to what degree are µ → eγ, µ → eēe and µ → e conversion on nuclei (µA → eA) complete and complementary? These processes are selected because the current experimental bounds on the branching ratios are restrictive ( 10 −12 ) [1][2][3][4], and experiments under construction aim to improve the sensitivities to 10 −14 → 10 −16 [5][6][7][8][9] or better [10]. Ideally, a complete list of complementary processes should be independent of the theoretical formalism used to establish it -for instance, it should not apply only to some models, or depend on a choice of operator basis in an Effective Field Theory (EFT). Nonetheless, let us start by restricting to "heavy" NP models, where the new particle masses are at a scale Λ NP m W . Such models can be parametrised below Λ NP in an EFT framework, which allows to separate "known physics" (the data and the Standard Model), from the theoretical speculations above Λ NP . So LFV is described by operators constructed out of Standard Model(SM) fields and respecting SM gauge symmetries. 1 In addition, Λ NP is assumed large enough to justify retaining only a few terms in the 1/Λ n NP expansion, which in the case of LFV, starts at dimension six ∝ 1/Λ 2 NP . In this EFT context, a complete list of observables should contain at least as many members as there are operator coefficients -otherwise there can be combinations of coefficients that are not probed (sometimes called "flat directions"). Unfortunately, for µ → e flavour change, there are many flat directions: current data on µ → eγ, µ → eēe and µ → e conversion impose 12 to 14 [11] bounds on the coefficients of the 80→90 operators listed in section 2.1.
A first step, is nonetheless to explore whether µ → eγ, µ → eēe and µ → e conversion are sensitive to the coefficient of each operator in the basis of section 2.1. Such a study is basis-dependent, and corresponds to calculating "one-operator-at-a-time" bounds, or "sensitivities"(the word used in this manuscript): a coefficient smaller than such values cound not have been seen. This is consistent with the original aim of identifying where LFV is not excluded, and the results are given in the tables of appendix A. However, the coefficients can be larger than these sensitivities, by sitting along various flat directions, which are discussed in section 2.5.
Loop effects are the backbone of this discussion. This is because the contact interactions which are induced at tree level in models, may not mediate the processes which are stringently constrained by experiment. Consider, for instance, a NP model that induces a tree-level contact interaction (bγ α b)(eγ α µ). This mediates the decay Υ → eµ [12], and contributes at one-loop to µ → eēe. But µ → eēe will have better sensitivity, both because the bound on the BR is more restrictive, and because the muon lifetime is longer, since it decays via the weak interactions, whereas the Υ decays electromagnetically. SM loop corrections to the contact interactions, are discussed in more detail in section 2.3.

JHEP02(2021)172
because this facilitates matching onto the chiral SMEFT operators at the weak scale, and because in the lepton sector, the electron from muon decays is relativistic, so ≈ chiral, implying that negligeable interference between operators involving e L vs e R . The operators are added to the Lagrangian as , the operator subscript Lor gives the Lorentz structure and chirality of the fermion bilinears, and the superscript ζ gives the flavour indices. Since this manuscript only considers µ → e transitions, all the operators containē and µ (the µ + → e + processes are described by the +h.c.), and the eµ indices are suppressed from the superscript.
The 22 four-lepton operators are: where l ∈ {e, µ, τ }, X, Y ∈ {L, R}, and X = Y . Then there are 50 operators withtwo leptons and a quark bilinear: where q ∈ {u, d, s, c, b}. And finally, there are 18 operators with two leptons, which include the dipoles and operators with two photons or gluons The running of the muon mass in the dipole operators is neglected here. The dimension seven operators O GG,L , O GG,R were included in µ → e conversion in [19], and the other gluon operators will not be further considered here. The µ → eγγ rate due to the various two-photon operators was calculated in [20], and the Crystal Box experiment [21] set the constraint BR(µ → eγγ) ≤ 7.2×10 −11 . The O F F,Y operators also contribute to µ → e conversion [22], and the SINDRUMII search for µAu → eAu currently has the best sensitivity to C F F,L and C F F,R .

Experimental bounds
We now want to relate the coefficients of these operators to experimental decay rates. We restrict to the bounds on µ → eγ, µ → eēe and µA → eA because the current experimental upper bounds BR 10 −12 [1][2][3][4] are restrictive, and will improve by orders of magnitude in coming years [5][6][7][8][9]. Furthermore, the branching ratios compare to a weak decay, so BR 10 −12 probes a new physics scale Λ NP 100 TeV.
The Branching Ratio for µ → eγ [23], and the current experimental bound [1] are so the dipole coefficients at the experimental scale should be inside the circle in coefficient space given by eq. (2.5), that is, should separately satisfy C D,X ≤ 1.05 × 10 −8 . The branching ratio for µ → eēe is [23][24][25] BR(µ → eēe) ≤ 10 −12 where ln mµ me = 5.35, so e 2 (64 ln mµ me − 136) 204.8e 2 18.78. Measuring the polarisation of the muon and the angular distribution of the electrons [24,25], (and even the polarisation of the electrons), could allow to discriminate among these various contributions.
Combining the µ → eēe and µ → eγ bounds in a covariance matrix allows to obtain separate constraints on the dipole and vector coefficients (see [11]). Since the current bound on BR(µ → eγ) is restrictive, this amounts to imposing the bound (2.5) on the dipole coefficients, and then neglecting them in (2.6): where the coefficients are evaluated at the experimental scale. The conversion of a muon to an electron in nuclei is a sensitive probe of µ → e flavour change in the presence of quarks. The µ − is captured into the 1s state of the nucleus, and can then convert to an electron by interacting with the nucleons or electric field of the nucleus. The SINDRUMII experiment at PSI, with a continuous muon beam, searched for µ → e conversion on Titanium and Gold [3,4], setting bounds BR(µA → e + A) 10 −12 . The theoretical rates for (Spin Independent) conversion on many targets are given in [26], and can be written [11]: where the Branching Ratio is normalised to the capture rate µA → ν µ A on the same nucleus, and is expressed in terms of the coefficients {C} of operators constructed with a JHEP02(2021)172 nucleons. Several comments: 1. The coefficient subscript gives the Lorentz structure (V or S), then the chiral projector of the lepton current, but the nucleon current is not chiral because its more useful in the non-relativistic limit to use a scalar(S), pseudoscalar(P ), vector (V ), axial vector (A) and tensor (T ) basis, where the P, A, T components contribute to the Spin Dependent rate [27,28] so are neglected here. The non-chiral coefficients can be written in terms of chiral coefficients as, eg 2. The formula after the last equality is given in the normalisation of [11], where the coefficients of nucleon operators have been assembled in vectors (2.10) (and similarly for C R ), and the overlap integrals of Kitano, Koike and Okada [26] for target A have been assembled in unit-normalised "target vectors" . Since these vectors are misaligned, Titanium and Gold can measure independent combinations of coefficients [11].

Translating µ → e conversion bounds onto quark operators
The nucleon coefficients can be written in terms of the quark coefficients (that we wish to constrain) using the expectation values of quark-currents in the nucleus {G N,q Γ }, defined as, where Q ∈ {c, b}. For the scalars, the first term is the tree contribution of light quarks in the nucleon, the second is the one-loop contribution of heavy quarks to the gluon density (which contributes to the scalar nucleon current), so the explicit gluon contribution given in the last term only contains contributions from m W or above (such as the (tt)(ēµ) operator).  Table 1. This table is taken from [28]. It gives matching coefficients between nucleon and lightquark-flavour-diagonal operators. The parenthese gives the uncertainty in the last figure(s). The G N,q S were obtained via EFT methods [33,34], and an average of lattice results [35] is used for the strange quark. The heavy quark scalar G S s are from [36]. In all cases, the MS quark masses at µ = 2 GeV are taken as m u = 2. This allows to translate the upper bound on the Branching ratio on Gold BR Au to a bound on the quark-level coefficients at a scale ∼ 2 GeV:

JHEP02(2021)172
where v m t , and C bb S,R is at m b rather than 2 GeV. The same bound applies for L ↔ R. The SINDRUMII upper bound on the Branching ratio on Titanium, (see eq. (2.8)), is slightly less sensitive to individual operator coefficients, but constrains a different linear combination. In the formulation of [11], a target nucleus of charge Z can be viewed as unit vectorv Z in the space of operator coefficients, such that BR ∝ |v Z · C| 2 . This allows to write the direction probed by Titanium as the direction probed by Gold plus an orthogonal vector:v T i =v Au cos φ +v ⊥ sin φ where cos φ =v T i ·v Au . So one can solve forv ⊥ , and with a covariance matrix obtain that the SINDRUM bound on Titanium gives a constraint on |v ⊥ · C|: (2.14) and the same bound applies for L ↔ R.

Accounting for loop effects
In a Wilsonian sense, the coefficients in the branching ratios of the previous section are evaluated at a scale Λ low near the experimental scale. In order to estimate sensitivities to JHEP02(2021)172 coefficients at Λ NP , loop corrections on scales Λ low → Λ NP must be "peeled off". These loops corrections are evaluated here in an EFT perspective, where only SM particles are dynamical below Λ NP , and new particle propagators are expanded in p 2 /Λ 2 NP around a contact interaction. 2 These SM loops are then independent of the New Physics, and only need to be calculated once as Renormalisation Group running of operator coefficients, and possibly as matching coefficients as discussed below. As usual with the Renormalisation Group, a one loop calculation generates the (log /16π 2 ) n terms for all n, a 2 loop calculation generates the (log n−1 /(16π 2 ) n ) terms, and so on. Provided that Λ NP is sufficiently large, only a few terms in the 1/Λ 2 NP expansion are required. The log-enhanced loops included here arise from the one-loop Renormalisation Group Equations (RGEs) of QED and QCD, which can be written in matrix form as where the operator coefficients are lined up in the row vector C. The diagonal anomalous dimension matrix Γ s of QCD only renormalises the scalar and tensor operators involving a quark bilinear. Since the QCD coupling α s is large and runs significantly, its one-loop effects are resummed: [40]. Two-loop QCD effects are neglected -although this gives an uncertainty of O(10%) -because QCD is less interesting than QED: QED mixes operators, allowing to transform an operator that is difficult to probe experimentally, into one that is tightly constrained. Finally, Γ is the well-populated anomalous dimension matrix of QED, augmented by two-loop QED mixing 3 of vector operators into the dipole [18].
The RGEs are solved analytically as an expansion in α e , which is approximated not to run; we retain the O(α e log) terms, and some of the O(α 2 e log 2 ) and O(α 2 e log) terms. The aim is to include effects that are 10 −3 . The solution is formally a scale-ordered exponential, which can be expanded in α e by defining dt = dµ/µ, and substituting αs 4π Γ s = −iH 0 , αe 4π Γ T = −iH int , into the solution for the time-translation operator given in chapter 4 of [41]. At O(α e ), this gives: (2.16) where there is no sum on KJ in f KJ Γ e KJ , D s =diag{1, . . . λ −a S,T } describes the diagonal QCD running of scalar or tensor operators involving quarks with λ = αs(µ f ) αs(µ i ) , a T = −4/23, a S = 12/23, a J = 0 for J = S, T , and a d describes the running of Γ e KJ in the case where it includes a running QCD parameter (eg a d = a S for anomalous dimensions mixing quark tensor operators to the dipole). At O(α 2 e ), the scale-ordered exponential gives 10% QCD corrections to the mixing of scalars to the dipole ("Barr-Zee"), so these QCD corrections are neglected in the analytic solutions given later.

JHEP02(2021)172
The particle content of the EFT changes when the scale crosses a particle mass, so the operator basis changes too. At each such threshhold, the Greens functions in the theory above and below must be identical, which allows to match coefficients from the basis above onto below. (Matching at the weak scale is recently discussed in [43].) Frequently, this matching can be performed at tree level in both theories; in some cases, the leading contribution of an operator arises via loop diagrams of the theory above, in which case these are included (hopefully after checking that they are scheme-independent, which is not always the case [42]). So for instance, the one-loop matching of scalar b, c quark operators onto the gluon operators, at the quark mass scale, is included as given in eq. (2.12).
What level of accuracy is required? The aim is to include the dominant contribution of each operator (at Λ NP ) to each process. Since almost every operator contributes to each observable (see the tables of appendix A), the question is whether the largest contributions have been included. This is difficult to demonstrate, because the dominant contribution may not be the lowest order in every perturbative expansion (loops, couplings, scale ratios); an example is the two-loop Barr-Zee diagrams, which give the dominant contribution of flavour-changing Higgs couplings to radiative decays like µ → eγ. However, the selection of terms included here has been checked against various models [44,45]. Notice also that including effects 10 −3 does not mean the calculation is accurate to three figures; rather, the resulting constraints on operator coefficients have uncertainties 10%.
Missing from the results given here, are most strong interaction effects beyond Leading Order, below 2 GeV. For example, quark loop contributions to µ → eγ and µ → eēe are included at scales above 2 GeV, but should be replaced by pion loops or resonances from 2 GeV to m µ . An interesting study in this direction is [46]. The µ → e conversion calculation is also at Leading Order in χPT.

Constraints at m W
This section gives the experimental constraints of eqs. (2.5), (2.7) and (2.8) expressed in terms of operator coefficients at the weak scale.

µ → eγ
The MEG bound on coefficients at m W is where all the coefficients are to be evaluated at m W , but the masses are on-shell (= m ψ (m ψ )). Several comments are in order: 1. At one loop, all the tensor operators, and the scalars involving 3 muons or electrons, can mix to the dipole by closing same-flavour legs and attaching a photon. (However, the 3e scalar coefficient can be neglected, because µ → eēe sets a more stringent bound on this operator.) (a) The τ -tensor is in the second parenthese, the quark tensors are in the fourth (the light quark tensors are neglected, because they contribute at one loop to SI µ → e conversion, unsuppressed by m q /m µ ). QCD causes C qq T,XX to run, and also the quark mass that appears in the anomalous dimension mixing the tensor to the dipole is taken running. These QCD effects are described by λ = α s (2GeV)/α s (m W ) 2.18, and where a T = −4/23, a S = 12/23 are respectively the anomalous dimensions of tensors and scalars/masses in QCD, and λ a T 0.873 and f T D 0.862. In eq. (2.17), the quark masses are at low energy: m q (m q ).
(b) It is interesting that the coefficients of the dipole and tensor operators can be comparable, allowing for unexpected cancellations. Numerically, the bounds of eq. (2.17) is where the chirality of the electron outgoing from these operators must be the same as that of the dipole operator into which they mix -so the combination above mix into C D,X for X = Y .
The C 2loop term in eq. (2.17) is compact but not quite correct: in reality, the lower cuttoff of the logarithm should be m τ for the τ τ operators, m b for the bb operators, 2 GeV for the cc,ss,uu and dd operators, and m µ for the µµ and ee operators.
3. the τ and (heavy) quark scalar operators mix to the tensor at one loop, and therefore the dipole at two-loop. Since the tensor-to-dipole mixing is O(1) for heavy fermions because enhanced by the mass, the scalar-dipole mixing is included for the τ on the JHEP02(2021)172 third parenthese, and for the heavy quarks in the last parenthese. The QCD running is 10%, if the heavy quark masses are taken at low energy m Q (m Q ), so neglected. The lower cutoff of the logarithm for quark operators is approximated as 2 GeV to given a shorter formula, but should be chosen as a function of the quark flavoureg m b for the b quark, and 2 GeV for u, d, s.

µ → eēe
The decay µ → eēe constrains the magnitude of several combinations of coefficients, as given in eq. (2.6). At one-loop, there are no operators that mix into the scalar operators by closing the f legs and attaching a photon to the loop, then attaching the photon to e + e − . These "penguin" diagrams mix the combination into both C ee V,XR and C ee V,XL . The 3e and 3µ vector operators also mix to the 3e vectors via a second penguin diagram which closes a different selection of legs, and the 3e operators are renormalised by photon exchange between the legs. As a result, the constraints from µ → eēe (see eq. (2.7)), combined with the MEG bounds on µ → eγ, give the following bounds on coefficients evaluated at m W : where X ∈ {L, R} and X = Y . The overall logarithm multiplying C ping1,X is an approximation to shorten the formula; more correctly, the bb coefficients in C ping1,X should be multiplied by a logarithm cut off at m b , and the logarithm for the ee or µµ coefficients should end at m µ (the correct logs are implemented in the numerical bounds in the appendix).
(And as usual, there should be hadronic loops to replace the light quarks below 2 GeV.)

µ → e conversion
The SINDRUMII bound on µAu→ e+Au, given in eq. (2.8) for out-going left-handed electrons, can be expressed as the bound on coefficients at m W : where ln = ln(m W /m low ) should be inside the square brackets because m low ∈ {m b , 2 GeV, m µ } depends on the coefficient, the masses are on shell (m ψ (m ψ )) and the dipole contribution is given at the experimental scale (it can be written in terms of coefficients at m W with eq. (2.17)). A similar bound applies with L and R interchanged. The bound from Titanium on an orthogonal combination of coefficients, given in eq. (2.14), can be written with the same conventions as where ln is defined after eq. (2.23) and the masses are on shell (m ψ (m ψ )).

What is not probed?
This section gives an incomplete list of coefficient combinations that are not probed by µ → eγ, µ → eēe, and µ → e conversion. There are 45 µ → e L operator coefficients in section 2.1, upon which the current bounds on µ → eγ, µ → eēe, µAu → e L Au and µTi → e L Ti set 1+3+1+1 = 6 constraints [11]. Including µ → eγγ, which was searched for by Crystal Box [21] would give an additional constraint (on a γγ operator) [20,22], but there remain 40 unconstrained directions in coefficient space. The case of operators that produce an outgoing e R is the same and independent, because the Branching Ratios independently constrain both processes. This section does not give a list of three (or six) dozen flat directions. The subset of flat directions listed here are selected because they are relatively stable under RG running, or "natural", according to the notion introduced it [11] (for cancellations among operator coefficients). It assumes that model parameters at the New Physics scale Λ NP do not know the experimental scale, so coefficients should not cancel against logarithms of scale ratios -because the scale ratio could be chosen by the observer but the coefficient is determined by the theory.
An elegant way to implement this notion of naturalness, is to only admit cancellations which are RG-invariant. For instance, if one restricts to operators containing a quark bilinear, and only considers their one-loop diagonal QCD running, then cancellations among the vector coefficients, or among the scalars, or the tensors, are RG-invariant. This could be generalised to include the one-loop QCD and QED running, by diagonalising the full anomalous dimension matrix. However, the diagonalisation would be difficult, and the eigenvectors would be curious combinations of different Lorentz structures and external legs -as opposed to the usual "intuitive" basis of section 2.1, where the operators correspond to potentially distinguishable interactions of physical particles. Furthermore, the only JHEP02(2021)172 degenerate eigenvalues might be among vector operators involving fermions of the same electric charge (because scalars mix to tensors whose mixing to the dipole is proportional to the fermion mass).
A more pragmatic implementation of this notion of "natural" cancellations is to use the perturbative solution of the RGEs given in eq. (2.16), and allow cancellations among coefficients which 1) run with the same QCD anomalous dimension, and 2) multiply a similar α em ln µ f /µ i factor. The formulae for the constraints on coefficients at the weak scale, given in eqs. (2.17), (2.22), (2.23), and (2.24), are arranged to display these "natural" cancellations, or flat directions: each parenthese of the formulae satisfies the two conditions. So any combination of coefficients that causes a parenthese to vanish is a "natural" flat direction.
An example of flatish directions, or operators which are poorly constrained by µ → eγ, µ → eēe and µ → e conversion are the axial vector operators of the form which contribute to µ → eγ at 2-loop with a mass enhancement (see eq. (2.20)). There should be three combinations of these operators, orthogonal to eq. (2.20), which are "flat" to the order calculated here. (The axial vector operators with f ∈ {u, d, µ} are not in the list, because they contribute at one QED loop to µ → e conversion.) A second example would be any combination of vector operators where the {C f } are chosen orthogonal to C ping,1 (see eq. (2.21)), and also to C 2loop if one wishes large C f .
Reference [11]'s notion of naturalness precludes cancellations among the combinations of coefficients in parentheses in eqs. (2.17), (2.22), (2.23), and (2.24), so it could be interpreted as transforming the single constraint from µ → eγ into a constraint on each of the five parentheses, and the 2 bounds from µ → e conversion into eight. However, "unnatural" cancellations can occur, see for instance eq. (2.19).

Complementarity
The aim of this section is to show that the observables considered here (µ → eγ, µ → eēe and µ → e conversion) give independent information about models. To do this, an alternative basis is proposed for coefficient space. This basis is constructed from observables, spans the subspace they probe, and can be defined at all scales.
In an EFT framework, a model M can be represented by the vector C M (Λ) of operator coefficients it induces, 4 which is scale-dependent but relatively simple to calculate at the New Physics scale Λ NP . An observable can be represented as one or several combinations JHEP02(2021)172 of operators whose coefficients it probes. This also depends on the scale, and is relatively simple to establish at the experimental scale Λ expt -for instance, µ → e L γ probes the operator O D,R (m µ ). However, to obtain the rate for an observable, one must evaluate the matrix element of the operators and integrate phase space. So it is convenient to define observable-vectors v obs (Λ) to live in the vector space of the coefficients, such that the rate for an observable can be written at a common scale Λ. For instance, at the experimental scale, eq. (2.5) implies that the two observable-vectors for µ → eγ point in the dipole directions, so one can choose whereû D,R ,û D,L are unit vectors such that C ·û D,R = C D,R . Some observable-vectors for µ → e conversion at the experimental scale are already given in eq. (2.11).
By assumption, the dynamics below Λ NP is Standard Model, so the model and observable vectors can be translated in scale by the SM RGEs. Since models are legion and observables are few, it might be more efficient to translate the observable vectors to Λ NP , rather than calculating C M (Λ expt ), as is commonly done. Indeed, if the observable-vectors were known as a function of Λ NP , then the predictions of a New Physics model would be simple to obtain from eq. (3.1). Part of the translation v obs (Λ expt ) → v obs (Λ NP ) appears in section 2, where the combination C M (m W ) · v µ→e L γ (m W ) is given by eq. (2.17) with X = R, C M (m W ) · v µAu→e L Au (m W ) appears in eq. (2.23), and the appropriate inner product for the orthogonal constraint from µTi → e L Ti is given in eq. (2.24). Similar formulae apply for outgoing e R , but the focus below is mostly on e L . Matching onto the SMEFT, and running up to Λ NP is left for a later work.
The observable vectors corresponding to µ → eēe can also be constructed. For this, it is convenient to take an idealised theoretical perspective, imagining that the chirality of the four leptons can be observed. Combined with the electron angular distributions [24,25], this allows to define "observable" vectors for µ → eēe at the experimental scale to study complementarity at Λ NP , because we are looking for information about the New Physics. However, as an illustrative exercise, we start at the experimental scale.
Notice that the complementarity of observables is a theoretical question, unrelated to the magnitude of current experimental bounds. In this it differs from the correlation matrix that can be constructed from the experimental constraints, which defines the allowed ellipse in coefficient space; the direction and magnitude of the axes of the ellipse depend on the relative magnitude of the constraints. Complementarity is also model-independent, and unrelated to correlations between observables that could arise in some models. Such correlations would occur if the projections of C M onto the observable vectors have similar scaling with the parameters of model M: = independent of model parmeters .
The notion of "model discriminating power", present in the literature [19,47,48], is related to complementarity, but compares the projection of different observable vectors, onto different model vectors. The discriminating power of observables therefore depends on the models considered, whereas their complementarity is model-independent. At the experimental scale, µ → eγ, µ → eēe and µ → e conversion have different external particles, so naively appear complementary. However, the dipole contributes to all three processes, preventing the observable-vectors from being orthogonal even at the experimental scale.
The dipole also contributes to µ → e conversion, but with a smaller weight (see eq. (2.13)), such v µ R →e L γ (2 GeV) and v µAu→e L Au (2 GeV) are almost orthogonal (separated by an angle of ∼ 88 degrees). So at low energy, an approximately orthogonal basis for the observablevector subspace can be constructed with the two dipoles, and the remaining observablevectors with the dipole subtracted. For instance, which are called u to distinguish them from the observable-vectors v. This situation at low energy is illustrated in figure 1, where low energy is taken to be 2 GeV in order to use quark operators for µ → e conversion. The plot is in polar coordinates in the subspace probed by µ → e L γ, µAu → e L Au and µ → e L e L e L , where the vertical z-axis corresponds to the dipole u µ→e L γ (2 GeV), and the x and y axes respectively to u µ L →e L e L e L (2 GeV) and u µAu→e L Au (2 GeV). This subspace corresponds to a Lagrangian at the experimental scale of

JHEP02(2021)172
where O Au is the combination of quark operators probed on Gold, s θ = sin θ, c φ = cos φ and so on, and the radial coordinate is taken This is the complete low-energy Lagrangian for the considered observables; the only "toy" aspect is that chirality of the electrons is fixed to make the parameter space of plottable dimension. The resulting BRs are: .222 cos θ + 9.08 sin θ sin φ| 2 .
Setting the BRs equal to their experimental limits (given in section 2.2), gives a bound on | C| which is plotted in figure 1 as a function of θ, φ; the coefficient space below the curves is experimentally excluded. The angular dependence of the bounds illustrates that the observables are the complementary at the experimental scale. The BR(µ → e L e L e L ) vanishes at φ = π/2 in the second plot (where θ = π/2) because both the dipole and four-fermion contributions vanish; the BR does not quite vanish at θ = π/2 in the first plot because the subdominant four-fermion contribution is present at φ = π/4. One also sees that µ → e conversion gives the best bound on Λ NP , for comparable coefficients of all three operators. This is because the conversion process is coherent on the nucleus (so is enhanced at large atomic number), and because the scalar quark densities in the nucleon are large (∼ 9, see table 1 -the Branching Ratio due to a vector quark current would be ∼ 10 → 100 smaller). Consider now the complementarity of these observables at m W , m W being a simple, although inadequate, substitute for Λ NP . In their evolution from the experimental scale to the weak scale, the basis vectors, like observable-vectors, will rotate in coefficient space and change length (the non-unitary matrix that evolves coefficients in scale can be written as a rotation, a diagonal matrix, and another rotation). The degree of orthogonality between the vectors can be obtained by taking inner products; at the weak scale, the basis vectors for the processes plotted here (= dipoles, plus observable-vectors with their dipole components subtracted) are still orthogonal to three figures despite their significant rotations (this is understandable; the vector space has ∼ 90 dimensions). But their lengths change; in particular the dipole vectors grow by more than a factor two, due to the large loop contributions shown in eq. (2.17).
The observables nonetheless remain complementary at m W , as shown in figure 2. These plots are similiar to those of figure 1, but differ in that the operator coefficients are at m W (standing in for Λ NP ). The first plot shows Λ NP as a function of the angle θ between the model vector C(m W ), and the direction probed by the dipole u µ→e L γ (m W ) (at φ = π/4, that is a model with identical coefficients in the directionsû µ L →e L e L e L andû µAu→e L Au ). The of the model-vector onto the plane perpendicular to the dipole. Were this figure at Λ NP , rather than m W , and if the three Branching Ratios were observed, it would allow to extract the all the information these rates give about New Physics (via dimension six operators and leading order RGEs), at the New Physics scale, so without the blurring effect of SM loops. As expected, µ → eγ is maximal at θ = 0, π, and vanishes for θ = π/2, and µ → e conversion is larger when the model vector is more orthogonal to the dipole. The BR(µ → e L e L e L ) follows BR(µ → e L γ) due to the enhanced (by RG running) contribution of the dipole operator to µ → e L e L e L ; it does not vanish with C · u µ→e L γ at θ = π/2 due to the four-fermion contribution. Similiarly, BR(µAu → e L Au) vanishes at θ π where there is a cancellation between the negative C · u µ→e L γ and positive C · u µAu→e L Au .

Discussion and summary
Reconstructing New Physics (NP) from data is a dream for many phenomenologists. If New Physics is heavy, then Effective Field Theory can be a tool in pursuing this dream, because JHEP02(2021)172 Figure 2. The bound on Λ NP from current constraints on µ → e L γ (black; vanishes at θ = π/2), µ → e L e L e L (blue; dips atθ, φ ∼ π/2) and µA → e L A (red; vanishes atθ, φ ∼ π), as a function of coefficients at m W . Dotted lines indicate possible bounds from future experimental sensitivities. The operator coefficients are parametrised in spherical coordinates in the subspace probed by the experimental processes: the vertical axis is the direction probed by µ → e L γ (see eq. (2.17)) and the xy plane is spanned by u µAu→e L Au (m W ) and u µ L →e L e L e L (m W ), see eq. (3.5). The first plot is for φ = π/4, and the second for θ = π/2 (where the dipole vanishes). See discussion at the end of section 3.
it allows to separate what is known -the Standard Model(SM) and low energy data about NP -from the NP that is pursued. In particular, it allows theoretical travels in energy scale, from the experimental scale towards the NP scale Λ NP , because the dynamical degrees of freedom are by assumption in the SM. However, in this manuscript, the experimental constraints are translated only as far as m W ; reaching Λ NP is left for later work.
The processes considered here are µ → eγ, µ → eēe, and µ → e conversion, because current experimental constraints are stringent (BR 10 −12 ) and are expected to improve significantly in coming years (→ 10 −16 ). These are low-energy processes, occuring at an energy-transfer ∼ m µ , whose experimental Branching Ratios (BRs) are reviewed in section 2.2. Together, they set ∼ a dozen constraints on µ ↔ e contact interactions at the experimental scale.
We would like to know what these processes can tell us about NP, with as few assumptions about the NP as possible. A first assumption is that the new particles are heavy, with masses at a scale Λ NP m W (possibly m W ). Furthermore, they are assumed to generate some three or four-particle, µ → e flavour-changing contact interactions (a list is given in section 2.1).
With these assumptions, section 2 explores whether, if µ ↔ e flavour-changing NP exists, we should see µ → eγ, µ → eēe and/or µ → eγ? In order to reach the shortdistance NP interaction that mediates a low-energy µ ↔ e transition, the intermediatescale SM loop corrections must be peeled off, for instance using the Renormalisation Group Equations (RGEs) summarised in section 2.3. This scale evolution can transform one µ ↔ e interaction into another, so at short distances/high scales, the experimental constraints apply to lengthly combinations of coefficients. These are given in section 2.4 at the scale m W , and the appendix tabulates the sensitivities to individual operator coefficients (oneat-a-time bounds), both in the basis of section 2.1 and in the SMEFT. These results

JHEP02(2021)172
include one-loop and some two-loop effects, and should give the leading contribution of each operator to each process (see discussion in section 2.3). The tables show that at this accuracy, µ → eγ, µ → eēe and µ → e conversion are sensitive to (almost) all the operators or section 2.1 (the possible exceptions are a few µ-e-gluon-gluon operators -see eq. (2.4) and following discussion). Recall that the operator list represents all QED×QCD-invariant LFV three or four-legged contact interactions, in a chiral representation for fermions.
It is unclear if this is reassuring, because cancellations in a matrix element can occur among different coefficients. Such directions in coefficient space are sometimes refered to as flat directions. There are many such flat directions in the coefficient space for µ → e flavour change, after imposing the dozen bounds from µ → eγ,µ → eēe and µ → e conversion, because there are O(100) operators. Indeed, there are so many flat directions, some one could be motivated to introduce a limit on how much cancellation can be "naturally" allowed. The flat directions are discussed in section 2.5.
However, more interesting that the flat directions, are the directions that the observables do constrain, because these are what we can use to discriminate among models. So section 3 introduces vectors in coefficient space, which correspond to observables (these are patterned on the "target vectors" of [11,28]). At the experimental scale, there can be several observable-vectors for a give process, each selects a combination of coefficients who interefere in the rate: where C M (Λ expt ) is the vector of coefficients corresponding to model M, evaluated at the experimental scale Λ expt . For instance, for µ → eγ, whose BR is given in eq. (2.5), the two observable-vectors at the experimental scale are the unit vectors that select the coefficients of the dipole operators in the vector C M (Λ expt ). The vectors for the remaining observables are given in section 3. The observable-vectors are interesting, because they are scale-dependent, and can be translated to Λ NP with the Renormalisation Group Equations. Then there would be no need to run operator coefficients down to the experimental scale for every model; rather, rates could be computed from coefficients at Λ NP , by dotting them into v obs (Λ NP ) This is left for future work in the lepton sector. In section 3, the observable vectors are used to study the complementarity of experimental processes. The vector(s) corresponding to a given observable indentifies the subspace of coefficients that the observable probes, so the misalignment angle between the observable-vectors of different processes, quantifies the complementarity of the processes. When the vectors evolve in scale, this misalignment angle can grow or shrink, indicating that the observables become more, or less, complementary at high scales. For the purpose of learning about New Physics, clearly it is desirable for observables to be complementary at Λ NP . This manuscipt only reaches the weak scale in RGE evolution, so figure 2 illustrates the complementarity of µ → eγ, µ → eēe and µ → e conversion at the scale m W .

A Tables of sensitivities
The tables in this appendix give the sensitivities of various processes (listed in the first row of the tables) to the operator coefficients given in the left column and evaluated at m W .

JHEP02(2021)172
The "sensitivity" of an (unobserved) process to a coefficient is calculated by allowing only that coefficient to be non-zero at m W . Coefficients smaller than these values are too small to have been observed, but larger coefficients could be allowed, if their contributions cancel against other coefficients.
Tables 3, 4, 8 and 9 contain the QED×QCD-invariant operators relevant below the weak scale, and listed in section 2 for µ → e flavour change The operators are added to the Lagrangian with the normalisation given in eq. (2.1): where the subscript is the Lorentz contraction, and the superscript ζ represents the flavour indices. Notice that in the normalisation used here, all operators annihilate muons (or τ s) and create electrons; the reverse process is assured by the +h.c. So the coefficients of an operator and its conjugate have the same magnitude.
At the weak scale m W , the low energy operator basis can be matched onto the SMEFT, so the sensitivities can be expressed in this basis. Tables 5 to 7 (for µ ↔ e) and 10 (for τ ↔ e) apply to operators in the SMEFT basis, added to the Lagrangian as: where the sum is over all dimension six operators and all flavour indices. In the usual SMEFT formulation, the +h.c. is not included for hermitian operators; here there is +h.c. for all operators but the hermitian operators are defined with a factor 1/2 (which gives the usual normalisation for coefficients, because the operator and its conjugate contribute to the Feynman rule.). The SMEFT convention of summing all flavour indices causes some four-lepton operators to appear several times: for ∈ {e, µ} and l ∈ {e, µ, τ }. So the coefficients of these identical operators are also identical, and the bounds in table 6 apply to the appropriate sum of coefficients:

JHEP02(2021)172
Tables 3, 4, and 5 to 7 give the sensitivity of µ → eγ, µ → eēe and µ → e conversion to µ ↔ e operators. Since the dipole receives loop contributions from almost all operators, it is sensitive to most coefficients. It also contributes to both µ → eēe and µ → e conversion; when the only contribution of a coefficient to these processes is via the dipole, the sensitivity is in parentheses in tables 3 and 4,.
The µ ↔ e sensitivities are given to three figures to mitigate rounding uncertainty; modulo misprints and factors of 2, they are expected to have a ∼ 10% uncertainty due to the one-loop QCD running, and a ∼ 50% uncertainty for scalar, tensor and GG operators in µA → eA, where the lattice and χPT determinations of the scalar quark current in the nucleon differ by ∼ 50% [32,33].
For comparaison, tables 8, 9 give estimated sensitivities of a selection of LFV τ decays to operator coefficients in the QED×QCD-invariant basis. In addition, table 10 gives the sensitivities of some hadronic τ decays to 2 2q operators of SMEFT -this illustrates, in the case of SMEFT operators involving quark doublets, the cancellations that can arise between the u and d quark contributions in τ decays to isotriplet mesons.
Note added. During the completion of this manuscript, appeared a comprehensive study of LFV in hadronic τ decays [49]. The authors calculate decays to a variety of final states using χPT with resonances, and perform some loop matching calculations in order to include the interesting operators O GG,Y [47]. The results in this manuscript are less precise, include τ → γ and τ → 3 but fewer hadronic decays, and include QED loops in the RG running up to the weak scale. (QED effects are interesting, because they can change the external legs and Lorentz structure, transforming difficult-to-detect operators into well-constrained ones.) Husek et al. include the numerically significant QCD running via HEPfit [50], and obtain constraints an SMEFT coefficients, and a correlation matrix.

A.1 Including a selection of tau decays
A few τ decays are included here, to illustrate the differences between LFV involving τ s and the µ ↔ e sector. Firstly, the experimental bounds are more restrictive for muon decays: BR(µ → eX) 10 −12 , to be compared with BR(τ → X) 10 −7 . So there is sensitivity to smaller coefficients in the muon sector. Furthermore, the loop contributions of some operators involving heavy fermion (ψ ∈ {c, τ, b}) bilinears, can be enhanced by the mass ratio m ψ /m µ,τ . This, for instance, enhances the sensitivity to O cc T,LL of µ → eγ with respect to τ → eγ. Both these effects can be seen in the tables.
An advantage of LFV τ decays is the multitude of different hadronic final states, which each constrain a specific combinations of quark operator coefficients. This differs from µ → e conversion, where many quark operators contribute in interference. There are therefore fewer "flat directions" for τ → e than for µ → e. This advantage is exploited in [49], but not here, where only a few decays are considered. Also, only τ → e results are listed; the sensitivities to τ → µ operators can be obtained by rescaling, as given in eq. (A.5).
The limits here are a first attempt to include QED loop effects in some LFV τ decays, allowing to estimate the sensitivity of these processes to operators that contribute via logenhanced loops. Previous works have included a wider range of τ decays, but focus mostly JHEP02(2021)172 on operators that contribute at tree level. The prospects for discriminating among operator coefficients by studying asymmetries and angular distributions was studied in [51,52] for leptonic decays (see also [48,53]) and using hadronic decays was considered in [47]. Black et al. [54] give the sensitivity of various rare decays to a subset of SMEFT coefficients that contribute at tree level; a more complete study in the SMEFT was performed recently in [49]. The intermediate-state contribution of heavy quark mesons to leptonic LFV τ decays is considered in [17].
Only the ∆F = 1 decays such as τ → eēe are included here; the operator basis and anomalous dimension matrices for such processes differ from the µ → e case only by permutation of lepton flavour indices. ∆F = 2 decays (eg τ → eμe) would require additional operators and a dedicated RGE analysis.
Three ingredients are required to calculate the sensitivities tabulated here: the experimental upper bounds on the BRs, the theoretical formulae for the BRs as a function of a complete set of operator coefficients at the experimental scale, and the matrix which accounts for loop contributions by expressing the coefficients at the experimental scale in terms of coefficients at the weak scale. The experimental BRs are in table 2, and the remainder of this appendix gives theoretical formulae for the BRs. The matrix is obtained by solving the RGEs for the LFV operator coefficients (between m W and m τ ), and since we restrict to ∆F = 1 operators, these RGEs are the same as for µ → e decays and conversion [18,45] (with some index changes). A dedicated analysis would be required to obtain reliable sensitivities, and "observable-vectors" for τ -LFV.
The masses of the final state leptons are neglected in the rates. This simplification means that the theoretical BRs and RG evolution are identical for the τ → µ and τ → e sectors, after interchanging µ and e indices. Therefore in the tables, only the four-fermion operators describing τ → e transitions are listed; for operators with two lepton indices, the sensitivities to τ → µ coefficients can be obtained by rescaling: where . . . are the indices corresponding to the final state X, and the Lorentz structure subscripts should be identical for both coefficients. ICI In the case of four-lepton operators, all the µ and e indices should be interchanged, e.g. C µτ ee Lor → C eτ µµ Lor , and so on. The decays τ → µēe and τ → eμµ are mediated by the τ → e and τ → µ operators considered here, but are not included due to temporary discrepancies in the tensor contribution, between my calculation and [51].
The calculation of τ decays to mesons is pedagogically introduced in [54], and a careful study considering many final states has recently appeared [49]. The decays considered JHEP02(2021)172 here are τ → e{π 0 , η, ρ}; the results for τ → µ{π 0 , η, ρ} can be obtained as in eq. (A.5). These mesons are interesting because they probe complementary combinations of operator coefficients at tree level. The decays to π 0 mesons probe axial vector/pseudoscalar operators, in the isospin=1 combination u − d. In the notation of [62], where with f π 92.2 MeV, the Branching Ratio in the presence of (axial) vector operators is where J µ A− = 1 2 (uγ µ γ 5 u − dγ µ γ 5 d), and the coefficient of 2 . RG mixing vanishes for the axial current, but it is renormalised (when attached to a chiral current) so at m W this becomes the "constraint" The correct constraint, and "observable-vector"(s) for this decay, could be obtained from the expression for the BR in terms of all coefficients that can contribute (including e.g. the pseudoscalar operators). However, eq. (A.8) allows to calculate sensitivities, and see cancellations, such as between C uu V,XL and C dd V,XL , due to which the SMEFT operators O EQ and O LQ1 do not contribute to τ → π 0 at tree level.
It is interesting to also include LFV τ decays to the isosinglet η, because there is a contribution from s quarks, and not a cancellation between the us and ds. Still in the notation of [62], with the approximation f η ∼ f π [64,65], one obtains the contribution

JHEP02(2021)172
At m W , this becomes Pseudoscalar operators can also contribute to the decays τ → π 0 , η. The operator expectation values give a contribution of pseudoscalar coefficients to the Branching Ratios of where in the normalisation of eq. (2.1), the coefficients of the operators of eq. (A.10) are . QED loops can mix tensor operators into (pseudo)scalars, so this will give some sensitivity to the u, d, s tensor operators.
Finally, decays to the vector ρ meson are normalised to BR(τ → ρν), assuming ρ → ππ (as in [63]; see [49] for a more sophisticated solution) and with the usual factor of 2 for the normalisation of neutral and charged particles: In the second expression, the contribution of the dipole operator (analogous to the dipole contribution to τ → eēe) is neglected, because the current experimental bounds on τ → eγ and τ → eρ 0 are comparable. QED loops mix vector operators of different quark flavour via penguin diagrams, giving this decay some sensitivity to the coefficients at m W of vector operators with a heavy quark current: