Z lepton flavour violation as a probe for new physics at future $e^+e^-$ colliders

In this work we assess the potential of discovering new physics by searching for lepton-flavour-violating (LFV) decays of the $Z$ boson, $Z\to \ell_i \ell_j$, at the proposed circular $e^+e^-$ colliders CEPC and FCC-ee. Both projects plan to run at the $Z$-pole as a Tera Z factory, i.e., collecting $\mathcal O\left(10^{12}\right)$ $Z$ decays. In order to discuss the discovery potential in a model-independent way, we revisit the LFV $Z$ decays in the context of the Standard Model effective field theory and study the indirect constraints from LFV $\mu$ and $\tau$ decays on the operators that can induce $Z\to \ell_i \ell_j$. We find that, while the $Z\to \mu e$ rates are beyond the expected sensitivities, a Tera Z factory is promising for $Z\to \tau\ell$ decays, probing New Physics at the same level of future low-energy LFV observables.


Introduction
The overwhelming evidence for neutrino oscillations has shown that lepton family numbers are not conserved and the Standard Model (SM) needs to be extended to account for neutrino masses. Hence, there is no fundamental reason why lepton flavour violation should not occur in processes involving charged leptons only, such as i → j γ (i > j). However, the rates of charged lepton-flavour-violating (LFV) processes induced by loops involving neutrinos and W bosons are suppressed way below the sensitivity of any conceivable experiment by factors proportional to ∆m 2 /M 2 W 2 ≈ 10 −49 , where ∆m 2 denotes the squared mass differences of the neutrino mass eigenstates. As a consequence, LFV transitions are among the cleanest and most striking signals for physics beyond the SM (BSM), in fact beyond its minimal extensions accounting for neutrino masses. Fortunately, an intense experimental activity is ongoing, with a particular focus on low-energy LFV decays of leptons and mesons. For a review on LFV, including future experimental prospects, see [1]. In this article, we will instead focus on a class of processes that is only accessible at high-energy colliders: lepton flavour violating decays of the Z boson (LFVZD), Z → i j , i = j.
In Table 1, we report current limits and future prospects on the three LFVZD modes. As we can see, the old limits obtained by the LEP experiments have been recently superseded Mode LEP bound (95% CL) LHC bound (95% CL) CEPC/FCC-ee exp.
by searches performed at the LHC. However, the sensitivity of these latter searches is limited by background events following from Z → τ τ decays to such an extent that at most an improvement by one order of magnitude can be expected after the completion of the future high-luminosity run of the LHC (HL-LHC) (since HL-LHC will collect up to 3000/fb [8] of integrated luminosity, and the LHC limits in Table 1 were obtained with 20/fb − 140/fb). Leptonic colliders instead provide a much more suitable environment to tame such background: for example, the LEP limit on Z → µe obtained in [2] was based on a sample of only 4 × 10 6 Z decays in a background-free situation, i.e., no candidate events were found. 1 These processes are therefore an ideal target for future leptonic colliders. In particular, both proposed projects of circular e + e − colliders, CEPC [9,10] and FCC-ee [11,12], plan to run for several years at a center-of-mass energy around the Z pole, thus acting as a "Tera Z factory", i.e., collecting O 10 12 Z decays, about six orders of magnitude more than LEP experiments. The last column of Table 1 shows the expected sensitivity of future e + e − colliders as estimated in [7] assuming 3×10 12 Z decays (corresponding to 150 ab −1 ). As we can see, at least for the Z → τ modes, CEPC/FCC-ee could improve on the present LHC (future HL-LHC) bounds up to 4 (3) orders of magnitude, due in particular to an expected excellent momentum resolution (0.1% at 45 GeV) of the planned detectors. 2 Given the outstanding expected sensitivity of future e + e − colliders on the Z → i j decays, an obvious question is whether this will be sufficient to test or discover new physics (NP) scenarios. Indeed, in presence of new physics leading to Z → i j , low-energy LFV processes are unavoidably induced by the virtual exchange of the Z itself: i → j Z * → j ff , where f is a SM quark or lepton. These processes are subject to strong constraints, which then translate into indirect bounds on LFVZD rates [13][14][15][16]. In this work, we plan to reassess the maximal possible LFV effects in Z decays in a model independent way, in view of the improved present 1 A major advantage of a leptonic collider is the knowledge of the momenta of the colliding partons, so that the constraint on the invariant mass of the two leptons m 2 i j = m 2 Z can be precisely implemented up to the beams energy spread, in contrast to the LHC where it is limited by the (large) width of the Z. 2 The sensitivity on BR(Z → µe) is limited to ∼ 10 −8 by backgrounds from Z → µµ with one of the muons releasing enough bremsstrahlung energy in the ECAL to be misidentified as an electron [7]. Only in presence of improved electron/muon separation methods a sensitivity down to 10 −10 could be achieved. and future searches for muon and tau LFV decays. From the theoretical point of view, the absence of signals in direct searches for the production of new particles at the LHC suggests an energy gap between the electroweak scale and the scale where new physics inducing LFVZD may exist, prompting us to work within the context of an effective field theory (EFT), i.e., introducing a set of higher-dimensional gauge-invariant local operators to be added to the usual SM Lagrangian. These operators, built out of SM fields and suppressed by inverse powers of the new physics scale, can parameterise the effects of any kind of NP models as far as the experimentally accessible energies are lower that the actual NP energy scale. Such an effective theory is known as the Standard Model Effective Field Theory (SMEFT) [17,18] (for a recent review see [19]) and provides the optimal framework for a model-independent analysis. In this context, LFV processes could be induced by dimension 6 effective operators, as discussed in detail in [20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38]. In this article we will employ the SMEFT framework to study low-energy constraints on LFVZD. This type of decays have been also studied within several UV-complete models, such as heavy sterile neutrinos [39][40][41][42][43][44][45], supersymmetry [22,46], leptoquarks [47], or in scenarios with extended gauge sectors [48][49][50][51]. They have also been previously explored in the context of SMEFT in [16,24,25,28,33].
The outline of this article is the following. In Section 2 we describe the effective field theory setup that has been employed for our analysis. How the SMEFT operators can induce LFVZD is shown in Section 3. In Section 4, we discuss how indirect constraints on Z → i j arise from low-energy LFV observables. The next section contains the results of our analysis, where we assume that the UV physics induces a single dominant operator (Section 5.1) or multiple operators that could possibly interfere (Section 5.2). We summarise and conclude in Section 6 while a number of useful analytical formulae and results are shown in the Appendix.

Lepton flavour violation in the SMEFT
Throughout this work, we will assume that the new particles related to the NP scale Λ responsible for LFV effects are quite heavy (Λ m W ) and that there are no other particles in between these scales.
In such a scenario, it is convenient to work in the SMEFT framework, where the basic idea is to parameterise the low-energy effects of the high-energy theory in terms of higher dimensional operators and the associated Wilson coefficients. More specifically, the Lagrangian will consist of that of the SM extended with a tower of higher-dimensional operators suppressed by inverse powers of Λ: where L SM contains renormalizable operators up to dimension-4, Q a are the effective operators of dimension-n and the C (n) a represent the corresponding Wilson coefficients (WCs) which

4-lepton operators
Dipole operators Lepton-Higgs operators , while B µν and W I µν are the U (1) Y and SU (2) L field strengths, and τ I with I = 1, 2, 3 are the Pauli matrices. Flavour indices are not shown. depend on the renormalization scale µ. In a given UV-complete model, the Wilson coefficients at the scale Λ can be determined by integrating out the heavy particles. In our modelindependent approach, however, we will consider the C (n) a (Λ) as independent free parameters. The dominant contributions to our LFV processes are then expected to come from dimension-6 operators, 3 hence we do not consider higher dimension operators in our analysis. Out of the 59 (without counting the combinations of flavour indices) dimension-6 operators [17,18] that can be constructed from SM fields and respect SU (3) C × SU (2) L × U (1) Y invariance (and baryon-number conservation), only a subset is relevant to us, namely the operators which contribute to LFV processes at tree-level or at 1-loop level. These are (i) 4-lepton operators, (ii) leptonic dipole operators, (iii) lepton-Higgs operators, (iv) 2-lepton 2-quark operators. A complete list of all 6-dimensional operators that can induce LFV processes as discussed in [25] is displayed in Table 2. In particular, among these operators, we highlight in blue those that modify the Z boson couplings to leptons and can therefore contribute to Z → i j decays at the tree level, that is, • the dipole operators C ij eW Q ij eW and C ij eB Q ij eB ; • the lepton-current Higgs-current operators C where we have explicitly shown the notation for the WCs and flavour indices that we are going to adopt throughout the work 4 . In order to impose experimental constraints on the coefficients of higher dimensional operators, one needs to evaluate the Renormalisation Group (RG) running from the scale Λ to the energy scale relevant for a given experiment. As customary, this is done by solving the RG evolution (RGE) for the higher-dimensional operators, which could not only change the value of a given WC, but also induce mixing between different operators. This means generating at low energies operators whose coefficients were set to zero at the scale Λ. For a LFV observable, this procedure may consist of up to three steps. First, the SMEFT RGE, known at one loop [53][54][55], will run the WCs from Λ to the electroweak scale ∼ m Z,W . This running already allows us to extract information relevant for the LFV decays of the Z or Higgs bosons, however it is not enough for energies below the electroweak scale. Therefore, in the second step the SMEFT is matched to the so-called Low-Energy Effective Field Theory (LEFT) [56,57], after integrating out the top quark and the W ± , Z and Higgs bosons. The final step would then be the QED×QCD running of the LEFT operators to the low-energy scale, m µ or m τ , according to the experimental observable that has to be evaluated. We notice that the QED running plays an important role for low-energy LFV observables as it can induce large operator mixing [32]. We implement all these steps for each of the LFV observables entering our numerical analysis with the help of the wilson [58] and flavio [59] packages.
Note that the above lepton-Higgs operators that could be responsible for Z → i j do not induce couplings of the physical Higgs particle with leptons, thus they cannot give rise to LFV Higgs decays, h → i j . The latter processes can be induced at the tree level only by the last lepton-Higgs operator of Table 2, Q eϕ3 , see e.g. [60,61]. Such an operator can be generated by the RG running if at least one of the five operators mentioned above is induced by the UV theory. Nevertheless this effect is proportional to at least one power of a small lepton Yukawa coupling (due to the necessary flip of the chirality of the leptons) and thus their contributions are substantially suppressed [54,55]. 5 Similarly, these five operators are not generated by running effects controlled by C eϕ3 (not at one loop). Hence, in the context of the SMEFT, Z and Higgs LFV effects are practically decoupled and we are not going to discuss the latter in this work. 4 In the case of 4-fermion operators, there might be redundant flavour combinations, e.g., C ijkl , C klij , C ilkj , C kjil giving rise to the same operator (directly or using Fierz identities). We will work in a non-redundant flavour basis and consider only one of these WCs. 5 For example, we find that C µτ eϕ3 (m h ) ≈

Lepton flavour violating Z decays
The effective interactions involving the Z boson and the SM leptons, including those responsible for LFV effects, are given by the following Lagrangian [22] where are the SM couplings of Z to right-handed (RH) and left-handed (LH) lepton currents respectively, with s w (c w ) being the sine (cosine) of the weak mixing angle. New physics effects are encoded in the effective couplings δg V /T , which at the tree level match the SMEFT operators as follows where the WCs have to be evaluated at the scale µ = m Z . The branching ratios of the Z decays into leptons, in particular of the LFV modes, are then given by the following expression [22,25] where Γ Z = 2.4952(23) GeV is the total decay width of the Z boson [62], and we summed over the two possible combinations of lepton charges, ± i ∓ j . As anticipated above, only five SMEFT operators (those highlighted in Table 2) can induce the LFV Z → i j decays at the tree level. Actually, as we can see from Eqs. (4,5), only three independent combinations of the corresponding WCs contribute: As already mentioned, these WCs are to be evaluated at µ = m Z . On the other hand, the SMEFT running induces mixing of various operators and therefore the LFVZD will be sensitive to more C i (Λ) beyond those five operators explicitly entering in Eq. (6). In our numerical analysis, we have implemented the full one-loop SMEFT running by means of wilson [58]. Nevertheless, given the large number of dimension-6 operators, it is helpful to identify beforehand those WCs that will be more relevant for the LFVZD after the one-loop RG evolution.
For the sake of the following discussion, let us focus on the vectorial couplings in Eq. (4), and consider only the gauge terms in the one-loop RGE [55], thus neglecting terms controlled by the Higgs self-coupling and Yukawa couplings. In that case, the running can induce LFV terms of any of the three Higgs-lepton operators from 24 additional WCs, although not all of them contributing with the same strength. Schematically, we can write it as with and the matrices γ ij encoding the anomalous dimensions, whose explicit form can be obtained from Ref. [55]. Despite the fact that the RGE will mix these 27 operators, from the structure of Eq. (8) it is clear that the contribution from the fourteen WCs in C 3 and C 4 will be negligible for the LFVZD, as they do not generate directly any Higgs-lepton operator. The vectorial contributions to the LFV Z decays will then come mainly from three WCs in Eq. (9) and one-loop suppressed contribution from ten WCs in Eq. (10). A similar classification can also be done for the dipole operators C eW and C eB . In this case, we can identify a set of 10 operators forming a closed RGE system, including themselves and other operators involving gauge and Higgs fields, not shown in Table 2. As we will show in Section 5.1, however, dipole operators are so constrained by low-energy LFV observables that result completely irrelevant for the sake of LFVZD, and we therefore refrain from discussing further their RGE effects.

Indirect constraints from low-energy LFV observables
The operators contributing to the LFV Z decays also contribute to low-energy LFV observables, as shown in Fig. 1. Among the most relevant observables, we find the µ to e conversion in nuclei, the radiative decays i → j γ, the leptonic decays i → j k¯ m , and the semileptonic decays τ → M , with M being a pseudoscalar or vectorial meson. All these processes are severely constrained experimentally, see Table 3, and thus set indirect constraints on the WCs we are considering and on the maximum allowed rate for the LFV Z decays. 6 As we will be dealing with low-energy observables, we need to use the full QFT machinery described in Section 2, i.e., we need to run the WCs at µ = Λ down to the electroweak scale using the SMEFT running, then match them to the LEFT operators and finally use the LEFT running down to µ = m τ or m µ . A complete LEFT basis, the one-loop RGE and the tree-level matching to SMEFT operators can be found in Refs. [56,57]. In our case, the most relevant LEFT operators will consist of four fermions with at least one spinor combination consisting of two different lepton flavours, and the photon dipole operator. Schematically, the former will take the form where f α,β refers to any light quark or lepton, P X,Y to left or right projectors and Γ A to the scalar, vector or tensor operators, i.e., Γ S = 1, Γ V = γ µ or Γ T = σ µν . On the other hand, the relevant photon dipole operator is given by which directly matches, at tree-level, to the orthogonal combination of C eZ , cf. Eq. (7), i.e., The case of the dipole has been extensively studied also beyond leading order, showing the relevance of performing the matching at one-loop [27] and including the leading two-loop anomalous dimensions [32]. Nevertheless, we will see that the naive analysis including only the tree-level matching and the one-loop anomalous dimensions already suggests that the role of the dipoles will be suppressed in the LFVZD due to the strong constrains from the radiative decays. Therefore, we will not include any of those higher order contributions in our numerical analysis.
LFV obs.  Once again, all this procedure is included in our numerical analysis using wilson and flavio. Nevertheless, and with the aim of improving our understanding of the numerical results, we give in Appendix A the analytical expressions for these observables in terms of the C i (Λ) after tree-level matching, but neglecting the RGE effects. We checked that these formulas provide a good approximation for the leading contributions, e.g. the contributions of C ϕe or C u to τ → ρ , and will describe the overall behaviour of most of our numerical results. Note however that they neglect RGE-induced mixing that might be important in some cases, such as the role of C in τ → ρ , which is a consequence of the SMEFT running in Eq. (8) and the equivalent LEFT one. We will discuss the relevance of these latter contributions while presenting our numerical analysis.
Finally, it is worth mentioning that the operators shown in Table 2 can also be probed at high-energy colliders, in particular the 2-lepton 2-quark operators that we are going to consider in the next section. Nevertheless, current LHC sensitivities for the operators we are interested in are not competitive with those at low-energy experiments [78].

Results and discussion
In Section 3 we have seen that there are five SMEFT operators, in three independent combinations displayed in Eq. (7), that directly generate LFV Z decays at the tree level and therefore they will give the dominant contributions to our observables. For this reason, it is interesting to analyse their effect when only one of these operators is present, assuming other operators are absent at that time, and to compare the potential of the LFVZD to that of other low-energy LFV observables in the search for new physics. Nevertheless, it is reasonable to expect that a given UV theory will generate several SMEFT operators at the same time. In such a case, possible interference between different operators could distort the results obtained considering each operator individually. In order to illustrate this idea, we will consider two non-zero operators at a time, paying special attention to possible flat directions, i.e., fine cancellations among contributions stemming from different operators. We will see how other operators beyond those in Eq. (7) could play an important role in this case, and we will assess whether they could suppress some of the low-energy LFV decays, thus allowing larger LFV Z decay rates. 7

Single operator dominance
We want to analyse the effect of the five relevant SMEFT operators for tree-level LFVZD, i.e., the two dipoles and three Higgs-lepton operators, assuming that only one of them is present at a time, which would be approximately the case if the underlying UV dynamics mostly matches to a single operator while inducing others at a substantially suppressed level. This hypothesis, however, needs to be defined at a given scale µ.
We start considering a single non-zero WC at µ = m Z , as this is the relevant scale for the Z decays. Table 4 shows the results we found for the maximum allowed LFVZD rates after requiring that all the the low-energy LFV processes that are induced by the same operators lie below their current experimental limits shown in Table 3. Notice that, in this case, the RGE effects do not affect the LFVZD, which can be directly computed by means of Eq. (6). This is not true for the low-energy observables, for which the proper matching and RGE to µ = m τ or m µ need to be taken into account, as discussed in Section 2. From Table 4 we see that dipole operators are extremely constrained by the radiative decays i → j γ, which translates to LFVZD rates beyond current and future experimental sensitivities. At this point, one could be tempted to tune C eB and C eW in such a way that the photon dipole in Eq. (15) vanishes, avoiding the bounds from i → j γ and maximizing LFVZD via the Z dipole. Nevertheless, this choice would have several drawbacks. Firstly, it seems very unlikely to have a UV model that leads to a vanishing C eγ (m Z ). Notice that even 7 In what follows, we assume real WCs for simplicity. In fact, the total decay rates for the LFV Z decays can not be affected by the phase of the WCs, and, even in case of discovery, we expect that the number of LFVZD events would be insufficient to perform an analysis in search of CP-violating effects.

Operator
Indirect Limit on LFVZD Strongest constraint if the UV model generated only C eZ at the NP scale µ = Λ, the RGE would induce a non-zero photon dipole at µ = m Z . This means that a huge fine-tuning between C eγ and the radiative effects would be needed to have C eγ (m Z ) = 0. Secondly, a vanishing photon dipole would only suppress the tree-level contributions to i → j γ, however higher order terms would still be important [25,27,32]. Although not included in Table 4, we have estimated the size of these higher order effects following [27] and found that the radiative decays would still impose strong bounds even in the extreme case of vanishing C eγ (m Z ), setting indirect limits on dipole mediated LFVZD beyond future sensitivities. On the other hand, Higgs-lepton operators, which do not generate i → j γ at the treelevel, 8 are less constrained and larger LFVZD are allowed. In the µ-e sector, the strongest current bounds are imposed by µ-e conversion in nuclei. This translates into an indirect bound of BR(Z → µe) 10 −13 , which unfortunately is still beyond the reach of future experiments, see Table 1. The results for the tau sector are however more optimistic. In this case, they are currently mostly constrained by τ → ρ decays, which imposes indirect limits of the order BR(Z → τ ) 10 −8 − 10 −7 . While still below the reach of current LEP/LHC bounds (as well as the expected HL-LHC sensitivity), these decay rates could be probed at a future Tera Z factory. Next, we consider as input the WCs at the NP scale µ = Λ, as it would be the scale at which we could integrate out the new heavy fields related to the UV dynamics and generate the SMEFT operators. In order to compare better the sensitivity reach of future high-and low-energy LFV experiments to such operators, we focus our discussion on what would be the largest NP scale that we could probe in each case. Under our hypothesis of switching on a single operator at a time, the LFV observables under consideration would scale as C(Λ) 2 /Λ 4 , up to some O(log µ/Λ) corrections from the RGE. This means that, for a given experimental upper limit on an observable, the maximum NP scale that we could be probing corresponds to the maximum value for the WC at that scale. Notice that higher NP scales would require nonperturbative WCs, while smaller scales would always be allowed provided that the WC is small enough. We show in Figure 2 the sensitivities for Λ for our most relevant operators and from different observables, where we have assumed that C(Λ) ≤ 1 from perturbativity arguments. In this case, and opposite to Table 4, we choose the Z dipole operator as input, which implicitly assumes C eγ (Λ) = 0, since this hypothesis is still challenging but more plausible at µ = Λ. We also show Q From Figure 2 we can see that current sensitivities (solid bars) are always worse in the case of the LFVZD than from low-energy observables, in agreement with our findings in Table 4, and especially in the case of the dipoles. Despite we chose to switch on only the Z dipole and not the photon one at µ = Λ, the RGE generate a photon dipole at low energies, providing a better sensitivity to NP from low-energy observables even in this extreme case. Unfortunately, the situation will be similar for future experiments, therefore we can conclude that LFVZD will not be competitive probing the dipole operators and ignore them for the rest of the analysis.
On the other hand, the Higgs-lepton operators are again more promising, in particular in the τ -sectors. We can clearly see how the huge improvement in sensitivities at the Tera Z factories will boost the potential of the Z → τ decays, reaching values that are competitive with low-energy observables. This result is remarkable and very promising, since the future limits for LFV τ decays shown in Table 3 are based on the most optimistic assumption that Belle-II searches will be background-free, which does not necessarily need to be the case. Therefore, being able to access the same NP scale with an independent high-energy observable would be of great value. This will not be the case, however, in the µ-e sector, where the future reach on Z → µe lies below the current bounds from low-energy processes, even considering the most optimistic sensitivity BR(Z → µe) ∼ 10 −10 as we did in Figure 2. Indeed, the outstanding future reach of µ → eee and µ → e conversion experiments will probe scales almost two orders of magnitude above those accessible at Tera Z factories through Z → µe.
In summary, we have seen that if a single operator dominates, future Tera Z factories searching for LFVZD could probe NP at the same level as low-energy experiments, in particular for the τ -sector of the Higgs-lepton operators. Notice that so far we have only considered the leading order operators for LFVZD, however we have seen that other operators can also be relevant due to RGE effects, in particular those of Eq. (10). Under this single operator dominance hypothesis, we can expect that these new operators lead to small LFVZD rates, as they are constrained by processes that they do generate at tree level. Nevertheless, they could still induce new interesting effects when combining several operators at a time, as we are going to explore in the following.

Interference of multiple operators
We now move to consider possible interference effects arising when multiple operators are induced by the UV dynamics at the scale Λ. As argued above, in the case of dipole operators constraints from LFV radiative decays are so difficult to overcome that we cannot expect these operators to be a substantial source of LFVZD. In the following discussion, we therefore focus on the three remaining Higgs-lepton operators.
As we have seen, the combination C (m Z ) contributes to Z → i j . Furthermore, we notice that, after decoupling the Z boson, the very same combination, at the same scale, also matches to the LEFT 4-fermion operators relevant for low-energy LFV processes, as can be seen by inspecting the formulae in the Appendix A. As a consequence, even in presence of a cancellation between C ϕ are practically identical for both LFVZD and low-energy processes. Hence, in the following, we will only show results for C (1) ϕ . Notice also that the RH operator C ϕe cannot interfere with LH Higgs-lepton operators neither for LFVZD, see Eq. (6), nor for low-energy LFV processes, cf. Appendix A. Therefore, if the UV physics mostly induces Higgs-lepton operators no cancellations are possible and the indirect limits reported in Table 4 are all non-vanishing. In order to study the possibility of non-trivial cancellations, we then have to consider simultaneous presence at the scale Λ of a non-zero coefficient of a Higgs-lepton operator (we start with C (1) ϕ ) and one or more 4-fermion SMEFT operators, such as 4-lepton or 2-lepton 2-quark operators (cf. Table 2) that are relevant for the low-energy constraints.
µ − e sector. In Figure 3, we plot contours of BR(Z → µe) on the plane C . We choose Λ = 1 TeV for this and the following Figures, nevertheless our results are qualitatively the same for other scales, following the scaling C i (Λ)/Λ 2 with small logarithmic modifications from the RGE. The lighter coloured regions are currently allowed by the µ → eee (blue) and µ → e conversion (orange) constraints, while the corresponding darker regions show how the allowed parameter space will shrink given the future expected sensitivities reported in Table 3.
From the examples in these plots, we can see that a value of the 4-fermion WC at the scale Λ = 1 TeV approximately three orders of magnitude larger than the value of C (1) eµ ϕ (Λ) and its RGE running from Λ to m Z , cf. Section 3. However, we can also see that such cancellations would occur for values of the parameters that are already excluded by the combination of the bounds from µ → eee and µ → e conversion.  The plots of Figure 3 also show how difficult it is to tune the parameters in such a way that LFVZD effects are enhanced relative to the muon LFV processes. From the top panel we can see that there is a flat direction cancelling the conversion rate of µ → e in nuclei for , as a consequence of the coherent nature of the conversion process and the consequent interference of the various contributions shown in Eq. (62). This numerical result is well reproduced by the analytical approximate formulae collected in the Appendix A.5. From these expressions it is easy to check that the cancellation requires if only these vector operators are involved as in our example. This also shows that a tuning of the parameters to obtain CR(µ → e, N ) ≈ 0 for a certain nucleus would not lead to an exact cancellation of the conversion rate in another nucleus, for which the overlap integrals V (n) and V (p) have different numerical values [79]. 9 Moreover, µ → eee is clearly unaffected by such a tuning as shown by the figure. In the displayed example, this implies an indirect bound BR(Z → µe) < 10 −12 even along the flat direction, which, although being an order of magnitude better with respect to the single operator analysis in Table 4, still lies beyond the expected reach at the Tera Z factory. The bottom panel of Figure 3 shows what happens if the dominant operator generated at high energies alongside C (1) eµ ϕ is a 4-lepton operator. An exact cancellation of µ → e conversion in nuclei is still possible, albeit it is due to radiative effects in this case. Indeed, the RGE running below the electroweak scale induces 2-lepton 2-quark operators through loops involving the 4-lepton operator, with coefficients of the order (α/4π) C eµee e log(m Z /m µ ). This explains why the flat direction of the orange region corresponds to C eµee e (Λ) ≈ 10 2 ×C (1) eµ ϕ (Λ). We can also see that µ → eee is only mildly suppressed (along another direction), which implies a combined indirect constraint of about BR(Z → µe) 10 −14 . More in general, no exact cancellation can occur in presence of only two dominant WCs contributing to µ → eee. The reason is that, as one can see from the formulae in Appendix A.2, a strong suppression of µ → eee would require a simultaneous cancellation of the coefficients in Eqs. (23,24), which is only possible if the coefficient of the RH Higgs-lepton operator C (1) eµ ϕe is also tuned to cancel out C eµee e in Eq. (24). Needless to say, yet another tuning would be required to suppress µ → e conversion, such as in the top panel of the figure. 10 In summary, a simultaneous cancellation of µ → eee and µ → e conversion in nuclei able to enhance the maximum possible BR(Z → µe) would require a fine tuning involving at least four operators and thus it looks extremely unlikely in the context of any UV-complete model. Therefore, we can conclude that we do not expect Z → µe to be observed at the CEPC or FCC-ee even if the optimal sensitivity shown in Table 1 (10 −10 ) could be reached. Incidentally, note the astonishing sensitivity of the upcoming low-energy LFV experiments Mu3e, Mu2e and COMET, which will constrain the coefficients of certain operators down to even 10 −7 -10 −8 if LFV effects are induced by TeV-scale new physics.
τ − sector. In Figures 4 and 5, we show our results for Z → τ e and the LFV τ decays τ → eee, τ → πe, and τ → ρe, bearing in mind that the corresponding plots for the τ − µ sector are virtually the same once the obvious changes (e.g. of flavour indices) are taken into 9 Notice that indeed the displayed current and future bounds follow from muon conversion in different nuclei, respectively Au and Al, cf. Table 3. However, the quantity appearing in the right-hand side of Eq. (16) is numerically similar in the two cases (0.19 for Au, 0.17 for Al), such that the slight difference in the flat direction is difficult to appreciate in a logarithmic plot as those in Figure 3.
10 Moreover, the above-discussed operators involved in µ → eee would destabilise the cancellation of µ → e conversion and a fine adjustment of the flat direction shown in the top panel would be needed. As in the case of the µ − e sector, we observe that the WCs could conspire to cancel out BR(Z → τ e), a possibility that is however incompatible with the low-energy constraints.
Conversely, Figure 4 shows that τ → πe and τ → ρe can be suppressed, but in general not simultaneously. This is due to the different dependence of the amplitudes of the two semileptonic decays on the WCs. For instance, from the top panel, we see that τ → πe vanishes for C Both numerical values are well accounted for by the approximate formulae in Appendices A.3 and A.4, see Eqs. (44,54). This result implies that it is not possible to tune the parameters to cancel both semileptonic τ decays and enhance the LFVZD effects. This is also shown clearly in the bottom panel of the figure, where τ → ρe features a flat direction while τ → πe does not. Notice that the latter is actually almost independent of C eeeτ , since the RGE contributions of the vectorial 4-lepton operators are very suppressed for this channel involving a pseudoscalar meson.
For both cases depicted in Figure 4, we can see that suppressing the purely leptonic LFV decays such as τ → eee would require introducing and tuning additional operators, hence it is difficult to envisage a situation where BR(Z → τ e) can attain values much larger than the limits we obtained in Table 4 for the single operator analysis. Still, note that, once the combined constraints are fulfilled, BR(Z → τ e) can be as large as ≈ 10 −7 , a rate which is largely within the sensitivity of a Tera Z factory, ≈ 10 −9 , see Table 1. Consequently, the plots show the nice complementarity in testing the NP parameter space of LFVZD searches at a Tera Z and the prospected sensitivity on LFV τ decays at Belle-II.
So far we have considered the LH Higgs-lepton operator C (1) eτ ϕ (the results for C (3) eτ ϕ being almost identical). Figure 5 shows that the outcome of our analysis does not qualitatively change if we consider operators involving RH currents: in this example, C eτ ϕe and the 2-lepton 2-quark operator C eτ uu eu . The different behaviours of τ → πe and τ → ρe are even more striking due to a relative sign difference of the C eτ ϕe contributions, see Eqs. (45,55), accentuating the complementarity between different LFV observables. As before, we find maximum allowed LFVZD rates well within the expected sensitivities at the Tera Z factory.
To summarise, also for the observables involving τ leptons, we found that accidental can-cellations of either Z → τ or τ LFV decays are extremely unlikely. Moreover, Figures 4  and 5 show that, if the UV physics induces several operators simultaneously with coefficients of comparable size (let's say, 10 −3 /(1 TeV) 2 ) including one or more lepton-Higgs operators, the capability of testing new physics of a Tera Z factory through searches for Z → τ is comparable to or better than the one of Belle-II through τ LFV processes.

Summary and Conclusions
In this work we assessed the potential of the proposed circular e + e − colliders CEPC and FCC-ee to discover new physics by searching for the LFV decays Z → i j with data collected operating at a center-of-mass energy ≈ m Z , and under the assumption that an integrated luminosity corresponding to O(10 12 ) visible Z decays ("Tera Z") will be reached. As shown in Table 1, a Tera Z factory has the potential to improve the present limits by more than three orders of magnitude. It is however not obvious whether such an excellent sensitivity is enough to test or discover new physics effects, since any new dynamics responsible for LFV interactions of the Z bosons would unavoidably induce LFV decays of muons and taus, which in turn set indirect limits on the maximum possible rates of Z → i j . We evaluated these indirect constraints within the model-independent framework of the SMEFT. The main findings of our analysis can be summarised as follows.
• Of the five SMEFT operators that can induce LFVZD at tree level shown in Table 2, only three, belonging to the class of the lepton-Higgs operators can lead to observable effects. The other two, that is, the dipole operators, are too tightly constrained by the low-energy LFV radiative decays i → j γ, see Table 4 and Figure 2.
• Combining the current bounds on BR(µ → eee) and on the µ → e conversion rate in atomic nuclei results in a strong indirect constraint on BR(Z → µe). The distinct dependence of the two processes (and also of the conversion rates in different nuclei) on LFV 4-fermion operators makes accidental cancellations that could enhance the allowed Z → µe rate very unlikely. This would indeed require a careful tuning of the coefficients of at least four independent operators. As a consequence, in a realistic UV model, the current (future) low-energy bounds imply BR(Z → µe) 10 −12 (10 −16 ) (see Figure 3), way beyond the expected sensitivity of Tera Z factories.
• Also in the case of τ − transitions, the interplay of distinct LFV τ decays (in particular τ → ρ and τ → π that, despite involving the same fermions, have a prominently different dependence on LFV operators) exclude the possibility that low-energy constraints can be realistically tuned away, cf. Figs. 4 and 5.
• Despite the constraints from low-energy LFV processes, BR(Z → τ ) is still allowed to be as large as ≈ 10 −7 , more than one order of magnitude below the current LHC bounds (and thus beyond the reach of HL-LHC too), but well within the sensitivity of future Tera Z factories, ≈ 10 −9 , cf. Table 1. This is clearly a consequence of the fact that τ LFV processes are comparatively less constrained than the muon ones.
• Even considering the most optimistic future bounds from Belle-II, we showed in Figs. 4 and 5 that Tera-Z searches for LFV Z decays will be able to test new physics as effectively as LFV τ decays (and up to scales of the order of 20-30 TeV as shown in Figure 2). This is the case even if the UV new physics gives rise not only to the lepton-Higgs operators responsible for Z → τ , but also directly to 4-fermion operators, as long as the coefficients of the latter are of the same order as those of the former (or smaller). On the other hand, in order for the LFV Z decays to be observable, the coefficients of the dipole operators need to be somewhat suppressed (which is certainly possible, since they are unavoidably generated at loop level).
In summary, we have shown that a future Tera Z factory has the potential to probe new LFV physics in the τ − sector at the same level of other low-energy experiments. They can thus not only provide new independent insight into these phenomena, but also play a major role if Belle-II does not reach the most optimistic background-free environment. Finally, it is worth stressing the importance of the Z → µe searches even if they are disfavoured by our analysis, as their experimental observation would imply a new physics scenario even beyond the SMEFT framework.

A.1 Radiative decays i → j γ
The branching ratio of radiative lepton decay processes can be expressed in terms of two form factors F T L and F T R [25]: At tree-level, these form factors match to the SMEFT photon dipole defined in Eq. (15) as Additional contributions beyond tree-level have also been computed, which could be important when the photon dipole is not generated at tree level. See for instance Refs. [25,27,32] for further details.

A.2 Leptonic 3 body decays i → j k¯ m
Given some small discrepancies between different available computations [21,22,25,32], we compute again the decay rates for i → j k¯ m . Starting from the LEFT Lagrangian [56], the relevant terms for the tree-level 3 body decays are ē j γ µ P L e i ē k γ µ P R e m + C SRR ee ē j P R e i ē k P R e m + h.c.
The expressions for the decays depend on the flavor combinations of the final leptons, as they could involve new possible contractions and symmetry factors. For the µ → eee, τ → eee and τ → µµµ decays we get, with the matching conditions given by and the photon dipole matching given in Eq. (15). Similarly, for the decays τ → eµμ and τ → µeē, we find, with Finally, for τ → eeμ and τ → µµē, with Our results are in agreement, when comparison is possible, to those of [21,22,32].

A.3 Semileptonic τ decays τ → V
The expressions for τ → V , with = e, µ and V = ρ, φ a vectorial meson, in terms of the LEFT Wilson coefficients can be found in Ref. [34]. In order to express them in terms of the SMEFT parameters, we use the tree-level matching relations from Ref. [56]. Here we report the resulting expressions. The branching ratio can be expressed as with λ(a, b, c) = a 2 + b 2 + c 2 − 2(ab + ac + bc), and the amplitude M receiving contributions of both vectorial and tensorial leptonic currents: From here, the squared amplitude is given as with and the interference After matching, the vectorial current receives contributions from the 2 lepton-2 quark vectorial operators and from the Higgs-lepton ones, while the tensorial current gets contributions from the dipoles and the tensorial Q equ in Table 2. The final expressions depend on the meson and for V = ρ, φ they are given by ϕl + C with f (T )V the (transverse) decay constants of the vectorial mesons.
A.4 Semileptonic τ decays τ → P Equivalently to the vectorial mesons, we take the expressions for the decays to a pseudoscalar meson P from Ref. [34] and apply the tree-level matching relations from Ref. [56]. The branching ratio for τ → P , with = e, µ and P = π 0 , K 0 , is given again by Eq. (38) after replacing V by P. The squared decay amplitude in this case can be expressed as, with g τ P L = g τ P SL − m g τ P V L + m τ g τ P V R , g τ P R = g τ P SR − m g τ P V R + m τ g τ P V L .
The expressions for the scalar and vectorial couplings depend on the meson and for P = π 0 , K 0 they are given by with f P the decay constants of the pseudoscalar mesons.

A.5 µ → e conversion in nuclei
The conversion rate can be expressed as [79], with Γ capt being the muon capture rate in a nuclei N and D, S (p/n) and V (p/n) the dimensionless overlap integrals, whose numerical values depend on the nuclei and can be found in Ref. [79]. After tree-level matching [56], we obtain that the dipole form factors are given by the vector form factors by with C (u) and finally the scalar form factors by Notice that a more appropriate description of the the scalar form factors, using the effective Lagrangian at the nucleon level and including also the effective interaction with gluons, should be considered. Nevertheless, the scalar operators do not enter our analysis and therefore we refer to Ref. [80] for further details.