$\tau \to \mu\mu\mu$ at a rate of one out of $10^{14}$ tau decays?

We present in a full analytic form the partial widths for the lepton flavour violating decays $\mu^\pm \to e^\pm e^+ e^-$ and $\tau^\pm \to \ell^\pm \ell'^{+} \ell'^{-}$, with $\ell,\ell'=\mu,e$, mediated by neutrino oscillations in the one-loop diagrams. Compared to the first result by Petcov in [1], obtained in the zero momentum limit $\mathcal{P}\ll m_{\nu} \ll M_W$, we retain full dependence on $\mathcal{P}$, the momenta and masses of external particles, and we determine the branching ratios in the physical limit $m_\nu \ll \mathcal{P} \ll M_W$. We show that the claim presented in [2] that the $\tau \to \ell \ell' \ell'$ branching ratios could be as large as $10^{-14}$, as a consequence of keeping the $\mathcal{P}$ dependence, is flawed. We find rates of order $10^{-55}$, even smaller than those obtained in the zero momentum limit, as the latter prediction contains an unphysical logarithmic enhancement.


INTRODUCTION
It is reported by several experimental collaborations, e.g. by CMS [3], ATLAS [4], LHCb [5], BABAR [6][7][8][9] and Belle [10], that the branching ratios for the charged lepton flavour violating (CLFV) decays τ ± → ℓ ± ℓ ′+ ℓ ′− , with ℓ, ℓ ′ = e, µ, can be as large as 10 −14 in the Standard Model extended with either a Dirac or a Majorana mass term for neutrinos.This follows from a claim by Pham in [2] that for these decays the GIM mechanism [11] produces a suppression of the form | i U ℓi U * Li log x i | 2 , where x i = m 2 νi /M 2 W , i = 1, 2, 3, m νi and M W are the masses of the three neutrinos and the W boson, and U is the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) mixing matrix [12,13].The result in [2] is in sharp contrast with the first evaluation by Petcov [1], which showed that these CLFV decays are instead power suppressed by | i U ℓi U * Li x i log x i | 2 , so that the smallness of the ratios m νi /M W crushes the branching fractions well below 10 −54 , far beyond the sensitivity of any foreseeable experiment.
The calculations in [1] and [2] differ as follows.Ref. [1] employed for the evaluation of the one-loop diagrams the zero-momentum-limit (ZML) approximation, which assumes vanishing momenta and masses of the external particles while it retains the dependence on the internal masses of neutrinos and the W boson.The ZML implicitly assumes the mass scale hierarchy (ZML) P ≪ m νi ≪ M W , where P generically stands for any of the external particle momenta and masses, e.g.P ∼ m L or P ∼ m ℓ .This approximation, even if far from the physical situation, allows a substantial simplification of the one-loop integrals, as in this way they depend only on x i .Ref. [2], on the contrary, argued that once the external momentum dependence is taken into account, the GIM cancellation in L → ℓℓ ′ ℓ ′ , with L = τ or µ, becomes actually much milder, with a suppression only of the form , which leads to branching ratio values of the order of 10 −14 .
If the prediction in [2] were true, it would imply, with current values for neutrino mixing angles and mass splittings [14], that the branching ratio of µ → eee could reach 10 −17 , for a lightest-neutrino mass of the order of 10 −10 eV or smaller.This would be just around the corner for the Mu3e experiment currently under development at the Paul Scherrer Institute in Switzerland, which aims to reach a sensitivity of Br(µ → eee) ∼ 10 −16 [15].For the tau, the rates would be in the range 10 −16 −10 −13 , still several orders of magnitude smaller than the current world averages, Br(τ → ℓℓ ′ ℓ ′ ) 10 −8 [16], and the expected sensitivity of Belle II, Br(τ → ℓℓ ′ ℓ ′ ) 10 −10 [17], and the HL-LHC, Br(τ → ℓℓ ′ ℓ ′ ) 10 −9 [18][19][20], but more than forty orders of magnitude larger than the prediction of [1].
Therefore, the question we address in this letter is if the branching ratios of L → ℓℓ ′ ℓ ′ can really change so dramatically once one assumes the physical limit (PL), i.e. the hierarchy (PL) m νi ≪ P ≪ M W , instead of the ZML, and keeps full dependence of external momenta and masses in the loop diagrams.
Figure 1.One-loop diagrams contributing to the CLFV decay L − → ℓ − ℓ ′+ ℓ ′− in the unitary gauge: the Z penguins (a,b), the photon penguin (c) and the box (d).Wave function corrections (not depicted here) must be included as well.
In [2] it is argued that in the Z-penguin diagram shown in Fig. 1a there are two propagators of nearly massless fermions, which give rise to a log x i when the momentum q of the Z boson approaches q 2 = 0.This argument is supported by a computation of the Z-penguin as an expansion in [2] concludes that (q 2 /M 2 W ) log x i must dominate the branching ratio in the x i → 0 limit.This conclusion is flawed.First of all, in [2] the loop integrals written in terms of Feynman parameters are computed via a simple Taylor expansion of the denominator appearing inside.Such approximation is not legitimate for arbitrary values of the Feynman parameters, so it does not lead to an expansion of the integral itself.A proper asymptotic expansion of a Feynman integral can be obtained, for instance, via the expansion-byregions method [21,22], in which one divides the whole integration domain into various regions and then performs different Taylor expansions in each region.Only the sum of all regions' contributions eventually yields the desired asymptotic series.In addition to that, even if such q 2 /M 2 W expansion were performed correctly, the calculation presented in [2] implicitly assumes also the hierarchy q 2 ≪ m 2 νi .Therefore, the series expansion f 0 (x i )+ (q 2 /M 2 W )f 1 (x i )+ . . .does not reproduce the correct x i → 0 limit at fixed values of q 2 , since this limit lies beyond the validity range of q 2 ≪ m 2 νi .Recently, Ref. [23] presented a calculation of L → ℓℓ ′ ℓ ′ in the PL, in which the one-loop diagrams are numerically evaluated with full dependence on P. At variance with [2], they found branching ratios compatible with those in the ZML or smaller.However, the authors of Ref. [23] neglect the contribution from γ-penguins (as in Fig. 1c) and therefore their results are gauge dependent.Indeed, in processes with flavour changing neutral currents the gauge cancels entirely only in the sum of boxes, Z-and γ-penguins [24,25].Even if Ref. [1] retained only the logarithmic enhanced term x i log x i arising only from the Z penguins and the boxes, the omission of γ-penguins is not legitimate anymore as soon as one departs from this approximation.So we are still left with the doubt if the branching ratios in [23] are smaller as a consequence of calculating in the Feynman gauge or if there is a deeper physical meaning.
In this letter we present the decay widths of L → ℓℓ ′ ℓ ′ in the PL, fully analytic in M W , m ν and external momenta and masses.We compute them by making a systematic asymptotic expansion in P/M W and m ν /P of all Feynman diagrams by means of the expansion by regions.We will show that the neutrino mass dependence We will give an explanation of this exchange of mass scales in the logarithm by analysing the effective operators mediating the decay once the Z and the W bosons are integrated out.

DETAILS OF THE CALCULATION
Let us consider the SM extended with neutrino masses of either a Dirac or Majorana nature.The flavour eigenstates of the left-handed neutrino fields ν ℓL entering in the weak interactions become linear combinations of the three mass eigenstates ν i with masses m νi : where ν iL is the left-handed component of ν i and U is the PMNS matrix.The decay of a heavy lepton L = µ, τ into three lighter charged leptons ℓ, ℓ ′ = µ, e, with masses m L , m ℓ and m ℓ ′ , respectively, proceeds then via three classes of one-loop diagrams shown in Fig. 1: the boxes, the Z and γ penguins.We neglect diagrams with the exchange of a Higgs boson as they are further suppressed by two extra powers of 1/M 2 W due the Yukawa interaction.The partial width given in [1] was obtained in the ZML.In this approximation, the one-loop integrals depend only on the ratio x i .To leading order in x i , the amplitudes of the three classes of diagrams are:  I. Branching ratio for the CLFV decays L → ℓℓ ′ ℓ ′ in the ZML and the PL for normal ordering (NO) and inverted ordering (IO) of neutrino masses.The ratio between the two is also reported.In the ZML we assume m1 = 0 (NO) or m3 = 0 (IO).
where G F and α are the Fermi and fine structure constants, respectively, and sin 2 θ W is the sine of the Weinberg angle.Retaining only the terms enhanced by log x i ∼ 50, which appear in the boxes and Z penguins, one obtains the prediction for the rate in the ZML [1]: and ). Eq. ( 6) is obtained by taking the limit m 1 → 0 and assumes normal neutrino mass hierarchy, i.e. m 1 < m 2 < m 3 .For inverted mass hierarchy, the subscript '1' must be substituted with '3' and i = 1, 2. Eq. ( 6) also neglects subleading m ℓ,ℓ ′ /m L corrections from phase space integration.The values of the branching ratios in the ZML with normal and inverted mass hierarchy are reported in Tab.I. Current PDG values are employed for the lepton masses, neutrino mass splittings and neutrino mixing angles [14].
Let us now describe our calculation performed in the PL.We generated the complete set of diagrams in the Feynman gauge, and their relative counter-terms, using FeynArts [26] with a modified version of the SM file to account for Dirac neutrino masses and lepton flavour mixing.The amplitudes were reduced to one-loop tensor integrals using Form [27], via the FormCalc package [28], keeping the complete dependence on M W , P and m νi .The setup was independently checked by a second implementation based on FeynCalc [29].
Nowadays, lengthy expressions for the tensor integrals could be obtained in principle in an analytical form with full dependence on m νi , M W , m L , m ℓ,ℓ ′ and the invariants s ij = (p i + p j ) 2 , with p 1−3 the momenta of the three outgoing leptons, however their use is prohibitively cumbersome.It is therefore more helpful to compute them as series in the small parameters P 2 /M 2 W and m 2 νi /P 2 .To this end, we employed the method of expansion by regions (for an introduction see e.g.[22]).For all oneloop diagrams, we divided the integration domain into different regions and, for each region, we performed a Taylor expansion with respect to the parameters that are considered small there.Afterwards, by integrating every expanded integrand over the whole domain, and by summing the contributions from all the regions, we obtained the desired asymptotic expansion of the original one-loop diagram.The advantage of this method, compared for instance to an expansion of the full result after integration, is that the integrals arising in each region can be handled much more easily than the initial one, as typically they depend on just one or two mass scales.
We performed first an expansion assuming P ∼ m νi ≪ M W , without distinguishing at this point the two scales m νi and P. In a second step, the integrals arising from the first stage are further expanded in the limit m νi ≪ P. The total amplitude is then obtained by retaining from this expansion only the leading dependence on m νi , while higher order terms further suppressed by P 2 /M 2 W or m 2 νi /P 2 , or terms independent on m νi , are discarded.We performed several numerical checks at different stages of the calculation as a sanity check.To this end we took advantage of Mathematicas arbitrary-precision numbers and Package-X's analytic expressions of one-loop integrals [30], available in any kinematic configuration.We verified that our approximated expressions for the tensor integrals became increasingly accurate both by including higher order terms in the expansion as well as by taking the limit M W → ∞ and m νi → 0, at fixed values of P.

RESULTS
The partial widths are given by integrating the squared amplitude over the three-particle phase space of L → ℓℓ ′ ℓ ′ .These massive phase space integrals depend on two variables s ij , plus two or three masses of the external particles.By employing the expansion by regions one more time, we computed the phase space integrals as series in m ℓ ( ′ ) /m L , retaining only the leading terms in the final expressions for the rate.We obtain for normal neutrino mass hierarchy: where Eqs. ( 7) and ( 8) depend only on the neutrino mass splittings ∆m 2 ij and not on the value of the lightest neutrino's mass.The expressions with inverted neutrino hierarchy are obtained similarly as for the ZML.Higher order terms not included in Eqs. ( 7) and ( 8) are suppressed by m 2 L /M 2 W or m 2 νi /m 2 L .These corrections would arise by further expanding the squared amplitude.In addition, there are also subleading m ℓ /m L terms from the phase space integration.
At variance with the results presented in [2], Eqs. ( 7) and ( 8) are power suppressed by and yield values for the branching ratios of the order of 10 −55 , see Tab.I.Moreover, compared to Eq. ( 6) in the ZML, they do not have a logarithmic enhancement log 2 x i ∼ 2500.On the contrary in its place we get only log 2 x µ ∼ 176 or log 2 x τ ∼ 58, which are of comparable size with respect to other terms appearing in Eqs.(7) and (8).For this reason, the branching ratios in the PL turn out to be about one order of magnitude smaller than those in the ZML.
Note also that the presence of the singular terms log x ℓ or log x ℓ ′ is not in contradiction with the Kinoshita-Lee-Nauenberg theorem [31,32] and the cancellation of mass singularities for inclusive observables.In fact Eqs.(7) and ( 8) are valid strictly in the PL, i.e. when m νi ≪ m ℓ ( ′ ) .The limiting case of vanishing chargedmasses and non-zero neutrino masses violates the assumptions of our derivation and therefore is not a meaningful limit of our expressions.
In the case of purely Majorana masses (with no Dirac component), the neutrino weak eigenstates are still a combination of three Majorana fields, however the PMNS matrix takes on two new CP violating phases.These phases (1, α 21 , and α 31 ) multiply the three columns of the Dirac mixing matrix so that U ℓ1 → 1 • U ℓ1 , U ℓ2 → e iα21/2 U ℓ2 , and U ℓ3 → e iα31/2 U ℓ3 for all ℓ.In Eqs.(6)(7)(8) the PMNS elements show up in pairs as U * Li U ℓi , with i always the same in the two factors, so the Majorana phases cancel.Finally, we can understand the mechanism that converts the x i log x i in the ZML into a x i log x L in the PL by looking at the underlying effective theory arising after integrating out the Z and the W bosons.For simplicity, let us concentrate only on the operators associated with such logarithmic enhancement in the box 1d.Both in the ZML and in the PL, we can shrink the two W propagators to a point-like interaction and match the amplitude onto the the following dimension-six and dimension-eight operators: The first two operators correspond to the usual Fermi interaction mediating µ and τ leptonic decays.They contribute to L → ℓℓ ′ ℓ ′ via the one-loop diagram in Fig. 2a.The third operator in ( 9) is necessary to renormalize the effective theory, i.e. to cancel the UV divergence from the diagram 2a.These operators' Wilson coefficients are: We can imagine performing the matching between the SM and the effective theory at a scale µ = M W , and evolving the coefficients to a lower scale via the renormalization group.The coefficient C 8 explicitly depends on the renormalization scale µ and this dependence reveals the difference between the PL and the ZML.In the ZML, the evolution of C 8 can proceed down to a scale µ ∼ m νi , at which point we can integrate out the neutrinos and remove the operators O Lℓ ′ 6 and O ℓℓ ′ 6 which contain the neutrino field.We are then left with an effective theory with only O 8 , whose Wilson coefficient is frozen at C 8 (m νi ), i.e. it contains a log(M 2 W /m 2 νi ) (compare with Eq. ( 5)).On the contrary, in the PL, C 8 can run only until the scale µ ∼ m L is reached.In this case, all operators in (9) are still active at the scale m L , however C 8 produces only a milder log(M 2 W /m 2 L ) enhancement.Similar considerations can be applied as well to the Z penguin 1a.Therefore, the ZML overestimates the values for the branching ratios as it allows an unphysical evolution of these operators between M W and m νi , while in reality the running stops at the physical intermediate scale m L where the process happens.

CONCLUSIONS
Several experimental collaborations reported that the branching fractions of L → ℓℓ ′ ℓ ′ can be as large as 10 −14 , following the observation in [2] that the GIM cancellation for these decays is not so severe and takes the form of In this letter we showed that this conclusion is wrong.
We calculated and presented for the first time the branching ratios in the PL by performing a series expansion of all one-loop diagrams in the small parameters P/M W and m νi /P.Our fully analytic expressions prove that the GIM suppression in these decays is power-like | i U ℓi U * Li x i | 2 , similar to that found previously in [1] for the ZML, so that the claim from [2] must be rejected.We predicted the branching ratios in the Standard Model including neutrino masses to be of the order of 10 −55 , even smaller than those obtained in the ZML, as the latter prediction contains an unphysical logarithmic enhancement.
In the end, we remark that since the GIM suppression is solely governed by the underlying effective description of the process, i.e. the hierarchy of the internal mass scales and the convergence properties of Feynman integrals, external momentum effects could not have affected the GIM cancellation to such a large extent, as claimed in Ref. [2], compared to the finding in the ZML.

Figure 2 .
Figure 2. Example of two diagrams mediating L → ℓℓ ′ ℓ ′ in a low energy effective field theory description.