Next-to-leading order prediction for the decay μ→ee+e−νν¯$$ \mu \to e\kern0.22em \left({e}^{+}{e}^{-}\right)\;\nu \kern0.2em \overline{\nu} $$

We present the differential decay rates and the branching ratios of the muon decay with internal conversion, μ→ e (e+e−) νν̄, in the Standard Model at next-to-leading order (NLO) in the on-shell scheme. This rare decay mode of the muon is among the main sources of background to the search for μ → eee decay. We found that in the phase space region where the neutrino energies are small, and the three-electron momenta have a similar signature as in the μ→ eee decay, the NLO corrections decrease the leading-order prediction by about 10− 20% depending on the applied cut. ar X iv :1 61 1. 03 72 6v 2 [ he pph ] 6 F eb 2 01 7


Introduction
Lepton flavour is not a conserved quantity in nature, indeed neutrino oscillations have unveiled that the Standard Model (SM) must be modified to include neutrino masses and mixing. Processes with charged lepton flavour violation (CLFV) can occur, therefore, through neutrino mixing in the loops; however, their rates are extremely small because of the suppression given by neutrino masses. For example the µ → eγ branching ratio, which is proportional to (m ν /M W ) 4 , is estimated to be at the level of 10 −50 or smaller, far beyond the sensitivity of any foreseeable experiment. For this reason, an experimental observation of CLFV would be a bright evidence of new physics (NP) beyond the SM.
Muon and tau decays with CLFV have been studied in a model independent way in the framework of higher dimensional operators [1][2][3] and in a general effective field theory description of the weak interactions at low energies [4]. Many scenarios of NP introduce additional sources of mixing between the lepton families that can easily lead to strong CLFV contributions. CFLV violation can be realized, for example, in minimal see-saw type extensions of the SM [5,6], in the MSSM via flavour non-diagonal SUSY breaking terms [7][8][9][10][11], or in the context of two-Higgs-doublet models [12][13][14][15] and composite Higgs models [16] with generic flavour structures in the lepton sector.
Despite the tau lepton has the advantage of a mass greater than the muon one and thus, from the theoretical point of view, a better sensitivity to NP effects, muons can be produced in a much larger quantity and can be measured with much better sensitivity. In the muon sector, the upcoming Mu2e [17], DeeMe [18] and Comet [19] experiments will search for µ → e conversion of muons bounded to a nucleus, while Meg [20] and Mu3e [21] at PSI are dedicated to the SM forbidden decays µ → eγ and µ → eee, respectively.
The current limits on the µ → e transitions are very stringent due to the constraints from Meg and Sindrum collaborations: B(µ + → e + γ) < 4.2 × 10 −13 90% C.L. [20], (1.1) B(µ − → e − e + e − ) < 1.0 × 10 −12 90% C.L. [22]. (1.2) These bounds will be further improved in the future: Meg-II is planned to measure the µ → eγ branching ratio down to 4 × 10 −14 [23] while Mu3e is expected to reach the level of 10 −16 for µ → eee [21,24,25]. Such an unprecedented precision requires an equivalent level in the control and suppression of the background. Radiative muon decay, µ → eγνν, constitutes an important source of background to µ → eγ searches; secondly it provides a tool for calibration, normalization and quality check of the experiment [26]. It was studied at next-to-leading order (NLO) accuracy in [27,28]. The background in µ → eee searches originates from the accidental coincidence of normal muon-decay electrons and positrons, that within the detector resolution show the characteristics of the decay signal, and from the muon decay with internal conversion, µ → e (e + e − ) νν, which is indistinguishable from the signal except for the energy carried out by neutrinos. This SM background can be suppressed only via an excellent momentum resolution and a precise reconstruction of the three-electron total energy, which must be as close as possible to the muon mass.
While the SM decay rate of µ → e (e + e − ) νν was studied at the leading order (LO) in [29][30][31][32], radiative corrections are currently missing in the literature. The precision goal of Mu3e experiment and the expected momentum resolution, about 0.5 MeV [24,25], call for the estimate of the NLO corrections, especially in the narrow region of the phase space where the missing energy is small and the three-electron momenta have a similar signature as the µ → eee decay. It is the aim of this paper to present the first calculation of such NLO corrections. 1 We begin our analysis in section 2 reviewing the SM prediction for the differential decay rate at LO. The technical ingredients employed in our calculation of virtual and real corrections are presented in section 3. Our NLO predictions for the branching ratios are reported in section 4, where we discuss also the impact on CLFV searches. Conclusions are drawn in section 5.

Conventions and LO decay rate
In this section we introduce our conventions and we discuss the tree-level decay rate. In the SM the decay of a muon into three electrons and two neutrinos proceeds through the emission of an off-shell photon with subsequent internal conversion into e + e − , as shown in figure 1. Since we are interested only in the leading contribution in G F , photon radiation from the W boson and electron-pair production from heavy bosons will not be considered in the following.
Let us consider specifically the decay of a negative muon: In the SM, the tree-level decay rate for an unpolarized muon is, in its rest frame, [34] is the Fermi constant, defined from the muon lifetime, and α = 1/137.035 999 157 (33) is the fine-structure constant [35,36]. Calling m µ and m e the masses of the muon and the electron (neutrinos are considered to be massless), we define the ratio r = m e /m µ . The four-vectors P , p 1 , p 2 and p 3 are the momenta of the muon, the two electrons and the positron, respectively. We denote with p 123 = p 1 + p 2 + p 3 the sum of the electrons and the positron momenta and with m 2 123 = p 2 123 their invariant mass squared. Also, we define m 2 12 = (p 1 + p 2 ) 2 and t 2 = (t 1 + t 2 ) 2 = (P − p 123 ) 2 the invariant masses squared of the two electrons and the two neutrinos, respectivelyt 1 and t 2 are the neutrino and anti-neutrino momenta. Note that the dependence on the (undetected) neutrino momenta has been integrated out analytically in (2.2).
Solid angles and momenta labeled with the superscript ' * ' are in the center-of-mass system (c.m.s.) of the two electrons and the positron, where p * 123 = (m 123 , 0), while those with ' * * ' are in the c.m.s. of the two electrons, where p * * 12 = (m 12 , 0). The z-direction of the solid angles dΩ * 3 and dΩ * * 1 are given by the direction of p 123 and p 12 = p 1 + p 2 , respectively. The dimensionless quantity G LO is a rational function proportional to the Born squared matrix element -it is thus Lorenz invariant -and it depends on all possible scalar products built with the momenta P, p 1 , p 2 and p 3 . We introduced such a function in order to factorize out the physical constants α, G F and m µ from the squared matrix element. Two indistinguishable e − are present in the final state; the symmetry property of the matrix element assures that G LO and the decay rate (2.2) are symmetric under the exchange The (lengthy) explicit expression of eq. (2.2) is provided as an ancillary file of this paper. For further details, we refer the reader to our Appendix.

NLO corrections: details of the calculation
In this section we will consider the SM prediction for the differential decay rate of (2.1) at NLO in α. Virtual and real corrections are evaluated in the Fermi V -A effective theory of weak interactions: The Fermi Lagrangian is where ψ µ , ψ e , ψ νµ , ψ νe denote the fields of the muon, the electron and their associated neutrinos, respectively; P L = (1 − γ 5 )/2 denotes the left-hand projector operator. Under this approximation tiny term of O(αm 2 µ /M 2 W ) ∼ 10 −8 due to the finite W -boson mass are neglected -they are even smaller than the NNLO corrections of O(α 2 ).
A Fierz rearrangement of the four-fermion interaction (3.2) allows us to factorize the amplitudes of virtual and real corrections into the product of spinor chains depending either on the neutrino momenta or on the muon and electron ones. In this way the neutrino phase space integral can be carried out analytically, as was done at tree level, so as to decrease by two units the dimensionality of the integral that must be performed numerically.

Virtual corrections
The one-loop amplitudes, shown in figure 2, are reduced to tensor integrals and subsequently decomposed into their Lorentz-covariant structure by means of the algebra manipulation program Form [37] and the Mathematica package FeynCalc [38,39]. For the numerical evaluation of the tensor-coefficient functions we employed the LooopTools library [40,41].
Ultraviolet (UV) divergences are regularized via dimensional regularization; UV-finite results are obtained by renormalizing the theory (3.1) in the on-shell scheme. Indeed, as shown long ago by Berman and Sirlin [42], to leading order in G F , but to all orders in α, the radiative corrections to muon decay in the Fermi V -A theory are finite after fermion mass and wave function renormalization. A small photon mass λ is introduced to regularize the infrared (IR) divergences, while the finite electron mass m e regularizes the collinear ones.
The contribution to the rate coming from the hadronic vacuum polarization, which is not calculable at low energy in perturbative QCD, is quite small in the muon case since the invariant mass of the electron-positron pair never exceeds the pion threshold. However this kind of correction starts to be relevant if one considers, instead of the muon, the tau decays τ → ( + − )νν. These effects can be taken into account expressing the hadronic vaccum polarization, Π had (q 2 ), in terms of e + e − → hadrons cross section data: The normalization factor 4πα(s) 2 /(3s) is the tree-level cross section of e + e − → µ + µ − in the limit s 4m 2 µ -note that σ(e + e − → hadrons) does not include initial state radiation or vacuum polarization corrections. The optical theorem connects R had (s) to the imaginary part of hadronic vacuum polarization: The vacuum polarization can be then obtained by means of the dispersion relation. In this work, we made use of the package alphaQED [43][44][45][46] for the evaluation of the functions R had and Π had . Figure 2. One-loop diagrams for µ − → e − (e + e − ) ν µνe decay. The muon and the electrons are drawn with bold and thin lines, respectively; the dots symbolize the Fermi interaction (for simplicity the neutrinos are not drawn). For each diagram, a symmetric one with p 1 ↔ p 2 interchanged must be considered.

Real photon emission
The rate of the bremsstrahlung process, the decay (2.1) where an additional photon is produced, blows up when the photon energy becomes small, leading in the phase space integral to the well-known logarithmic IR singularity. In order to handle the IR singularity we adopted a phase-space slicing method: we introduced a small photon energy cut-off ω 0 and we divided the real emission contribution into a soft and a hard part.
The soft part, which contains the IR singularity, comes from the phase space region where the photon energy is below the threshold ω 0 . As in the case of the virtual diagrams, the IR singularity is regularized by a small photon mass λ. By taking advantage of the soft photon approximation for the amplitude of the bremsstrahlung process, it is possible to perform the integral with respect to the photon momentum analytically; the processindependent result, which was derived in [47] (see also ref. [48]), depends only on the charges and momenta of external particles in the corresponding Born process. We checked that the IR divergence cancels out once the soft part is added to the one-loop diagrams.
Although the addition of the soft photon emission to the virtual corrections is sufficient to obtain a finite differential width, in general it is not adequate for real experiments since they cannot provide a photon energy threshold ω 0 small enough for the validity of the soft photon approximation, which neglects terms of order ω 0 /m µ . Therefore it is necessary to include the hard part as well, i.e. the contribution to the rate due to photons with energy greater than ω 0 . The soft and the hard parts must be properly merged to assure the ω 0 -independence of the final physical observables: the numerical integration of hard part yields a log ω 0 /m µ enhancement that must cancel against the explicit ω 0 -dependence in the soft part for sufficiently small values of ω 0 . The value of ω 0 can be fixed once the prediction for the NLO corrections reaches a sort of "plateau", i.e. when a further reduction of ω 0 cannot be resolved anymore within the numerical error.

Branching fraction
The branching ratio B of the µ − → e − (e + e − ) ν µνe decay can be obtained integrating the differential decay rate (2.2) over the allowed kinematic ranges, and multiplying it by the muon lifetime τ µ = 2.1969811 (22) × 10 −6 s [49]. We recall that the Particle Data Group (PDG) defines the Fermi constants of weak interactions from the muon lifetime evaluated in the Fermi V -A theory [49]; its definition is given by where F (x) = 1 − 8x + 8x 3 − x 4 − 12x 2 ln x is the phase space factor while δ µ incorporates the QED correction evaluated in the Fermi V -A theory: the corrections of virtual and real photons up to O(α 2 ), as well as the contribution of the decay (2.1) at tree level [50][51][52][53][54][55][56][57][58]. Note that it is possible to make use of eq. (4.2), instead of the experimental value of τ µ , for the normalization of width; in such way the dependence on G F and m 5 µ is removed from the branching ratios.
The analytic integration over the kinematic ranges (4.1) of the LO differential rate in eq. (2.2) yields [31] (1 + ln 2) ln 2 2 where with I(r 2 ) = 9.47056, and the kernel function K(u) is given in eq. (11) of ref. [57]. As discussed in [31], the integral (4.4) behaves like ln n r for non-negative integers n ≤ 3, so that it is singular in the limit m e → 0. It contains also vanishing terms in that limit, but since the original calculation of K(u) neglected them, the result in eq. We present in the first row of table 1 the LO and NLO branching ratios, the O(α) correction coming from virtual and real diagrams, denoted with δB NLO , and the K-factor, which is the ratio between the NLO and the LO prediction. They are computed taking into account the full dependence on the mass ratio r. The numerical integrations were performed with Monte Carlo methods by means of the Cuba library [59]; the results were tested with different numerical integration methods. In table 1 the uncertainty due to numerical errors is labeled with the subscript "n". Moreover we also estimated the theoretical uncertainty associated to the renormalization scale variation; they are denoted with the subscript "µ". It is quantified by converting the renormalization scheme for α from the on-shell scheme to the hybrid MS adopted in the NNLO calculation of the muon lifetime [54,57]. In this scheme the electron loop in the photon vacuum polarization is renormalized in MS while all other fermion loops in the on-shell scheme. The quoted uncertainty is the difference between the NLO prediction evaluated at the renormalization scales µ = m e , corresponding to the on-shell scheme, and µ = m µ .
It is interesting to note that a Born-virtual interference term yields a null contribution after phase space integration. Indeed, let us consider the one-loop diagrams where the positron-electron pair is produced by two virtual photons, corresponding to the boxes and the pentagons in the last two rows of figure 2. If we regard as a cut diagram the intereference between these one-loop amplitudes and the tree-level ones in figure 1, we recognize that there is a closed fermion loop attached to three photon lines. The Furry theorem assures that the contribution after phase space integration is zero. However, such interference term cannot be neglected in the differential decay rate, in fact the cancellation happens between couples of phase space points related by an exchange of the p 1 and p 3 momenta (or p 2 and p 3 ). We explicitly verified that this interference term vanishes within the numerical error but we neglected it in our Monte Carlo integration in order to speed up the numerical convergence.
The branching ratio of (2.1) was measured long ago by the Sindrum experiment [60], This measurement agrees with our theoretical prediction in table 1; new more precise results are expected in the future by the Mu3e experiment [21].

Impact on CLFV searches
The relative magnitude of radiative corrections were also studied in the specific final-state configuration of the decays (2.1) where the neutrino energies are very small and the total energy of the three electrons is close to m µ . As already mentioned in the introduction, this phase-space region is of particular interest to µ → eee searches because the muon decay (2.1) can mimic the three-body decay mode with CLFV. The upper panel of figure 3a shows dB/dm 123 , the normalized NLO differential rate as function of the three-electron invariant mass m 123 , close to the end point region m 123 = m µ . The local K-factor is drawn in the lower part. The rate, evaluated at fixed value of m 123 , is fully inclusive in the bremsstrahlung photon.
Beside that, we calculated also the branching fraction applying a cut on the missing energy, in analogy to the analysis done in ref. [32] at LO. Let us define B( / E max ) to be the integral of the differential decay rate over the phase space region satisfying This constraint is fulfilled at LO applying, under the condition / E max < m µ /2, the following integration limits: where m min 123 = m 2 µ − 2 / E max m µ + t 2 ; the other limits of integration are left unchanged. In the calculation of the NLO corrections to B( / E max ) we assumed the maximum missing energy / E max to be smaller than the photon detection threshold, and much lower than the muon mass. In figure 3b we show the branching ratio B NLO ( / E max ) versus the cut on the missing energy, in the upper panel, and its relative magnitude with respect to the LO prediction, in the lower panel. Error bands depicted in figure 3 are the assigned theoretical error due to renormalization scale variation. They are evaluated as in the case of the inclusive branching ratio. Errors due to numerical integration are typically smaller than the first.
In addition, we report in table 1 the branching ratios for different missing energy cuts: / E max = 1, 5, 10, 20, 50 and 100 m e ; our LO results are in good agreement with the values given in ref. [32]. 2 In ref. [32], where the LO decay is considered, a fit of the branching ratio in the endpoint region was presented: , with κ LO = 2.99 × 10 −19 . (4.8) We performed a similar fit employing as input data the NLO branching ratios for / E max = 1, 2, . . . , 10 m e . Taking into account the numerical error of B( / E max ), we obtained the following value for the constant κ at NLO accuracy: In (4.8) the exponent of / E max is fixed; relaxing such constraint and assuming it also to be a free parameter, i.e.

Discussion and Conclusions
In this work we studied the SM prediction of the differential rates and branching ratios of the muon decay with internal conversion µ → e(e + e − )νν at NLO. Virtual and real corrections were computed using the effective four-fermion Fermi Lagrangian plus QED and QCD, taking into account the full dependence on the mass ratio r = m e /m µ . We employed the library LoopTools for the numerical evaluation of the coefficients appearing in the Lorenz-covariant decomposition of tensor one-loop integrals. Real corrections were calculated with a phase space slicing method. For photon energies below the ω 0 threshold, the photon phase space integral is worked out analytically by taking advantage of the soft photon approximation for the bremsstrahlung amplitude. Above the IR cut-off, we employed the complete amplitude of the real photon emission process. A sufficiently small ω 0 is then chosen in the calculation of the physical observables to assure in the final results the cancellation of the ω 0 dependence between the soft and the hard part. Virtual and real corrections can be obtained as a Fortran code from the authors.
The branching ratio at NLO accuracy was presented in table 1. Our prediction is in agreement with the old measurement by the Sindrum experiment. In addition to this, we studied the decay rate as a function of the three-electron invariant mass, at the end-point region, and as a function of the cut on the missing energy.
From our results we can observe that while the global K-factor is close to unitythe shift δB NLO gives a correction of 2 × 10 −3 -locally, in the configuration where m µ − E 1 − E 2 − E 3 → 0, the relative size of radiative corrections is as large as 10 − 20% of the LO. Such an enhancement is given by the smallness of the / E max cut, which forces the bremsstrahlung photon to be emitted in the soft-collinear region, where the corrections can behave as (α/π) ln(m e /m µ ) ln( / E max /m µ ). Radiative corrections reduce the LO prediction of the width; the effect can be visualized in figure 3 shifting the curves downward. We estimated the theoretical uncertainty due to numerical errors in the Monte Carlo integration and renormalization scale variation.
Our results can be employed also for the evaluation of the tau decays τ → e (e + e − )νν and τ → µ (µ + µ − )νν by properly substituting m µ → m τ and m e → m , with = e, µ.

Acknowledgments
In the last phase of this work we learnt about an ongoing independent calculation of the NLO corrections to µ − → e − (e + e − ) ν µνe decay carried out by G. M. Pruna, A. Signer and Y. Ulrich. We would like to thank them for valuable discussion and for sharing the results of their work, which are in good agreement with ours in table 1.
Moreover, we would like to thank M. Passera for participating in the early stages of this work. We are also grateful to N. Berger for bringing the problem to our attention and for useful discussion and correspondence. The authors acknowledge the support from the Swiss National Science Foundation.

A Phase space decomposition
In this appendix we discuss the phase space parametrization employed in the Monte Carlo integration of the LO, virtual and real differential widths.

A.1 Phase space: LO and the virtual corrections
The generic n-body phase space element of a particle with momentum P decaying into n particles with momenta labeled by p 1 , . . . , p n is dΦ n (P ; p 1 , · · · , p n ) = δ 4 (P − it can be decomposed as nested sequence of two-body pseudo-decays applying recursively the splitting formula dΦ n (P ; p 1 , · · · , p n ) = (2π) 3 dq 2 dΦ j (q; p 1 , · · · , p j ) dΦ n−j+1 (P ; q, p j+1 , · · · , p n ), By means of eq. (A.2), we can express the five-body phase space shared by the LO and the virtual in the following way: dΦ 5 (P ; p 1 , p 2 , p 3 , t 1 , t 2 ) = dΦ 2 (t 12 ; t 1 , t 2 ) (2π) 3 dt 2 × dΦ 2 (p 12 ; p 1 , p 2 ) (2π) 3 dm 2 12 × dΦ 2 (p 123 ; p 12 , p 3 ) (2π) 3 dm 2 123 × dΦ 2 (P ; t 12 , p 123 ), where all momenta and invariant masses in (A.3) have been defined in section 2 (see also figure 1). The two-body phase space of a generic q → q 1 q 2 decay is, in the c.m.s., where q 2 = M 2 . The energy and the momentum of q 1 and q 2 in their c.m.s. is fixed by the masses and the invariant mass of the system: Employing the parametrization (A.4) in order to express each dΦ 2 in eq. (A.3) in its own c.m.s., and performing the trivial integration over the angles dΩ 123 and dφ * 3 , we get (A.7) Solid angles and momenta labeled with the superscript ' * ' are in the c.m.s. of the two electrons and the positron, where p * 123 = (m 123 , 0), while those with ' * * ' are in that one of the two electrons, where p * * 12 = (m 12 , 0). The allowed kinematic ranges in eq. (4.1) follow from energy conservation in each recursive splitting. Moreover, the condition on the missing energy (4.6) written in terms of the auxiliary momenta t 12 and p 123 , leads to the set of restricted integration limits in eq. (4.7).