Global Lepton Flavour Violating Constraints on New Physics

We perform a global analysis of the bounds from charged lepton flavour violating observables to new physics. We parametrize generic new physics through the Effective Field Theory formalism and perform global fits beyond the common one-operator-at-a-time analyses to investigate how much present data is able to constrain the full parameter space. We particularly focus on leptonic and semileptonic operators with light quarks, identifying unbounded flat directions, detailing how many are present and which operators are involved. The analysis is performed in the general LEFT formalism, which contains all possible low-energy effective operators relevant for lepton flavour violation, as well as in more restricted scenarios, when operators come from a SMEFT completion. We find that flat directions play no role in the fully leptonic four-fermion operators. Conversely, they significantly hinder the ability to derive global bounds on semileptonic operators, with several flat or at least very poorly constrained directions preventing to fully constrain the parameter space. These results are particularly affected by the proper inclusion of uncertainties in the parameters describing µ − e conversion, which decrease the number of well-constrained directions in operator space when treated as nuisance parameters in the fit. While present data is able to provide global constraints on all operators only in the more restricted scenarios we investigated, very strong correlations among the parameters must exist to avoid conflict with the different observables. We provide correlation matrices approximating our full results as a useful tool to compare present data with particular UV completions.


Introduction
The Effective Field Theory (EFT) formalism constitutes an extremely useful tool to parametrize new physics beyond the Standard Model (SM) of particle physics and to derive modelindependent constrains on its existence.In particular, when the new particles and interactions of the theory beyond the SM are characterized by a mass scale that is not yet achievable in our present searches, it is very natural to integrate these new heavy degrees of freedom out of the theory.Then, their indirect effects at low energies are instead represented by a tower of effective operators of dimension larger than 4 involving the light degrees of freedom and suppressed by corresponding powers of the heavy mass scale through their Wilson coefficients.Among all the observables testing the validity of the SM and probing for new physics, flavour changing processes are extremely suppressed in the SM due to the Glashow-Iliopoulus-Maiani (GIM) mechanism [1] and, as such, provide one of our best windows to explore the physics beyond.This fact is particularly true for charged lepton flavour violating (cLFV) processes, where the GIM cancellation is controlled by the negligibly small neutrino masses.Taken at face value, processes such as µ → eγ or µ − e conversion in nuclei are constraining the mass scale characterizing flavour-changing d = 6 effective operators to Λ > 10 3 TeV [2,3], far beyond the reach of collider searches.
The SMEFT [4,5], that is the EFT that can be built with the SM particle content and respecting its symmetries, is characterized by a very large number of operators, particularly when a general flavour structure allowing for flavour violation is allowed.In this scenario, 2499 parameters are necessary to describe the SMEFT at d = 6 [6], the lowest dimension inducing cLFV transitions.Given the complexity of the problem, one of two common simplifying assumptions is usually made when studying the constraints existing in the EFT operators.The first is to keep only terms linear in the Wilson coefficients of the effective operators, that is, keeping only the interference with the SM contributions [7][8][9].This approach allows to perform comprehensive global fits and identify possible flat directions that may avoid constraints.On the other hand, this approach is not appropriate when dealing with processes that are not possible or very suppressed in the SM, such as charged lepton flavour changing processes.In these cases a common simplification is to consider only one operator at a time [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27].While this allows for a very straightforward derivation of constraints for the corresponding Wilson coefficients, these bounds may be regarded as too aggressive, since they will miss any possible flat directions between the different effective operators contributing to a given observable.
In this work we explore this issue beyond the one-operator-at-a-time approach aiming to investigate what is the present status of the constraints that can be placed on effective operators violating charged lepton flavour.We perform a global analysis of both low-energy EFT (LEFT) and SMEFT frameworks, extending previous results for the LEFT µ−e sector [28,29] and SMEFT τ − ℓ sector [30,31].As a simplifying first step, for the operators that contribute to cLFV with an hadronic component, we do not consider flavour change in the quark sector and we focus only on interactions with the lighter quarks (u, d and s), relevant for the most important observables.Moreover, we will derive constraints on the Wilson coefficients of effective operators directly at the low-energy scale relevant for the observables.
In section 2 we first list the LEFT operators that may contribute to the cLFV observables and discuss how the matching with the SMEFT affects how many independent operators may appear in total.In section 3 we describe the fully leptonic cLFV observables and conclude that, in this case, the one-at-a-time approach would be equivalent to a more elaborate global fit involving several operators.In section 4 we describe the relevant observables constraining semileptonic operators and detail how many independent combinations can be probed, highlighting the importance of nuclear uncertainties when exploring the µ − e sector.In section 5 we discuss if any semileptonic LEFT operators are already constrained by present data when cancellations among different operators are allowed and discuss the existence of flat directions.In section 6 we study how these results are affected by the matching to the SMEFT.In section 7 we make the final simplifying assumption of considering only first generation quarks in the interactions and discuss the conditions under which bounds for all operators may be derived.We summarize our results and present our conclusions in section 8. Finally, we provide the details of our analysis as well as the correlation matrices in the appendices, so our global constrains can be easily incorporated to particular UV-complete scenarios.

EFT framework
In the SMEFT, the SM is extended with a tower of new operators, each of them constructed with the SM particle content and respecting its symmetries, but suppressed with inverse powers of the scale Λ at which the new degrees of freedom are expected, with C n the Wilson coefficients (WCs) associated to each of the dim-n effective operators.At lowest order, cLFV observables are introduced at dimension 6 by dipole, lepton-Higgs and 4-fermion operators, as summarized in Table 1 in the Warsaw basis [5].
On the other hand, the most relevant observables to constrain charged lepton flavour change are decays of light mesons and charged leptons.As such, the appropriate EFT to describe these low-energy processes is the LEFT [32], sometimes also referred to as WEFT or WET, built from the particle content of the SM minus the heavy degrees of freedom (top, Higgs, W and Z) and respecting the unbroken QED and QCD symmetries.The lowest order cLFV operator in the LEFT is the dimension-5 dipole operator, At dimension 6, the operators of interest would be 4-fermion operators containing (at least)  LEFT operators for processes with charged lepton flavour violation (α ̸ = β) extracted from Ref. [32].Also included the d = 5 dipole operator.
two charged leptons, where , with G F the Fermi constant, and the sums run over the quarks q, the lepton chiralities X, Y = L, R and the Lorentz structures x = V, A, S, P, T for quarks and y = V, S for leptons.The lepton flavour indices α, β, γ, δ = e, µ, τ are again assumed to be LFV (at least) in the α-β sector, while the quark sector will be assumed to be flavour diagonal and composed only of u, d and s quarks, as discussed in the introduction.We will also assume all WCs to be real.
The operators of our LEFT basis are listed in Table 2, taken from Ref. [32] but with small modifications 1 .In particular, we choose a basis with definite chirality for the charged leptons, since chirality flips are suppressed by their mass and can be largely neglected, but separating vector from axial and scalar from pseudoscalar structures for the quark bilinears, since in this way it is most straightforward to connect each operator with decays of and into pseudoscalar or vector mesons, as well as with spin independent or dependent contributions to µ − e conversion in nuclei.
When these LEFT operators are matched to the SMEFT, the SU(2) L structure imposes non-trivial correlations and constraints altering the counting of relevant independent operators [32].In particular, the scalar and tensor structures are much simpler, since several of them are not generated by the SMEFT at d = 6 (they are QED invariant, but not U (1) Y invariant), and those that are generated match to a single operator 2 .More precisely, 1 Notice that we consider O V,γδR αβL and O V,γδL αβR = O V,αβR γδL as two different operators, since the former/latter introduces a left/right cLFV current.This will be relevant, for instance, when matching to the SMEFT.
2 Notice that, these relations will hold matching the two theories at the electroweak scale.We neglect the modifications that may come from running down to the low scale relevant for the observables.
Consequently, the scalar and tensorial LEFT sector is much more general than that of the SMEFT.In Section 6, we will make use of this feature to further restrict our global analysis after presenting the status for the LEFT framework.
Conversely, and in contrast to the scalar and tensor operators, the 4-fermion vector and axial operators generally receive contributions from two SMEFT counterparts.Indeed, not only the 4-fermion SMEFT operators coupling two lepton currents are present, but also operators inducing a coupling between a (flavour changing) lepton current and the Z boson after the latter is integrated out.In particular, the combination C (1) Hℓ for left-handed currents and C He for right-handed.All in all, keeping only the leading SMEFT terms, the tree-level matching reads: (1) (1) (1)

C
(1) (1) where, for simplicity and since we are focusing on the flavour-change in the lepton sector, we have set the CKM matrix to the identity.Notice that the effect of the CKM would be to include additional SMEFT operators involving the s and b quarks together with the d providing a (CKM-suppressed) contribution to the d mass eigenstate.Since these operators would induce flavour violation both in the quark and lepton sectors, they would be constrained with additional observables.
From the above relations, we see that 14 different SMEFT operators (although in 13 independent combinations) match into only 12 LEFT operators.This means that low-energy cLFV observables cannot be enough to fully constrain the vectorial sector of the SMEFT.Fortunately, higher energy observables where the Z is not integrated out allow to disentangle among the several SMEFT operators that contribute to a given LEFT one.
In particular, C Hℓ ≡ C (1) Hℓ and C He will induce cLFV decays of an on-shell Z-boson.These decays have been searched for and are strongly constrained by the LHC [33,34] and directly probe the operators of interest [14,22]: where we have defined Since both vertex corrections contribute incoherently, it is possible to extract bounds on each of the WCs, which in turn would allow to disentangle and constrain the SMEFT 4-fermion operators in the r.h.s. of Eqs. ( 7) to (18), provided the corresponding LEFT coefficient in the l.h.s. is bounded.Finally, and for completeness, the matching for the dipole operator is given by: Thus, assuming that LEFT operators originate from the low-energy contribution of the d = 6 SMEFT implies, overall, a significant reduction on the amount of free parameters given the strong restrictions imposed in the scalar and tensor sectors.In Section 5 we will derive general model-independent constraints through the LEFT paradigm and then, in Section 6, we will restrict our study to the subset of LEFT operators generated from the low-energy d = 6 SMEFT.Disentangling among the different SMEFT contributions would simply require to consider the constraints from LFV Z decays on the corresponding operators (see Eq. ( 19)) and then run down and match to the LEFT to combine with the low-energy bounds we will derive in this work and present in the following sections.

Radiative and three-body cLFV decays
Dipole and 4-lepton operators can induce cLFV radiative decays ℓ α → ℓ β γ and three body decays ℓ α → ℓ β ℓ γ lδ .See Table 3 for the limits on these processes.The expressions for their respective branching ratios can be found in the literature, see for instance Refs.[11,12,16,22].Therefore, we will refrain from explicitly showing them, as the crucial observation to be made is that each operator contributes incoherently.This implies that the extraction of bounds on the WCs under the assumption that only one operator is present at a time is equivalent to performing a global fit.We will thus present the bounds on the cLFV LEFT operators that can be constrained with the aforementioned low-energy observables, while also listing the operators that are unconstrained.
The bounds on the WCs are shown in Fig. 1.All cLFV dipole operators can be simultaneosuly constrained by the three radiative decays.This family of operators can also mediate three-body decays and other cLFV processes we will discuss later, such as µ − e conversion in nuclei.However, due to the very stringent bounds imposed by the radiative decays, we will neglect them for the rest of our discussion.
As for cLFV 4-lepton operators, the situation is slightly more involved and depends on the flavour structure of the operators.Considering for instance the Lorentz structure with two left-handed currents and for α ̸ = β, there are 15 independent operators.Out of those, 7 are constrained by three body decays, leaving 8 operators unbounded.The exact same counting holds for the operators with two right-handed currents.In the case of the vector structures Table 3: Set of low-energy cLFV current bounds relevant for our analysis.coupling L and R currents, out of the 33 operators, 18 can be constrained, leaving 15 operators without bound.Finally, for scalar structures 18 out of 36 are bounded and thus 18 remain unconstrained.
The flavour structures of the unconstrained 4-lepton operators involve flavour combinations in which a decay is kinematically forbidden, since the heaviest lepton field involved in the c Figure 1: Current 95% CL upper bounds on 4-lepton and dipole LEFT operators.The same bounds hold for the corresponding (L ←→ R) operators.In the case in which the LEFT is matched as the low-energy realisation of the SMEFT, the bounds on scalar operators do not apply, since they are not generated, but the rest still hold.The bounds shown in the plot are collected in appendix B. operator appears several times, and thus are much less straightforward to probe.In particular, there are 8 such flavour combinations, which can violate flavour in either one or two units, ēµμµ, ēτ τ τ, μτ τ τ, ēµτ τ (∆F = 1) , ēµēµ, ēτ ēτ, μτ μτ, ēτ μτ (∆F = 2) , and come in all the possible quiralities and Lorentz structures3 given in Table 2. Interestingly, the operators with 2 muons and 2 electrons, which simultaneously violate L µ and L e in 2 units, are not entirely unconstrained, since they induce muonium-antimuonium oscillations M µ → Mµ .The rather strong constraints on this process [54] can be used to extract bounds on the corresponding operators [55] (see also Ref. [56] for a recent discussion on bounds for ∆F = 2 operators).Furthermore, operators with ∆F = 1 in Eq. ( 22) will mix with other ∆F = 1 operators for which bounds do exist (such as the μeēe structure).Nevertheless, we do not include this effect as in this work we always consider all WCs directly at the low scale relevant for the observables under study.All in all, no global fit is required for dipole and 4-lepton operators due to their incoherent contributions.Consequently, the LEFT bounds displayed in Fig. 1 can be regarded as global.These bounds also directly apply to the scenario where the LEFT operators arise as the lowenergy operators of the d = 6 SMEFT, which could be then translated to the SMEFT WCs using Eqs.( 4) to (18), with the caveat of those operators which receive contributions from the LFV coupling to the Z.In these instances, the bounds should be combined with LFV Z decays, as outlined above, to independently constrain the different contributions.

Semileptonic cLFV observables
We now shift the focus to cLFV 2ℓ2q operators, which induce, among other processes, µ − e conversion in nuclei, leptonic meson decays and semileptonic tau decays.Contrary to the previous section, the 2ℓ2q operators contribute in different coherent combinations and thus lead to a non-trivial analysis of constrained and flat directions.This will result in potentially different bounds depending on whether the approach is global or assumes only one operator at a time, as we will quantify later in sections 5, 6 and 7.In this section, we start presenting all relevant expressions for our analysis, given in terms of the LEFT operators of Table 2 at the relevant scale for each observable.

cLFV semileptonic τ decays
Operators involving flavour change between a τ and a lighter charged lepton with two quarks induce decays of the τ into the lighter charged lepton and one or two mesons.B factories are ideal to probe for these rare decays and stringent constraints exist for their branching ratios from BaBar and Belle (see Table 3).
Although the number of WCs to probe is high and it involves different chiral structures (see Table 2), the plethora of meson final states listed in Table 3 with available bounds offers great complementarity, allowing to constrain many directions in parameter space.For example, the decay τ → ℓπ 0 provides sensitivity to the axial and pseudoscalar WCs for both chiralities of the charged lepton [57]: whereas the decay τ → ℓπ + π − gives access to the vector/tensor and scalar WCs via the ρ and f 0 resonances, respectively.The computation of the latter requires the input of different combinations.For the c γδLV αβR and c γδRS αβR structures, however, there are additional ones which cannot be related to the above via Fierz identities, explaining the difference in the previous counting.hadronic matrix elements (see e.g.[15,20]).Neglecting interference terms, the parametric dependence of the BR in terms of the WCs is [20]: In addition to the spin of the final meson, different isospin4 structures also allows to constrain different combinations of WCs.For instance, while decays into the ρ-resonance allow to constrain the isovector combination of vector coefficients, the τ → ℓω decay can be used to constrain the isoscalar combination5 instead [57] BR Similarly, the τ → ℓη and τ → ℓη ′ decays [15] constrain the isoscalar and s-quark combinations of pseudoscalar and axial operators, complementing the constrain from τ → ℓπ 0 : Finally, the decay τ → ℓϕ can be used to access the vector and tensor involving s-quarks [57]:

Decay constant
Value Matrix element Value Table 4: Meson decay constants and nuclear matrix elements.Due to the lack of lattice data on the transverse decay constant of the ω-meson, we estimate f T,ω /f ω to be equal to f T,ρ /f ρ .
For the gluon element ãN there is currently no ab initio computation, so only estimations are available.As in Ref. [29], we use the FKS [72] estimate, whose 30% uncertainty is motivated by 1/N c corrections.Notice also the large uncertainty in g s,N S arising from a tension between phenomenological and lattice computations.
where here and in all the above expressions f M (f T,M ) stand for the (transverse) decay constants of the mesons, and the pseudoscalar matrix element h q η ( ′ ) is related to the axial (f q η ( ′ ) ) and gluonic (a η ( ′ ) ) elements through the Ward identity: Their numerical values are given in Table 4.
In analogy to the ρ resonant contribution being included in the full τ → ℓπ + π − decay, the τ → ℓKK process would receive contributions from the vector resonances, particularly from the ϕ.While this decay could provide complementary information to the other observables listed above, as it is sensitive to all isospin combinations [15,20], we will not include it in our analysis since the necessary form factors have large uncertainties.

µ − e conversion in nuclei
Operators containing a muon, an electron and two light quarks can induce µ − e conversion in nuclei, for which very strong bounds exist.In particular, this process has been searched for in four different elements (S, Ti, Pb and Au), as shown in Table 3.
Most studies [10,13,16,25,28] focus on the spin-independent (SI) conversion rate since this contribution is enhanced by the coherent sum of all nucleons and thus grows with the mass number of the nucleus considered.However, this SI conversion rate only receives contributions from vector and scalar operators6 as: where V (p),(n) and S (p),(n) are, respectively, the vector and scalar overlap integrals for the proton and neutron, which have been evaluated in [13,73,74].In our analysis, we will use the latest results of Ref. [74].Moreover, the proton and neutron couplings are given by: where the finite recoil contribution from tensor operators in the scalar combinations has been included [75,76].Even though this contribution is suppressed by m µ /m N , this finite recoil term can provide useful constraints.The nuclear scalar form factors can be evaluated from pion-nucleus scattering [77] and from isospin breaking corrections extracted from the protonneutron mass splitting [78,79].Their values are given in Table 4.
Eq. (31) shows that the overlap integrals play a critical role in defining the constrained and flat directions of WCs and, as such, their uncertainties can have qualitative repercussions on a global fit.This is in stark contrast to the τ − ℓ sector, in which the constrained directions are less prone to uncertainties since they arise from criteria such as isospin or CP.When all uncertainties are neglected, these overlap integrals define independent directions in WC space for each of the 4 nuclei for which bounds are available and thus constrain the 4 operator combinations in Eqs.(33).However, as pointed out in Ref. [80], once the nuclear uncertainties are accounted for (at the 5 − 10% level) some of these directions may become parallel in WC space.This has the devastating effect of reducing the number of constraints that are available, leading to new flat directions that were not present in the analysis without nuclear uncertainties.We will quantify this effect later when performing the global analysis of the µ − e sector.
On the other hand, axial, pseudoscalar and tensor operators can mediate spin-dependent (SD) µ−e conversion in nuclei with spin [17,18,29].While generally subdominant without the coherent enhancement over all nucleons, it is still useful to consider the SD contribution since the stringent experimental bounds on µ − e conversion would still translate into meaningful constraints for these kind of operators.
One of the main limitations of SD conversion versus SI conversion is that the former can only happen in nuclei with spin.Since nuclear pairing combines the nucleons in spinless pairs, nuclei with even mass number are spinless and therefore no bound from SD conversion can be derived.Moreover, the spin structure of heavy nuclei, such as Au or Pb, is challenging to quantify from the theory side.As a consequence, titanium is the only phenomenologically viable element so as to derive limits from SD conversion, having two isotopes with spin: 47 Ti and 49 Ti, with relative abundances of 7.44% and 5.41% respectively.
For an isotope of spin J and keeping only the first multipole 7 , the spin-dependent conversion rate reads [29]: where S τ ij , C τ 0 and C τ 1 are, respectively, the spin structure factors, and the isoscalar and isovector combinations for the transverse (τ = T ) and longitudinal (τ = L) modes.The coefficients are defined as: where the relative sign ± refers to X = L/R and δ ′ and δ ′′ are corrections arising from several nuclear effects [81].We assume the values δ ′ = −0.28(5) and δ ′′ = −0.44(4)from Ref. [29].
The index i = 0, 1 refers to the isoscalar and isovector basis, which is related to the proton and neutron basis as with the proton and neutron axial, pseudoscalar and tensor couplings given by in terms of the the nucleon matrix elements {g A , g 5 , f 1,T } q,N , whose numerical values are given in Table 4.The pseudoscalar matrix elements are related to the axial ones via the Ward identity: being ãN the gluonic matrix element.Finally, the spin structure factors S τ ij , when considering just the first multipole, read: 7 Higher multipoles can be safely neglected, as their contributions are much suppressed.
where the fit functions (u) are given in Ref. [29].Consequently, for each lepton chirality and regardless on how many isotopes we have data on, at most 4 directions can be constrained by SD conversion, namely the the 4 WC combinations in Eqs. ( 35)-( 36 Consequently, the parameter space of axial, pseudoscalar and tensor operators cannot be fully probed with only SD µ−e conversion and complementary probes, such as meson decays, are needed 8 .The 4 independent directions to which SD constraints can be sensitive to can in principle be bounded already through 47 Ti and 49 Ti data.Indeed, Eq. (34) shows that, for a fixed lepton chirality, each isotope constraints 2 directions in WC space, since each of the two modes (T and L) contributes incoherently.Consequently, the two Ti isotopes are able to constrain the 4 WC combinations to which SD conversion is sensitive to.
Summing up, µ − e conversion data would constrain on 8 different coefficient combinations per chirality, 4 from SI and 4 from SD conversion, when nuclear uncertainties are neglected.Assessing how this number is reduced in a more realistic scenario requires a careful and systematic treatment of the nuclear quantities that control the directions and their respective uncertainties.We will investigate this issue in scenarios for which bounds on all Wilson coefficients can, in principle, be derived.In particular, we will compare the results of global fits to the data when uncertainties are properly included as nuisance parameters to the naive fit without them so as to gauge their impact.

Meson decays: M → µe
Fully leptonic cLFV decays of light pseudoscalar mesons π 0 , η and η ′ (see Table 3) provide complementary probes of µ − e cLFV operators involving quarks, as they constrain different WC combinations than µ − e conversion.In particular, this kind of decays are sensitive to axial and pseudoscalar coefficients, accessible only via the subdominant SD contribution that, as we discussed above, does not provide enough information as to constrain them all.
Unlike the τ − ℓ sector, we will not assume isospin symmetry here so as not to artificially lift flat directions, since the scalar and pseudoscalar nuclear matrix elements contributing to µ − e conversion are extracted without that assumption.The expressions for the branching ratios as a function of the WCs are the following [29]: where the matrix elements are given in Table 4.
In principle, bounds on these three decays would allow to constrain 6 different combinations of WCs, given the different relative sign of the pseudoscalar contribution depending on the lepton charges.In practice, however, some of the experimental limits shown in Table 3 translate into extremely weak bounds.Even when considering only one operator at a time, the η and η ′ limits imply bounds larger than O(10) for the WCs of axial operators, since their contributions are chirality-suppressed.For the pseudoscalar operators the situation is marginally better and the constraint on the η decay translates into an O(10 −1 ) bound, while the corresponding η ′ constraint is still marginal.Thus, we will only add to the global fit the decays of the π 0 and the η.
Lastly, notice that none of the previous observables probing cLFV in the µ − e sector are sensitive to the vector operator involving the s-quark.This operator could in principle be constrained by cLFV decays of vector mesons such as the ϕ-meson: Unfortunately, the available constraint on this decay [46] leads to a rather weak bound O(10 2 ) on the corresponding WCs, thus leaving c sV µeX essentially unconstrained.Consequently, we will not consider this bound for the purposes of our global fit.

LEFT analysis
From the expressions for the cLFV observables shown in the previous section, we first consider only one of the LEFT operators in Table 2 at a time and derive upper limits on their WCs (or equivalently lower limits for the scale Λ they probe).We present these bounds in Fig. 2.
As expected, all the constraints in the τ − µ and τ − e sectors are similarly strong, around 10 −3 and 10 −4 for the WCs, since the constraints from BaBar and Belle on the different cLFV decays of the τ are all of the same order (see Table 2).Notably, the bounds on the pseudoscalar operators are slightly stronger since they circumvent the chirality flip suppression otherwise required for the decays into pseudoscalar mesons.In the µ − e sector there is a much larger disparity on the order of magnitude of the constraints depending on which is the observable that dominates the bound.Indeed, the strongest constraints between 10 −7 and 10 −8 correspond to the scalar and vector couplings since these induce the, resonantly enhanced, SI contribution to µ−e conversion.The tensor structures are also strongly bounded at the level of 10 −5 to 10 −6 since they may also mediate the SI conversion via their finite recoil contribution with an additional m µ /m N suppression.Bounds on axial and pseudoscalar operators rather stem from their contribution to the SD conversion, which lacks the coherent enhancement.Nevertheless, pseudoscalar operators are also bounded at the 10 −5 level given the enhancement with m N /m q in Eq. ( 39), while the bounds on the axial structure is order 10 −4 .For all Lorentz structures, the constraints on operators involving the s quark are between one and two orders of magnitude weaker given their suppressed nucleon matrix elements with no bound at all for the vector operator.This justifies the simplified scenario without s quarks studied in Section 7, since flat directions involving s quarks would need to overcome this matrix element suppression.
In this section we go beyond this one-at-a-time approach and attempt to derive global bounds on these LEFT operators discussing, both qualitatively and quantitatively, how the results change when all operators are considered at the same time.For this, and also for the results presented in the next sections, we build the χ 2 function adding the constraints from all the observables discussed in section 4, where we have assumed the best fits of all upper bounds to be at zero and a Gaussian distribution in all cases.Since the parameter spaces describing the different EFTs we will confront with the observables are rather sizable, we will explore them via MCMC sampling.Furthermore, given the existence of both flat directions and of Observable bounds on different operator combinations of very disparate orders of magnitude, we adjust the proposal function of the MCMC sampling to the directions along which we expect the most stringent bounds, based on the different observables, and with matching steps.Otherwise, it would be very easy to miss the extremely thin directions that are poorly bounded or even unconstrained and conclude that stronger bounds applied for all operators.We believe this might explain some differences between our results and others present in the literature.We then derive frequentist confidence intervals for each WC profiling over all others from the χ 2 values obtained after the MCMC exploration of the parameter space.This will be summarized as 95% C.L. bands in the different figures shown throughout the paper.We will also show the different 2-parameter depictions of the frequentist confidence intervals in triangle plots in Appendix A to better showcase the rather strong degeneracies as well as the impact of nuclear uncertainties in µ − e conversion constraints.Finally, we also provide approximate correlation matrices so as to be able to implement our results taking into account properly the very strong correlations found in our analysis.Finally, notice that all of the processes discussed in the previous section constrain simultaneously operators having both chiralities for the lighter lepton involved, since their eventual interference is chirally suppressed.Therefore, in the following we will present our discussion for a fixed chirality, keeping in mind that the qualitative arguments apply to operators of both light lepton chiralities.

τ − ℓ sector
As shown in Table 2, and for a fixed chirality of the lighter lepton ℓ, there are 5 operators per quark flavour, each with a different Lorentz structure (V, A, S, P, T).In our three flavour approach with operators only involving u, d and s, this means a total of 15 indepenndent operators.
All in all, the τ cLFV semileptonic decays under consideration probe 11 different combinations of WCs.Since the LEFT contains 15 operators that can contribute to these processes, there are 4 flat directions (F1, F2, F3, F4) that remain unconstrained: • F3, F4: two combinations of c qA τ ℓX q=u+d,s and c qP τ ℓX q=u+d,s .
The direction F1 can be easily understood since the only constraints on scalar WCs come from τ → ℓπ + π − , which only constrains the isoscalar combination and c sS τ ℓX , leaving the isovector combination unconstrained.The direction F2 arises from the fact that τ → ℓπ 0 is the only probe of isovector axial and pseudoscalar operators, therefore the combination that cancels Eq. ( 24) is left unconstrained.Analogously, F3 and F4 are a consequence that isoscalar and s-quark axial and pseudoscalar operators are only constrained by τ → ℓη and τ → ℓη ′ , thus leaving two unbounded directions along which both Eqs. ( 27) and (28) vanish at the same time.
Consequently, a global analysis is not able to simultaneously constrain all of the coefficients, since combinations along the flat directions are unconstrained.Nevertheless, the flat directions presented above are restricted to axial, scalar and pseudoscalar WCs.Thus, vector and tensor coefficients can all be unambiguously constrained through a global fit in spite of the flat directions.In particular, τ → ℓπ + π − and τ → ℓω respectively constrain the isovector and isoscalar combinations of vector and tensor WCs, while τ → ℓϕ independently constrains the s-quark vector and tensor.Additionally τ → ℓπ + π − also independently constraints the scalar s-quark operator.Summing up, a global analysis can constrain 7 out of the 15 LEFT operators involving flavour change in the τ sector.These constraints, result of the global fit, are shown in Fig. 3.
Finally, it should be noted that adding τ → ℓK + K − would help closing the F1 direction, however we did not include it since there are large uncertainties on the relevant form factor 9 .On the other hand, the flat directions involving axial and pseudoscalar operators are much more difficult to lift.This is because these operators are only involved in the decays with pseudoscalar mesons, and only the decays to π 0 , η and η ′ are available, providing constraints on just 3 independent combinations.

µ − e sector
In complete analogy to the τ − ℓ sector, there are 15 independent LEFT operators involving µ, e and two light quarks.In an idealized scenario, i.e. neglecting nuclear uncertainties, µ − e conversion data can probe 8 different combinations of WCs, as outlined in section 4.2, 4 from SI and 4 from SD µ − e conversion.Additionally, the π 0 and η decays discussed in section 4.3 can provide (weaker) bounds on 2 additional directions 10 .Thus, 15 − 10 = 5 flat directions are present in this scenario.Nevertheless, as for the τ case, some Lorentz combinations can be fully constrained even in presence of these flat directions.Indeed, SI µ − e conversion constrains the scalar and vector couplings of the proton and neutron.For the vector operators this translates directly into the bounds on the couplings to the u and d quarks with the same vector structure shown in Fig. 3.However, there is a flat direction (F1) corresponding to the s-quark vector operator, as none of the observables analyzed in the global fit has any sensitivity to it.In particular, in the SI conversion rate (which is the only observable sensitive to vector operators), the nucleon vector form factor for the s-quark vanishes due to the vector current conservation.
On the other hand, the scalar operators of proton and neutron also constrained from SI µ − e conversion receive contributions from all three scalar operators as well as from the tensor structures through their (suppressed) contribution from the finite recoil term (see Eq. ( 33)).As such, these two constraints are not enough to derive bounds on the scalar or tensor operators.
Regarding the pseudoscalar operators, 2 combinations are constrained from the decays of the π 0 and η → µe.Furthermore, the isoscalar combination of pseudoscalar operators can also be bounded from the isoscalar contributions to SD µ − e conversion (see Eqs. ( 35)-( 36)).Thus, global bounds can be found for all 3 pseudoscalar operators.Conversely, the 9 operators involving scalar, axial and tensor structures are bounded by the 2 constrains from the scalar couplings to the proton and neutron from SI µ − e conversion mentioned above and the 3 remaining independent combinations in Eqs. ( 35)-( 36) from SD µ − e conversion.All in all, 5 constraints for 9 operators, leaving another 4 flat directions unconstrained.Thus, the LEFT µ − e sector has 5 flat directions: • F2, F3, F4, F5: four combinations of c qx µeX q=u,d,s with x = A, S, T .
Summing up, through a global fit and even in absence of nuclear uncertainties, it is only possible to simultaneously constrain the pseudoscalar operators as well as the u and d vector structures.This is depicted in Fig. 3, showing the dramatic effect of performing a global analysis rather than considering just a single operator at a time.It is also interesting to see how the constraints that do survive in this global fit compare to the ones derived through the one-operator-at-a-time approach.While the bounds from the τ − µ and τ − e sectors are only slightly changed, if at all, those on the µ − e sector are relaxed by about 4 orders of magnitude.For the pseudoscalar operators the three necessary constraints come from SD µ − e conversion, π and η → µe decays and, as such, the profiled bound is dominated by the weakest of three, namely η → µe.Conversely, in the one-operator-at-a-time approach the bounds are dominated by the most stringent observable, in this case those from µ − e conversion, explaining the huge relaxation seen in Fig. 3.
For the vector operators, this is instead a consequence of the different operator directions being probed by the SI transition of the different nuclei characterized by their overlap integrals in Eq. ( 33) being very close to each other.The situation becomes even worse when uncertainties on their values are taken into account through the appropriate nuisance parameters.The role of the uncertainties will be fully taken into account and discussed more in detail in the following sections with more constrained scenarios where their impact is more easily investigated.

SMEFT analysis
As we have shown in the previous section, present constraints are not enough to derive bounds on all d = 6 LEFT operators given the several flat directions that remain unconstrained.Thus, it is interesting to consider how the situation changes when the LEFT is obtained as the low-energy description of d = 6 SMEFT operators.As discussed in section 2, the d = 6 SMEFT matching generates non-trivial correlations for the (pseudo)scalar and tensor WCs (see Eqs. ( 4)-( 6)): where the upper (lower) sign refers to X = L(R) operators.These relations reduce the number of independent operators, making them easier to constrain. Observable Table 6: Same as Table 5 but for when the LEFT operators are induced by low energy d = 6 SMEFT, which has less independent WCs (see text for details).In particular, the isoscalar and isovector tensor operator map to the same SMEFT contribution, constrained by both τ → ℓω and τ → ℓπ + π − .

τ − ℓ sector
Upon matching with d = 6 SMEFT, the tensor coefficients for down-type quarks vanish, while the scalar and pseudoscalar coefficients are identical up to a sign.This reduces the number of independent WCs from the 15 in LEFT to just 10 coefficients.On the other hand, these same correlations also reduce the number of independent constraints as some of them become redundant.In particular, since there is no tensor coefficient c dT τ ℓL , both decays τ → ℓπ + π − and τ → ℓω overconstrain the only tensor coefficient c uT τ ℓL .Analogously, the coefficient c sT τ ℓL vanishes, rendering its bound from τ → ℓϕ irrelevant.We display this new counting for the SMEFT schematically in Table 6.All in all, in the low-energy d = 6 SMEFT, 2 of the 11 constraints are redundant, reducing the number of independent constraints to 9.
Compared to the LEFT scenario, Eq. ( 49) implies that the constraint on c sS τ ℓX from τ → ℓπ + π − translates to a bound on c sP τ ℓX as well.Similarly, the bound on the isoscalar combination c uS τ ℓX + c dS τ ℓX translates to a bound on the isovector combination c uP τ ℓX − c dP τ ℓX .The isovector combination of axial WCs is then constrained through the bound on τ → ℓπ 0 .This leaves unconstrained the isoscalar and s-quark axial coefficients, as well as the isovector combination c uS τ ℓX − c dS τ ℓX , for which there are only two remaining constraints coming from decays to the η and η ′ mesons, so only one linear combination of these three operators remains unconstrained.Thus, in this context, the single remaining flat direction will be: The bounds resulting from the global analysis to this framework are shown in Fig. 4. The results seem very similar to those obtained in Fig. 3. Indeed, the bounds on the WCs that were bounded in Fig. 3 are essentially the same in Fig. 4 and no new bars for other WCs appear.This is because, even if the correlations among WCs implied by the d = 6 SMEFT allow to lift 3 out of the 4 flat directions present in the general LEFT scenario, the remaining flat direction involves all the WCs that were not previously bounded.Thus, even though the parameter space is in general much more strongly constraint and 3 flat directions have been lifted, when profiling over all WCs no individual bound is found for the WCs involved in the remaining flat direction.Thus, a better perspective on how constrained the parameter space is in this scenario is provided by Fig. 6 in appendix A through the triangle plot with all different projections for constrained directions in the parameter space.47)- (49).These global bounds, along with their correlations, are collected in appendix A.
Finally, we note again that this flat direction could in principle be lifted by including the τ → ℓK + K − channel which, contrary to τ → ℓπ + π − , is sensitive to all combinations of scalar operators (π + π − is not sensitive to the isovector combination), and thus allows for a constraint on all scalar WCs.This would help closing the whole SMEFT parameter space in the τ − ℓ sector.However, we do not include τ → ℓK + K − since the scalar hadronic form factors suffer from large theoretical uncertainties, as previously discussed.

µ − e sector
Analogously to the previous section, the number of independent operators inducing cLFV in the µ − e sector with light quarks is 10 when matching with d = 6 SMEFT.
We discuss first the simplified (and overly optimistic) scenario of neglecting the impact of the nuclear uncertainties.As discussed before, µ − e conversion data can at most constrain 4 + 4 = 8 combinations between SI and SD, and thus is not enough to cover the whole parameter space.Therefore, it is necessary to consider also meson decay bounds that, although substantially weaker, probe complementary combinations of WCs.This adds to a total of 4 + 4 + 2 = 10 constrained directions, which, at face value, seem enough to constrain all of the 10 WCs.However, as mentioned in the previous section, none of these observables receive contributions from s-quark vector currents, so the corresponding WC (c sV µeX ) remains unbounded, and only 9 Wilson coefficients can be constrained at most.
On the other hand, π 0 → µe is sensitive to the isovector combination of axial coefficients, which is already bounded more strongly by SD µ − e conversion.Indeed, an analysis of the correlation matrix of the 8 WC combinations probed by µ − e conversion reveals that the flat unconstrained direction corresponds to a combination of isoscalar and s-quark axial WCs.This direction is obviously orthogonal to the isovector combination probed by π 0 → µe and thus cannot be lifted by this limit.This means that the constraints derived for axial operators will be dominated by the η → µe bound, which, for axial operators, is very weak due to the chirality suppression.Consequently, even though it is technically possible to simultaneously constrain all these 9 WCs, the constraints derived for axial operators are very weak, ∼ O(10), and do not show up among the other light red bands for the range shown in Fig. 4. On the other hand, it is still possible to obtain meaningful global bounds for the remaining Lorentz structures, even if they are substantially weaker than those derived when considering one operator at a time.
We will now address the effects of adding nuclear uncertainties to the analysis.Accounting for the impact of all uncertainties entering the µ → e transition estimation requires treating as nuisance parameters many different quantities and the analysis can easily get out of hand and become numerically unfeasible.For this reason, we will implement uncertainties in a selection of nuclear quantities that can be particularly harmful, since they may effectively reduce the number of independent constraints.These are: • Nuclear overlap integrals: these define the directions probed by SI µ − e conversion.
The main source of uncertainty comes from the fact that lepton-nucleon interactions are computed at LO in χPT, see Refs.[18,66,82] for discussions about the size of the possible NLO corrections.We will consider 5% and 10% uncertainties for the overlap integrals of light and heavy nuclei, respectively.In particular, we will consider the parameters V (p),(n) and S (p),(n) in Eq. ( 31) as free and independent with Gaussian priors centered around their nominal value and a 5%/10% uncertainty.With this parametrization, there will be a confidence level at which the freedom allowed for the overlap integrals makes two of them parallel [80] and, therefore, redundant instead of complementary at this C.L.This has a dramatic impact in the analysis, since a new flat direction appears, loosening many constraints and changing significantly the correlations among the remaining ones, as depicted in Appendix A.
• Nuclear corrections δ ′ and δ ′′ to the axial contribution of SD µ−e conversion: these have relatively large uncertainties and, when they become equal, a new flat direction arises.This is due to the fact that, in this limit, both longitudinal and transverse modes in Eqs. ( 35)- (36) would become sensitive to the exact same combination of isovector axial and tensor coefficients11 .
• Gluonic matrix element ãN : the naive estimation of this parameter, which would suffer from a ∼ 30% uncertainty [29], leads to a cancellation of the isoscalar pseudoscalar contribution to the SD transition via the relation between the pseudoscalar and axial matrix elements in Eq. ( 40).Thus, taking into account this uncertainty has a sizable impact in the final results of the global fit.
The results of the global fit after including all these nuclear uncertainties are displayed in Fig. 4 with dark red bars.The comparison with the light red bars in absence of uncertainties is remarkable.The greatest impact is seen in the vector and scalar structures whose bounds mainly came from the SI contribution.Indeed, for a sufficiently high confidence level (somewhat below 95% where the bounds are depicted), the directions corresponding to SI µ − e conversion in Pb and S become parallel to those of Au and Ti, respectively.This effectively means that two constrained directions are lost in the fit, and thus µ − e conversion can only constrain 6 combinations out of the total 8 directions probed without nuclear uncertainties.These two lost directions must be then supplemented by the two directions probed by meson decays, which entail weaker constraints.This is best shown by Fig. 7 in appendix A with two very striking features.The first is the rather strong correlations between most observables, which reduces the allowed regions to thin lines in the parameter space.This is the consequence of the extremely different orders of magnitude of the bounds that each operator can set, with more than 8 orders of magnitude difference between those coming from SI µ − e conversion to those from the η → µe decay.The second remarkable feature in Fig. 7 is the abrupt change from the 1 to the 2 σ regions for several parameters.Some correlations are lost and the corresponding allowed regions become much larger when going from 1 to 2 σ C.L. In some parameters there is also a dramatic change in the overall constrain as shown by the abrupt jump of the profiled χ 2 depicted in the diagonal panels.This jump and the dramatic change in the correlations, happen close to the 95% C.L., which is when the nuisance parameters can vary enough so as to make the directions probed by Pb and S as well as Au and Ti parallel, as discussed above.

SMEFT with only first generation quarks
Given that it is not possible to constrain all WCs simultaneously even when considering only operators that match to SMEFT at d = 6, we consider the further simplification of operators involving only the first quark generation.This is motivated by the observables involved when constraining the µ → e transitions, where the dominant contributions are always from operators with u and d quarks.Thus, flat directions involving s quarks may be regarded as particularly fine-tuned, since they generally need to overcome this additional suppression.This is not the case for τ → e and τ → µ transitions, since the s content of several mesons involved in relevant τ decays is significant.Nevertheless, we also show results for the τ − ℓ sector in this simplified scenario for completeness12 .

τ − ℓ sector
As discussed in section 6.1, when all the constraints from the different possible cLFV τ decays are considered, only a single flat direction remains unconstrained in the low energy d = 6 SMEFT.This flat direction corresponds to a combination of the isoscalar axial and isovector scalar operators as well as the axial operator with the s quark, for which only two independent  47)- (49).These global bounds, along with their correlations, are collected in appendix A. constraints are available from τ → ℓη and τ → ℓη ′ .If we assume operators involving only first generation quarks, these two processes are now enough to independently constrain the isoscalar axial and isovector scalar operators and no flat directions remain.Thus, in this restricted scenario, present data are enough to unambiguously constrain all cLFV operators from d = 6 SMEFT involving first generation quarks and we display the resulting bounds from our global fit in Figure 5.

µ − e sector
In the low energy, d = 6 SMEFT with only first generation quarks, there are 7 relevant fourfermion operators.Thus, the dominant 8 different constraints from µ − e conversion are in principle sufficient to simultaneously constrain all of the WCs at hand, obtaining the bounds shown in Fig. 5 with light red color.This is no longer the case when nuclear uncertainties are properly accounted for.As previously discussed, uncertainties effectively reduce the overall number of constrained directions by µ − e conversion data.As such, µ − e conversion data will need to be supplemented by bounds coming from meson decays in order to derive constraints on all operators also in this simplified scenario.
The results of our fit taking into account nuclear uncertainties, are shown in Fig. 5 as dark red bars.The treatment of the nuclear uncertainties as nuisance parameters is equivalent to what is described in section 6.2.We again find that bounds on vector and scalar operators are the most affected, as they can no longer be solely constrained through SI µ − e conversion.However, one important difference is that now the constraints have suffered from a milder degradation with respect to the case of Fig. 4.This is a consequence of the fact that, in this simplified scenario, µ − e conversion can place constraints on 8 − 2 = 6 combinations after nuclear uncertainties are taken into account.Therefore, in order to completely constrain the 7 WCs, only one additional constraint coming from meson decays is necessary, contrary to the scenario with s-quarks, in which two extra meson constraints were necessary.As in the previous section, these constraints are very strongly correlated and are significantly different between the 1 and 2 σ regions, when the nuisance parameters are able to reduce the number of constrained directions, as show in Fig. 8 in appendix A.

Summary and Conclusions
Charged Lepton Flavour Violation is one of our best windows to probe for generic new physics beyond the Standard Model, since the GIM suppression through the tiny neutrino masses makes these searches virtually background-free.The EFT formalism is particularly suitable to study these processes in a model-independent way.Given the large number of free parameters introduced by the EFT formalism, the most common approaches to analize their constraints are either to consider only one operator at a time or to stay at linear order in the WCs so as to explore possible flat directions that may relax the former constraints.However, the contributions of the new cLFV operators to the observables will necessarily be quadratic, as there are no SM contributions to interfere with, and therefore going beyond the one-operatorat-a-time approach in these studies is challenging.
In this work we have focused on the impact of potential flat directions and on studying how the bounds obtained from a one-operator-at-a-time may be affected by them.After introducing the LEFT formalism and its matching to the d = 6 SMEFT, we analized the bounds from cLFV lepton and meson decays as well as from µ − e conversion in nuclei on the WCs directly at the low-energy scale relevant for the most important observables.
We find that, for the dipole and fully leptonic operators, there are no flat directions that may relax the bounds that can be derived simply through the one-operator-at-a-time approach from processes such as ℓ α → ℓ β γ and ℓ α → ℓ β ℓ β ℓ γ , respectively.While there are some operators that are not bounded by these processes and that we list, these do not hinder the constraints derived on the others which are summarized in Fig. 1.
Conversely, for semileptonic 4-fermion operators the situation is very different.It is also very different for cLFV involving the τ with respect to the µ − e case.The former is mainly bounded by searches of τ decays to a lighter lepton and a meson.The different isospin and nature of the meson gives sensitivity to operators with different quark content and Lorentz structure.Nevertheless, and despite the many different and complementary searches, we find that in the most general case there are several flat directions that relax some of the constraints that would be obtained through the one-operator-at-a-time approach.In particular, our LEFT scenario has 15 coefficients, 5 per Lorentz structure (tensor, vector, axial, scalar and pseudoscalar) and 3 per quark type (up, down and strange), that may contribute.While through a one-operator-at-a-time approach all WCs can be bounded at the level of 10 −3 −10 −4 (see Fig. 2), we find that these constraints apply in a fully global analysis for only 7 of the 15 operators (see Fig. 3), corresponding to the vector and tensor structures as well as the scalar one for the s quarks.Among the remaining 8 operators, 4 flat directions exist so that individual bounds on their WCs cannot be derived in a global analysis as one may freely move along the unconstrained flat directions.
Given this situation, we then analize the more constrained scenario of the d = 6 SMEFT operators at low energy.This situation is described with only 10 instead of the 15 parameters of the LEFT scenario.Nevertheless, as shown in Fig. 4, the only global bounds are still those derived for the vector and tensor structures as well as for the scalar operator with the s quark.Indeed, there is a single flat direction involving the other 5 operators remaining.Thus, even though 9 parameter combinations out of the 10 independent WCs are bounded down to ∼ 10 −4 , Fig. 4 does not display a bound for 5 of them since they are all involved in the unconstrained direction.Adding τ → ℓKK would close this remaining flat direction and lead to global constrains for all WCs, nevertheless a better handle on its form factors is needed in order to properly include it in the analysis.The final simplified scenario we consider is when only first generation quarks participate in the observables.In this case, as shown in Fig. 5, all WCs are constrained also in the global analysis with bounds ranging from 10 −3 to 10 −4 .
The µ − e sector for the leptonic operators is the most complex to analyze, since the would-be flat directions are determined by overlap integrals and nuclear parameters defining the particular combination of operators that contribute to the SI and SD transitions of each nuclei.Moreover, uncertainties on these quantities, when incorporated as nuisance parameters to a global fit, alter the constrained directions to the point that some become linearly dependent on others.Indeed, within uncertainties, the directions probed by Pb and S become parallel to those determined by Au and Ti respectively, reducing the 4 independent parameter combinations probed by the different elements to only 2. As such, additional data on more complementary nuclei would be very helpful [80].We notice, however, that SI transitions are fully characterized by only 4 operator combinations: the scalar and vector couplings to protons and neutrons.Thus, once 4 independent operator combinations corresponding to 4 different nuclei have been bounded, new nuclei cannot provide complementary information.The situation is similar for the SD contribution to µ − e conversion.Again only 4 operator combinations (isoscalar and isovetor transverse and longitudinal modes) contribute.Furthermore, present data on µ − e conversion in Ti already provide independent constraints on all 4, so that additional data will not allow to constrain new directions.Conversely, we find that cLFV meson decays such as π 0 → µe or η → µe do provide complementary information, although the bounds are much weaker.Thus, improving these constraints does have a significant impact in global fits when correlations and flat directions are fully accounted for.
With the one-operator-at-a-time approach and neglecting nuisance parameters, very stringent bounds down to 10 −7 are found for the LEFT vector and scalar operators that contribute directly to the SI µ − e conversion in nuclei.Tensor structures have somewhat weaker bounds of 10 −5 through their finite recoil contribution to the SI transition.The axial and pseudoscalar structures contribute to SD transitions instead, which lack the resonant enhancement of the SI and lead to bounds between 10 −3 and 10 −5 (see Fig. 2).None of the observables are sensitive to the vector operator involving the s quark, which remains unbounded.When analizing the LEFT scenario with its 15 free parameters, only bounds on the vector and pseudoscalar operators may be derived.For the vectors they are relaxed from 10 −7 to 10 −3 due to the very many flat directions present (see Fig. 3).The bounds on the pseudoscalar operators are now dominated by the η → µe process and relaxed by around 5 orders of magnitude.
When the low energy d = 6 SMEFT with its 10 operators is considered instead, all of these flat directions are lifted.However, there is a single direction involving the axial operators in an isoscalar combination as well as the s quark contribution which is only very weakly constrained, beyond the range of Fig. 4, by η → µe.The bounds on the other operators are also very degraded with respect to the one-at-a-time constraints since, even when independent, many of the directions constrained by µ − e conversion are almost parallel.Thus, from Fig. 2 to Fig. 4 the constraints weaken by ∼ 4 orders of magnitude in the global SMEFT analysis.The situation worsens when nuclear uncertainties are accounted for and SI µ − e conversion can effectively constrain only 2 directions.Therefore, the bounds on the vector and scalar operators become only order 1 or even larger for the former.
Finally, with the additional simplifying assumption of no operators involving s quarks, bounds on all WCs may be derived through the global analysis.Without nuclear uncertainties, these bounds on the µ − e sector are surprisingly similar to those on the τ − ℓ sector.While with the one-at-a-time analysis much stronger constraints in the µ − e were found, in a global fit scenario, they become diluted by the very strong correlations present between the different observables.Moreover, when nuclear uncertainties are included as nuisance parameters, the bounds on the µ − e sector become significantly weaker than those in the τ − ℓ, particularly the ones for the vector and scalar operators, as previously discussed.
All in all, we find that flat directions play no role in the bounds on fully leptonic operators and the naive one-operator-at-a-time approach leads to reliable constraints on the relevant WCs.Conversely, flat directions appear and lead to fully unconstrained parameter combinations for semileptonic cLFV 4-fermion operators.While these flat directions were not found in previous global scans of the parameter space, we believe this is due to our different scanning strategy as outlined in section 5. Remarkably, we find that when the operators are those induced by the low-energy d = 6 SMEFT and if only two independent couplings (one for up and one for down quarks) are considered for each Lorentz structure, present data allows to lift all flat directions present and obtain unambiguous bounds on all operators.However, in the case of the µ − e sector, these bounds are around 4 magnitude weaker than in the oneoperator-at-a-time approach and become even weaker when nuclear uncertainties are properly accounted for.Nevertheless, extremely strong correlations among them exist, reflecting the underlying directions that are much more stringently constraint.Thus, we also provide correlation matrices in a GitHub repository that contain all these nuances as a useful tool to incorporate cLFV constraints in particular UV-complete scenarios.

A Detailed fit results
Here we present quantitative information about our bounds and the correlations among the different parameters so as to allow the implementation of our constraints in specific scenarios.
Correlations are of particular importance for the global fit results, specially for the µ − e global bounds, which get substantially degraded with respect to the one-at-a-time scenario.This is due to the fact that some specific combinations of WCs are very poorly constrained (or not at all) together with extremely tight bounds in other very specific directions.However, if a particular UV completion does not align along these weakly constrained directions, the corresponding constraints will be much tighter than those directly inferred from the plots showed in the previous sections.
Usually, this information is easily conveyed through the covariance matrix.However, all the observables considered, being cLFV and therefore not present in the SM, depend quadratically on the WCs at the leading order.Furthermore, as there is no signal in any of the observables, the best-fit point corresponds to all WCs vanishing and the resulting χ 2 test-statistics will be a purely quartic polynomial.
Indeed, if one tries to approximate the test-statistics from its Taylor's series around the best-fit point, the first non-vanishing order will be the quartic: and thus the covariance matrix vanishes for our test-statistics, such that the information about possible correlations between the coefficients will actually be contained within the co-kurtosis tensor of 4th derivatives.Ideally, a change of basis to the second order polynomials in the WC relevant for the observables considered would allow to make the χ 2 a quadratic function of the parameters analized.Unfortunately, after doing this change of basis for the most relevant observables, it is not possible to express the remaining WC combinations necessary for the rest of the measurements as a function of the new variables in an univocal way.Thus, for the sake of presenting results that can be more easily implemented when deriving constraints on specific models and that approximate the results of our global fit, we construct a "proxy-covariance matrix" whose global bounds match those extracted from the full analysis of the true teststatistics at a given confidence level.
In particular, we construct a "proxy-test-statistics" by using as observables the square root of the branching ratio, instead of the branching ratio itself.This guarantees that our "proxytest-statistics" is a quadratic polynomial in the coefficients and consequently has a well-defined covariance matrix.We then normalise the covariance matrix so as to only retain information on how correlated are the WCs.Lastly, the normalisation of the variances of each of the WCs is set to the global bound obtained in the full analysis of the true test-statistics at some given confidence level.Thus, the information about the strongly correlated directions is preserved and the bounds agree with those of the actual global fit at the confidence level selected.Notice, however, that using this covariance matrix for different confidence levels would provide wrong results as the dependence on the WCs of the true test statistics (quartic) and the approximated one (quadratic) is different.
In particular, we will give our results in terms of the global bounds {σ i } on the coefficients {c i }, with a normalized "proxy-covariance matrix" ρ, such that the "proxy-test-statistics" can be easily constructed as follows: The bounds {σ i } and their corresponding inverse-covariances ρ −1 are available in text file format in the following GitHub repository.

A.1.1 Results with u, d and s quarks
As argued in section 6, when considering all light quarks, there is a single flat direction remaining, which we isolate by defining the following uncorrelated coefficient combinations: where the ∓ sign corresponds to X = L/R, respectively.At 95% CL, the global bounds read, which are mostly uncorrelated, as shown in Fig. 6.   6: Correlations between the different τ − operators, as extracted from our numerical analysis of the LEFT from d = 6 SMEFT at low energy, considering all light quarks.Red and blue lines correspond to 68% and 95% C.L., respectively.Notice the generally mild correlations between the coefficients, which explain why the global bounds (except in the presence of flat directions) are so close to the one-at-a-time bounds (see Figs. 4 and 5).A very similar plot is obtained for the τ − µ sector.

A.1.2 Results with first generation quarks
At 95% CL, the global bounds read, A.2 µ − e SMEFT results As argued in sections 6 and 7, all WCs under consideration, except for c sV µeX , can be simultaneously constrained in our global fit.We collect these global bounds here for the analyses including the nuclear uncertainties, as well as their correlations showed in Fig. 7 and Fig. 8.
We note that, contrary to the τ sector, the correlation between some of the operators is extremely strong.This fact explains why the global bounds are so much weaker compared with the the one-at-a-time scenario (see Fig. 4), pushing some of the profiled bounds to O(≫ 1).Nevertheless, these strong correlations imply that these very weak bounds can only be saturated in models that predict very specific cancellations.Saturating all of them would in general result in incorrect results, underestimating the cLFV bounds, and thus we instead provide the correlation matrix ρ −1 available in the GitHub repository in order to correctly include these bounds on specific setups.
We also remark the big difference between the 68% C.L and 95% C.L. bounds, being the former orders of magnitude stronger for some operators.The reason, as explained before, is that at higher C.L. two constraints from SI µ−e conversion are lost due to nuclear uncertainties in the overlap integrals defining the constrained directions (see Eq. ( 31)).Therefore, the 95% C.L. global bounds presented in this work cannot be naively translated into other confidence levels and the complete χ 2 must be used.
In the following we present bounds and correlations extracted for the operators in which the electron has L-chirality, but very similar results are found in the R-chiral case, which are available in the repository.

A.2.1 Results with u, d and s quarks and nuclear uncertainties
At 95% CL, the global bounds read, with the strong correlations shown in Fig. 7.

A.2.2 Results with first generation quarks and nuclear uncertainties
Global bounds at 95% CL: with the strong correlations shown in Fig. 8.

B Bounds on dipole and four-lepton operators
For the sake of completeness, we also show the bounds on the LEFT operators as extracted from radiative and three-body leptonic decays.As discussed in section 3, the incoherent contributions of each of the operators renders the following 95% C.L. bounds totally uncorrelated:  1.9 The same exact bounds apply to the WCs obtained from interchanging (L ←→ R).

v 2 Λ 2 Figure 2 :
Figure 2: Current 95% CL bounds on LEFT cLFV operators with quarks considering only one operator at a time (see Fig. 1 for the rest of operators).Missing bars indicate that there is no (relevant) bound for those operators at present.

v 2 Λ 2 Figure 3 :
Figure 3: Current 95% CL global bounds on LEFT cLFV operators involving the three lightest quarks.All operators are considered at the same time and their WCs are profiled over to obtain individual bounds.Missing bars indicate that there is no (relevant) global bound for those WCs.For easier comparison, we depict as empty bars the one-at-a-time constraints of Fig. 2.

v 2 Λ 2 Figure 4 :
Figure 4: Current 95% global bounds on the cLFV LEFT operators with the three lightest quarks induced by d = 6 SMEFT at low energies.Color code as in Fig. 3, but now darker red bars show the effects of including nuclear uncertainties in the µ − e analysis.Bounds on pseudoscalar operators are equal to the scalar ones due to the correlations in Eqs.(47)-(49).These global bounds, along with their correlations, are collected in appendix A.

v 2 Λ 2 Figure 5 :
Figure 5: Current 95% CL global bounds on the cLFV LEFT operators with only first generation quarks induced by d = 6 SMEFT at low energies .Color code as in Fig. 4. Bounds on pseudoscalar operators are equal to the scalar ones due to the correlations in Eqs.(47)-(49).These global bounds, along with their correlations, are collected in appendix A.

Figure
Figure 6: Correlations between the different τ − operators, as extracted from our numerical analysis of the LEFT from d = 6 SMEFT at low energy, considering all light quarks.Red and blue lines correspond to 68% and 95% C.L., respectively.Notice the generally mild correlations between the coefficients, which explain why the global bounds (except in the presence of flat directions) are so close to the one-at-a-time bounds (see Figs.4 and 5).A very similar plot is obtained for the τ − µ sector.

Figure 7 :
Figure 7: Same as Fig. 6 but for the µ − e sector, showing the very strong correlations between some of the coefficients.This explains the huge relaxation of the global bounds with respect to the one-at-a-time scenario in Fig. 4.

Table 2 :
The list of relevant d = 6

Table 5 :
Summary of combinations of LEFT WCs constrained by semileptonic τ decays, given in the isospin basis u ± d. Filled boxes indicate to which WCs each observable is sensitive to, with same color and number indicating coherent contributions between those WCs.Thus, different colors/numbers correspond to the 11 independent constraints in the τ − ℓ sector.