Flavor violating leptonic decays of $\tau$ and $\mu$ leptons in the Standard Model with massive neutrinos

We have revisited the computations of the flavor violating leptonic decays of the $\tau$ and $\mu$ leptons into three lighter charged leptons in the Standard Model with non-vanishing neutrino masses. We were driven by a claimed unnaturally large branching ratio predicted for the $\tau^-\to \mu^- \ell^+ \ell^-$ ($\ell=\mu, e$) decays \cite{Pham:1998fq}, which was at odds with the corresponding predictions for the $\mu^-\to e^- e^-e^+$ processes \cite{Petcov:1976ff}. In contrast with the prediction in \cite{Pham:1998fq}, our results are strongly suppressed and in good agreement with the approximation done in ref.~\cite{Petcov:1976ff}, where masses and momenta of the external particles were neglected in order to deal with the loop integrals. However, as a result of keeping external momenta and masses in the computation of the dominant penguin and box diagrams- we even find slightly smaller branching fractions. Therefore, we confirm that any future observation of such processes would be an unambiguous manifestation of new physics beyond the Standard Model


Introduction
Lepton flavor violating (LFV) processes are forbidden in the standard model (SM) [1] with massless neutrinos. However, the experimental evidence of neutrino oscillations [2] claims for an extended model with neutrino mass terms. For massive neutrinos, the mass matrix will be nondiagonal in the interaction (weak) basis, as occurs in the quark sector [3], and the mixing of three light neutrinos could be described through the 3 × 3 unitary Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix [4]. In such scenario, charged LFV transitions could arise, for instance, from one loop diagrams involving a couple of W ν vertices with different flavor neutrinos each. However, it turns out natural having a strong suppression for this class of processes owing to a GIM-like mechanism [5], just as it has been reported for the µ − → e − γ decay, with a prediction at an unobservable low rate: BR(µ − → e − γ) ∼ O(10 −55 ) [6][7][8], which is far away from the capacity of any current or foreseen experimental facility.
By way of contrast, the prediction for the τ − → µ − + − ( = µ, e) decays given by ref. [10] indicates that the GIM cancellation for these processes is much milder and a value of BR(τ − → µ − + − ) ≥ 10 −14 is reported. An updated evaluation using the amplitude derived in ref. [10], employing the latest global fit results for neutrino mixing [11,13] yields a branching fraction ∼ 4 · 10 −16 for the three muon channel. Both values are still far away from the PDG upper bounds, 1.5·10 −8 (for = e) and 2.1·10 −8 ( = µ) at 90% confidence level 1 . Similarly, we verified that using the values reported in refs. [11,13] for the neutrino mixing parameters, Pham's result [10] would predict a µ − → e − e + e − branching ratio of ∼ 2 · 10 −21 , larger than Petcov's prediction (∼ 10 −53 evaluated with updated neutrino masses and mixings input) [8] by at least some thirty orders of magnitude. Again, the current upper limit on this decay channel (1 · 10 −12 at 90% C.L. [11]) is still far from testing Pham's result [10]. This author claims that this unexpectedly large estimation is due to the presence of a divergent logarithmic term depending on the neutrino mass, which comes from a one-loop diagram that involves two neutrino propagators (diag. (d) in our fig. 1).
Certainly, considering effects or processes that arise from quantum corrections could involve divergent loop integrals. However, in any renormalizable theory, the possible divergences must vanish order by order (in the loop or effective field theory expansion) to be able to define (finite) observables. In fact, in a QFT the divergences can be classified into two types: ultraviolet (UV) and infrared divergences (IR). The former (UV) appear in the high-energy regime and they can be healed redefining the theory parameters, whereas the latter (IR) occur in the low-energy regime and can be classified in soft and collinear divergences, which cancel however in properly defined (IR-safe) observables [15]. We show that the seeming logarithmic divergent behavior of the LFV amplitude reported in ref. [10] is not present, as the vanishing momentum transfer approximation considered in that paper lies outside the physical region. Consequently, the rates of L − → − − + decays in the SM extended with massive neutrinos are extremely suppressed, in agreement with ref. [8]. It is worth noting that the LFV amplitudes must vanish in the limit of massless neutrinos. This requirement is satisfied by the result of Ref. [8], but it is not the case in ref, [10] which behaves as j U Lj U * j log(m j /m W ) for very small neutrino masses. Our result, as it will be shown below, satisfies the expected agreement with the SM.
In section 2 and 3 we discuss in detail our computation of these processes and compare it to those in refs. [8] and [10], showing explicitly why the approximation in [10] is unreliable, and reproducing the results of [8] in the approximation where masses and momenta of the external particles are neglected from the beginning. However, we also analize the numerical accuracy of this approximation. Finally, we state our conclusions in section 4. Several appendices complete technical details of our calculation.

Z-Penguin contribution emission from internal neutrino line
The L − → − − + decays can be induced through the diagrams depicted in fig. 1. Since the main purpose of this work is to falsify the existence of the logarithmic divergent term claimed in ref. [10], we first concentrate on the amplitude of the diagram (d). We have, however, verified the corresponding expressions for the loop integrals in ref. [8] for the particular process µ − → e − e − e + , when masses and momenta of external leptons are neglected in the computations. Particularly, in Ref. [8] it is shown that the corresponding branching ratio is completely dominated by those diagrams with two neutrino propagators, i. e. (d) and (e) in fig. 1, which contribute comparably.
In our analysis, we keep employing the convention used by ref. [10], in order to denote the masses and momenta (see fig. 1) of the external leptons, that is M and P (m and p) stand for the mass and momentum of the L − ( − ) lepton, respectively. In this way, the amplitude of the diagram (d) can be written as )v(p 2 ) 2 is independent of the loop integration, whereas the relevant part for the latter is given by the effective ZL transition as follows: where U im are entries of the PMNS mixing matrix. In the Feynman-'t Hooft gauge we have After making the loop integration using dimensional regularization in order to deal with the (logarithmic) UV divergences, the Lorentz structure for the Γ λ j factor can be written as follows, .., f ) are functions given in terms of the momentum transfer q 2 , and the neutrino mass squared (of course F k functions will also depend on the mass of the W gauge boson and external masses, but these have well-defined values).
At this point, it is worth to note that in the approximation where the momenta of the external particles are neglected in equation (2.3), such as it is done in ref. [8] for the µ → 3e decay, the computation is simplified considerably, as the only possible contribution is given by the F 0 a function, where we are using a superscript 0 in order to distinguish this approximation. In this simple case, the F 0 a function will not depend on q 2 and is given in terms of the Feynman parameters as follows Whereas in terms of PaVe functions it is given by Now, one analytical expression for the F 0 a function can be obtained in a straightforward way either integrating over the Feynman parameters in eq. (2.5) or using the definition of the B 0 (0, m 2 , m 2 ) scalar function in eq. (2.6). In such a way that after making an expansion around m 2 j = 0 we obtained From eq. (2.7) it turns clear that, in this approximation, the amplitude will be proportional to the neutrino mass squared, where the dominant contribution, due to the big gap between the neutrino and W boson mass scales, comes from the first term as it involves a relative factor log compared to the second one 3 , whereas the independent terms on neutrino mass will vanish by the GIM-like mechanism. Therefore, the structure of the matrix element for the contribution of the diagram (d) in fig. 1 in the approximation where masses and momenta of the external particles are neglected is given by where we have defined We verified that eq. (2.8) reproduces the result reported in ref. [8] considering only the first term in eq. (2.7) and the simple case of two families.
Returning to the general case (non-zero masses and momentum of the external particles), we also have obtained the F k functions using both Feynman parametrization (we will denote the corresponding expressions by F F k ) and the Passarino-Veltman (PaVe) technique (denoted by F P V k ) [16,17], employing FeynCalc [18]. In particular, we agree with the expressions previously reported in ref. [10] in terms of the Feynman parameters 4 , namely the F F k functions can be written as and D j is defined as (2.14) We have omitted in f a the term associated with the UV divergence since it is independent of m j and vanishes owing to the GIM-like mechanism.
On the other hand, the F k functions in terms of the PaVe scalar functions are given as follows where λ is the Kallen function λ(x, y, z) = x 2 + y 2 + z 2 − 2(xy + xz + yz), and the ξ k factors can be found in the appendix A. 5 Unlike the approximation made in ref. [8], the presence of masses and momenta of the external particles in the computation hinders the way for the derivation of analytical expressions for the integrals in eqs. (2.10) or (2.15) 6 . Nevertheless, we have done a numerical cross-check between both expressions, where we have employed the Looptools package [19,20] for the evaluation of the PaVe functions and a numerical Mathematica [21] routine for the evaluation of the parametric integrals (see fig. 2). We have found an excellent agreement between these two expressions for values of q 2 < 4m 2 j , which are, however, out side of the physical domain for the considered decays, since q 2 min = 4m 2 m 2 j . In this way, owing to the simpler integrals, we verified that a better precision is found in terms of PaVe functions than using Feynman parameters, this feature is illustrated, as an example, for the particular case of the Zτ µ transition in fig. 2 for the (dominant, as we will show) F a factor. 5 The cancellation of the UV divergences for the Fm functions in terms of the PaVe functions occurs again by the GIM mechanism. This can be verified easily owing to the fact that the sums over the coefficients of the different scalar B 0 functions, which contain an isolated divergent term, are independent of m j . That is 6 The analytical expressions of the first integrals over the y parameter in eqs. (2.11), (2.12) and (2.13) can be derived from the formulas reported in appendix B.
Black dashed line stands for the numerical evaluation of the Fa function in terms of the Feynman parameters depicted by FF a (eq. (2.10)), whereas the red line corresponds to the evaluation in terms of the PaVe functions represented by FP Va (eq. (2.15)). We have found some numerical instabilities for the evaluation of the FF a function in the region 0.01 GeV< mj < 0.1 GeV. On the other hand, a better precision is achieved in the evaluation of the FP Va function with the help of the Looptools package. In order to perform a comparison with the approximation done in ref. [8], we also show the complete At this point, we want to stress that we disagree with the approximation done in ref. [10] in order to estimate the relevant dependence on the neutrino mass of the F k functions. We highlight that we are studying a process where the momentum transfer q 2 must be non-vanishing and in principle is much larger than the neutrino squared mass, m 2 j , which comes from the loop computation. Therefore, using an expansion around q 2 = 0 in order to simplify the integration over the Feynman parameters keeping the terms proportional to m 2 j in the denominators of equations (2.11), (2.12) and (2.13), as it is done in ref. [10], modifies substantially the behavior of the original functions in the interesting physical region for the neutrino masses and, as a consequence, it gives rise to an incorrect infrared logarithmically divergent behavior of the F k functions when m j goes to zero, without any possible cure. In particular, the dependence on the momentum transfer, q 2 , plays a crucial role in the behavior of the F k functions. In this respect, we point out the presence of a small imaginary part in the F a function, which emerges for the physical values 4m 2 j < q 2 . As we mentioned before, the q 2 minimum in the L − → − − + decay is given by 4m 2 , which is much larger than neutrinos masses. This, together with the difficulties in obtaining analytical expressions directly for the F k functions suggests employing some numerical approximation to deal with the problem. Because of this, we approximate the F k functions in the physical region for the neutrinos masses by fitting the curves for the real and imaginary parts of the F k functions evaluated in terms of the PaVe function 7 . We have found a reasonably good fit of the form 7 Our fits for the F k functions are taken with the precision of the Looptools package considering a neutrino mass varying from 10 −15 GeV to the benchmark point mµ (me), for a fixed value of q 2 = 4m 2 µ (q 2 = 4m 2 e ) for the Zτ µ (Zτ e and Zµe) vertices.
where u = 1 for k = a, b and u = M for k = c, d, e, f and the respective values for the Q k = Q R k + iQ R I and R k = R R k + iR R I factors of all considered channels are given in appendix D. From eq. (2.19), it turns clear that the Q k factors will not contribute owing to the GIM-like mechanism, whereas the relevant contributions is given by the R k factors. Then, according to our numerical results, we find that the relevant factors of the F b , F c and F d functions are suppressed with respect to the F a factor. On the other hand, despite the respective factors of F e and F f functions are larger than those of the F a function, when the momentum transfer becomes smaller and smaller their helicity suppression makes them negligible. Therefore, we will concentrate on the contribution of the F a function. Furthermore, in order to justify our results, we have made an expansion for the PaVe functions involved in eq. (2.18), following the same strategy that Cheng and Li for the µ → eγ decay [6], that is expanding the loop integrals around m 2 j = 0 (more details of our expansions are given in appendix E), and with the help of Package-X program [24], we have been able to rewrite the F P Va contribution as follows, , and the f Q and f R factors can be found in the appendix E. We verified that our numerical fits for the Zτ µ and Zτ e vertex are in a very good agreement with eq. (2.22), whereas a deviation is found for the Zµe vertex, as can be seen in Table 9, we consider the results obtained from eq. (2.22) for the effective vertices as our reference ones.
In this way, we can approximate the amplitude for diagram (d) according to eq.
Now, in order to evaluate the respective branching fractions for the L − → − − + decays we considered the state of the art best fit values of the three neutrino oscillation parameters [11,13].
Neglecting for the moment the box contributions, we get the branching fractions reported in table 1.

Decay channel Our result
Ref. [8] Table 1. Branching ratios for the L − → − − + decays (neglecting the box and the subdominant penguin contributions with only one neutrino propagator), which are obtained using the current knowledge of the PMNS matrix. The last column values correspond to the approximation where external masses and momenta are neglected [8]. Our results are smaller than those by around one (two) orders of magnitude for the µ (τ ) decays.

Contributions of the box diagrams
Now, in order to make a complete comparison with the approximation done in ref. [8] we have also obtained the amplitude for the box diagram (e) in fig. 1. Note that unlike the penguin diagram (d), which involves two neutrino propagators of the same flavor, the box diagram (e) can involve two neutrino propagators with different flavors. Thus, in full generality, the amplitude can be written as follows where we defined and the relevant loop integral is given by (see fig. 1 (e)) Since we have written the equation (3.3) in terms of P , p 1 and p 2 momenta the integral must take the form The H k factors depend upon the kinematical variables s 12 = (p 1 + p 2 ) 2 = q 2 and s 13 = (p 1 + p) 2 , in addition of m i and m j . Anew, in the approximation where momenta of the external particles are neglected in eq. (3.3), the only contribution is given by the H 0 a function, which will not depend either on s 12 or s 13 . In such case, we obtained the following simplified expression Whereas, in terms of PaVe functions, H 0 a reads In the same way that F 0 a form factor, an analytical expression for H 0 a can be obtained easily from either eq. (3.5) or eq. (3.7). This time, making a double Taylor expansion, first around m 2 i = 0 and then around m 2 j = 0, we obtained that (3.8) , the amplitude -in this approximation-is given by with Again, we verified that taking into account the first term in eq. (3.8) and considering only two families, eq. (3.9) reproduces the expression reported in ref. [8] for the amplitude of the box diagram 1 (e) in the µ → 3e decay.
In the general case, we also obtained the H k (k = a, b, ..., j) functions in terms of both Feynman parameters integrals, H F k , and PaVe functions, H P V k . This time, the H k functions will depend on the squared masses of two different neutrinos, m 2 j and m 2 i , and on two independent phase space variables s 12 and s 13 . Using Feynman parametrization these functions read where (3.13) In the previous expressions, the denominator function is given by (3.14) Expressions are rather lengthy in terms of the PaVe functions, so that here we only present the expression for the dominant H a function, which can be written as and where χ k factors are reported in Appendix A.
As far as the general case is concerned, we can see that although there are additional contributions associated with the H k functions, with k = b, c, d, . . . j; they are expected to be suppressed, as they correspond to higher-dimensional operators, with respect to the H a function associated with a (V − A) × (V − A) operator. Therefore, we will concentrate on the H a function in order to estimate the box diagram contribution. We also have done a numerical cross-check between the expressions for the H a function given in terms of the Feynman parameters eq. (3.11) and the PaVe functions eq. (3.15), as can be seen in fig. 3. In this case, it turns very complicated and far away of the purpose of this work to obtain an analytical expression for the H a function in eq. making an expansion for the respective scalar PaVe functions, owing to the number of propagators involved and the dependence on two different neutrino masses. However, we can expect a good approximation through our numerical results, such as occurs with the penguin contribution. Thus, we estimate the relevant dependence on the neutrino mass for the H a function fitting the curve for the real and imaginary parts of the H a function evaluated in terms of the PaVe functions considering fixed values for the m i , s 12 , and s 13 parameters 10 . We obtained a good fit of the form where R Ha ≈ 1.5 + i0.007, for all different τ → − − + channels, whereas R Ha ≈ 1.5, for the µ − → e − e − e + channel. These numbers were obtained considering that ∆m 2 ij = 10 −3 eV 2 , and representative values for s 12 and s 13 within the corresponding phase space. Now we can evaluate the branching ratios for the L − → − − + decays using the previous results. We will first make a partial evaluation neglecting the penguin contributions (only box diagrams are considered), which yields the values in table 2.
Our final results, where the dominant penguin and box contributions are considered, are collected in table 3, where they are compared to those obtained using Petcov's results [8] with updated input. Our predictions are even smaller than Petcov's updated results, as a consequence of keeping external masses and momenta in our computations.

Conclusions
We have revisited the L − → − − + decays in the SM with massive neutrinos. We obtained expressions in terms of both Feynman parameters and scalar Passarino-Veltman functions for the relevant loop integrals of the (dominant) diagrams that involve two neutrino propagators considering non-vanishing masses and momenta of the external particles. Opposed to the previous calculation reported in ref. [10], we found that all the different amplitudes for these processes are strongly suppressed (as they are proportional to the neutrino mass squared). In the particular case of the penguin contribution with two neutrino propagators, we highlight that it is crucial to save the dependence on the momentum transfer in the Feynman integrals in order to evaluate the amplitude in the physical region for the neutrino masses. This fact avoids the incorrect divergent logarithmic behavior in the amplitude claimed in ref. [10]. As far as the box contribution is concerned, we found that the dominant term comes from H a function that is associated with a (V-A)×(V-A) operator, and it is in good agreement with the approximation done in Ref. [8].
Current and forthcoming experiments were approaching the limits predicted by ref. [10] on the SM prediction for the lepton flavor violating τ − → µ − + − ( = µ, e) decays due to non-zero neutrino masses. This prediction was at odds with ref. [8] corresponding computation for the µ − → e − e + e − decays predicting an extremely suppressed, unmeasurable branching ratio (as in L − → − γ processes). The most important result of our analysis is the confirmation (in agreement with ref. [8]) that any future observation of L − → − − + decays would imply the existence of New Physics.

A One-loop PaVe scalar functions
In this appendix we collect the ξ ij i=a,...,f ;j=0,..,5 factors entering our results in eq. (2.18): As far as the χ k factors entering the H P Va functions in eq. (3.17), they are given as follows

B Some useful integrals
As we mentioned in the text, analytical expressions for the double integrals in eqs. (2.11), (2.12) and (2.13) are not easy to obtain. However, the first integrals over the y-Feynman parameter can be derived from the following expressions where we have defined the functions that follow Because of the necessity of antisymmetrizing the amplitude when = , the total contribution for the sum of the penguin and box diagrams in this case is given by On the other hand, when = , there is only one penguin diagram since the neutral Z boson does not change flavor. Besides, we have to add the box diagram interchanging − (p) ↔ − (p 1 ). Therefore, we have where β Ha has been defined in the main text and In the Petcov's approximation, taking only the dominant term and since the contribution of the penguin and box diagrams have opposite sign, the dominant terms are given by the second terms in eqs. (C.1) and (C.2), respectively. Therefore, |M 2 | is given by The unpolarized differential decay width for the L − (P ) → − (p) − (p 1 ) + (p 2 ) decays is given by where N is the number of identical particles in the final state and s 12 = (p 1 + p 2 ) 2 = q 2 and s 13 = (p 1 + p) 2 . The corresponding integration limits are given by The numerical values for the Q k and R k factors involved in of our fits for the Zτ µ, Zτ e and Zµe effective vertices are given as follows Zτ e (q 2 = 4m 2 µ )   Table 7. Same as Table 6 but considering q 2 = 4m 2 e . E Expansion of the PaVe functions around m 2 j = 0 The scalar PaVe functions involve in eq. (2.18) calculation are defined as follows , where p 2 = m 2 , P 2 = M 2 and q 2 = (P − p) 2 = m 2 + M 2 − 2P · p.
Although there are infrared divergences ∆ I on the eqs. (E.12, E.13, E.15, E.16), the factor R a is free of them. Further, there is no dependence on the renormalization scale, and these results are in agreement with our numerical fits. Taking into account the imaginary part of the C 0 (m 2 , M 2 , q 2 , 0, m 2 W , 0) function, it is possible to derive analytically that the imaginary parts appearing in the last column of Table 9 are exactly π.