Fully differential NLO predictions for the radiative decay of muons and taus

We present a general purpose Monte Carlo program for the calculation of the radiative muon decay $\mu\to e\,\nu \bar{\nu}\gamma$ and the radiative decays $\tau\to e\, \nu \bar{\nu}\gamma$ and $\tau\to\mu\,\nu \bar{\nu}\gamma$ at next-to-leading order in the Fermi theory. The full dependence on the lepton masses and polarization of the initial-sate lepton are kept. We study the branching ratios for these processes and show that fully-differential next-to-leading order corrections are important for addressing a tension between BaBar's recent measurement of the branching ratio $\mathcal{B}(\tau\to e\, \nu\bar{\nu}\gamma)$ and the Standard Model prediction. In addition, we study various distributions of the process $\mu\to e\,\nu \bar{\nu}\gamma$ and obtain precise predictions for the irreducible background to $\mu\to e \gamma$ searches, tailored to the geometry of the MEG detector.


Introduction
The study of muon decays has a very long and rich history and will continue to play a prominent role in particle physics [1]. It allows for precise measurements of Standard Model (SM) input parameters and is a powerful tool in searches for physics beyond the SM. In this work we are concerned with the radiative decay of the muon, µ → e ννγ. The branching ratio (BR) for this process depends on the energy cut on the photon ω 0 that is required to make the quantity well defined. For the standard choice of ω 0 = 10 MeV the BR is roughly 1%. Given the vast number of muons that can be produced, it should be possible to study radiative muon decays with very good precision.
Despite this, the current experimental knowledge on this process is rather limited. The standard value that is usually quoted for the BR is from the sixties and has an error of about 30% [2]. There are preliminary results of the Pibeta collaboration [3] with much better precision and MEG has performed a measurement [4] for rather stringent cuts on the electron and photon energies. Apart from measuring the BR, the SM can also be tested by measuring Michel parameters of a general formula for muon decays [3,5,6]. Furthermore, if the neutrinos have very little energy, this process is an irreducible background to searches for the lepton-flavour violating decay µ → eγ.
Related decays of the tau, τ → e ννγ and τ → µ ννγ have also been studied. The measurements of Cleo [7] for ω 0 = 10 MeV (in the rest frame of the tau) have recently been improved upon by the BaBar collaboration [8] and the BR are now known to about 3%. Furthermore, a preliminary analysis concerning the measurement of Michel parameters is also available from Belle [9]. Apart from the usual tests of the SM these decays also offer the possibility to gain information on the anomalous magnetic moment of the tau [10].
The radiative decay of leptons is most easily computed in the effective Fermi theory with a four-point interaction. Corrections beyond the Fermi theory due to the W -boson propagator [11,12] turn out to be much smaller than the NLO corrections. The tree-level calculation within the Fermi theory has been considered by several authors a long time ago [13][14][15][16][17]. Due to the photon bremsstrahlung the helicity of the final-state lepton does not have to be left-handed [18][19][20][21]. After some partial results [17,22] a full next-to-leading (NLO) calculation for the BR was presented in [23,24]. As for a related calculation of the rare decays of leptons [25], the results presented in [23,24] allow to obtain the differential decay width at NLO with cuts on the photon and electron energy and angles between them. In this article, we generalize these results and present a fully differential Monte Carlo program. This enables the implementation of arbitrary cuts, allowing to mirror the experimental situation more closely.
A fully differential calculation has a large impact for radiative decays of both, the muon and the tau. Starting with the latter, it has been noted [23] that the BaBar result for the BR of τ → e ννγ [8] has a tension of about 3.5 standard deviations with the NLO result. This measurement was done using rather stringent cuts on the final state and then converted to the standard cut of E γ ≥ ω 0 = 10 MeV. As we will show, the NLO corrections are important for this conversion and potentially resolve this tension. Regarding the muon, taking into account the geometry of the MEG detector allows to obtain a tailored background study for µ → eγ searches.
After briefly discussing some technical aspects of our calculation in Section 2 and giving some details about the Monte Carlo code in Section 3 we will study the BR for the rare decays of the muon and tau in Section 4. Whenever possible we compare with previous results in the literature [23,24] and find full agreement. Making use of the flexibility of our computation, we also study the impact of the cuts used in the measurement of the BR of the tau. In Section 5 we present several distributions related to the radiative decay of the muon, adapted to the studies of MEG and Mu3e. Finally, we present our conclusions in Section 6.

Methodology
This calculation of the radiative decay µ → e ννγ follows very closely the related calculation of the rare decay µ → e νν ee presented in [26]. The starting point is the QED Lagrangian with the Fierz rearranged effective Fermi interaction where P L = (1 − γ 5 )/2 is the usual left-handed projector. The QED part of the Lagrangian contains electron ψ e and muon fields ψ µ but no quark fields. Adding the tau field, ψ τ and making obvious modifications in the four-fermion operator the same Lagrangian can be used for the radiative decays of the tau. Despite this being an effective theory, it was shown by Berman and Sirlin [27] that all QED corrections are finite after the usual QED renormalization and that G F does not get renormalized. As noted previously, we keep m e = 0 and perform the computation for an arbitrary polarization of the initial-state lepton.
At leading order the photon in the radiative decay µ → e ννγ can be emitted from the muon or the electron so that two diagrams need to be calculated. At NLO there are eight one-loop diagrams, four mass-counterterm diagrams and six real emissions. The virtual matrix elements are generated using FeynArts [28] and FormCalc [29] and evaluated with FORM [30] and the scalar integrals were computed with the library LoopTools [31]. We have confirmed the matrix elements with a modified version of GoSam [32][33][34]. The modification for GoSam necessary to compute polarized amplitudes in (2.1) are detailed in [26].
Since we keep all lepton masses different from zero, the phase-space integration over the real matrix element generates soft, but no collinear singularities. The former are dealt with using FKS subtraction [35,36]. Since FKS treats soft and collinear singularities separately, its application in this case is particularly simple.
We renormalize the lepton masses and the electromagnetic coupling α in the on-shell scheme. The wave-function and the coupling renormalization factors, as well as the explicit result of the virtual and real corrections are regularization-scheme dependent. We have used conventional dimensional regularization (CDR) and the four-dimensional helicity scheme (FDH) [37] and verified that the final physical result is scheme independent, as described in [38].
We have checked that the corrections beyond the Fermi theory are suppressed by m 2 µ /m 2 W ≈ 10 −6 for µ → e ννγ . In particular, as for the normal muon decay, loop diagrams of the full Standard Model do not introduce logarithmically enhanced contributions log m 2 µ /m 2 W . Hence, for a NLO calculation it is sufficient to work within the Fermi theory.

Monte Carlo
The program presented in this paper is a standard parton-level Monte Carlo code, or more precisely a Monte Carlo integrator, written in Fortran 90. Weighted events with tree-level kinematics (for the Born term and the virtual corrections) or with an additional photon (for the real corrections) are generated in the rest frame of the incoming muon or tau. These events can then be used to define arbitrary observables and implement arbitrary cuts. The numerical phase-space integration is done using VEGAS [39]. The code with instruction how to use it is available from the authors upon request.
In the case of the real corrections, following the FKS method, we also have to compute the soft counter event. The FKS procedure introduces an unphysical parameter, usually denoted by ξ cut [35]. This parameter determines the size of the soft region for which the soft subtraction is carried out. The dependence of the real corrections on ξ cut is cancelled by the corresponding dependence of the integrated counterterms. We have checked our implementation by verifying the independence on ξ cut of the full final result, virtual plus real corrections. Running a particular calculation with different values of ξ cut serves as an indirect check on the infrared safety of the observables.
Since we keep m e = 0 there are no collinear singularities. However, to ensure numerical stability we have to control pseudo singularities that arise due to the smallness of m e . To do so we choose a custom tailored phase-space generator over a recursive generation to align the pseudo singularities with the variables of the VEGAS integration. This optimizes the adaption of the grid for the integration and leads to very stable numerical results.
Running times obviously depend on the required numerical precision. The very precise results for the BR given in Table 1 require about 10 CPU hours per process (on a reasonably modern machine), whereas a distribution with very precise numbers for all bins, like the one given in Figure 5, takes of the order of 1000 CPU hours. Obviously, an arbitrary number of distributions can be computed in one go and very reasonable results can be obtained reducing the running time by a factor 5 − 10. Typically, about 50% of the time is spent on the virtual corrections and the (soft-subtracted) real corrections each, whereas the integrated counterterm and the Born term need about 5% and 1% of the total time, respectively.
The amount of Monte Carlo statistics necessary to obtain small numerical uncertainties increases for sharply falling tails of distributions. To improve the predictions in such cases we configured special runs that focus on this region. This has been done for the tail in panel (a) of Figure 3 where we dedicated a special run for the region E / ≤ 6 MeV. Furthermore, for the two-dimensional distribution presented in Figure 4 the various rows have also been computed in separate runs.

Branching ratios
In this section we present some results for the branching ratios with (l, ) ∈ {(τ, µ), (τ, e), (µ, e)} and Γ l = 1/τ l is the experimentally measured width of the lepton. To make the branching ratio well defined and avoid a soft divergence, a minimum photon energy ω 0 has to be imposed. The usual choice is ω 0 = 10 MeV. We have used the following standard values for the various inputs [40]:  Apart from the mean life of the tau they are all known to a precision that is well beyond what we need. Table 1 summarizes our results for the branching ratios for the three radiative lepton decays at leading order, B LO , and next-to-leading order. In the latter case we distinguish between the exclusive and inclusive case. For B excl (B incl ) we require exactly (at least) one photon satisfying the cuts (in this case energy E γ > ω 0 ).
Let us start by focussing on the first three columns in the top half of Table 1. There we give the branching ratios with ω 0 = 10 MeV and no other cuts applied. Our results agree well with an earlier calculation presented in [23]. The error given is only the numerical error from the Monte Carlo integration.
The size of the NLO corrections for B excl are 10%, 2.5% and 1.7% for (τ, e), (τ, µ) and (µ, e), respectively. This is in agreement with the parametric expectation for the NLO corrections, (α/π) ln(m l /m ) ln(ω 0 /m l ) B LO , which amounts to roughly 10% for (τ, e) and 3% for (τ, µ) and (µ, e). Accordingly, we expect the theoretical error due to omitted NNLO corrections to be about 1% for (τ, e) and considerably smaller for (τ, µ) and (µ, e) [23]. As we will see, these estimates are only valid for a sufficiently inclusive branching ratio. In the case of the tau, there is an additional error of about 0.2% due to input value of the tau mean life.
In order to compare with the recent measurement of MEG [4] we also present in the last column of Table 1 the µ → e ννγ branching ratio for the MEG cuts, i.e. ω 0 = 40 MeV and E e > 45 MeV. Once more, our results agree with an earlier calculation [24]. These cuts are very restrictive and reduce the branching ratio by more than five orders of magnitude. Furthermore, through these cuts the NLO corrections are increased from 1.7% to 5.6%. Of course, energy conservation forbids the production of two photons with E γ > 40 MeV together with E e > 45 MeV. Hence, there is no distinction between B excl and B incl in this case.
For the radiative decay of the muon, the results presented in Table 1 agree with the experimental values B exp [2,4], albeit with a large experimental error. The situation for the tau decays is less clear. The most precise measurements of these branching ratios have been made by BaBar [8,41]. As argued in [23] these measurements are to be compared with B excl . For the (τ, µ) decay the agreement is satisfactory but for the (τ, e) decay there is a 3.5 σ discrepancy [23]. In the remainder of this section we will revisit this issue, making use of our fully differential NLO computation to match the actual measurement as closely as possible.
Let us start with the most problematic case, the decay τ → e ννγ. For the BaBar measurement, tau pairs are produced through e + e − collisions at √ s = M Υ(4S) = 10.58 GeV. The event is then divided into a signal-and tag-hemisphere. In order to reduce background events, rather stringent cuts on the kinematics of the decay products e and γ in the signal hemisphere are applied. In particular, the following requirements are made: All the quantities are given in the centre-of-mass frame. These cuts can be easily implemented in our code. To this end, we generate taus in their rest frame, boost them to a frame such that they have energy √ s/2 and then apply the cuts (4.3) in this boosted frame. As we will see, the NLO corrections will have an important effect when 'undoing' the cuts, i.e. when extracting B excl (with only the cut E γ ≥ 10 MeV in the tau rest frame).
In order to illustrate this we have devised the following simplified scheme: Let N obs be the measured number of events including all cuts. To obtain the branching ratio this is multiplied by a factor exp again refer to the cut E γ ≥ 10 MeV and the cuts (4.3), respectively. More precisely, we require that exactly one photon passes the cuts. To assess the importance of NLO corrections when extracting B exp we write Thus, is a purely theoretical factor that describes the difference of using a LO or NLO computation in the determination of B exp . The results for the various factors described above are given in the first row of to · B exp = 1.704(50) · 10 −2 , in much better agreement with the theoretical NLO result B excl = 1.645(1) · 10 −2 . We recall that we expect an additional theoretical error of the order of 1% due to the uncomputed higher-order corrections.
Of course, the same procedure can be repeated for the τ → µ ννγ decay. In this case, some of the cuts applied by BaBar are A computation of the factor reveals that the effects here are more modest and amount only to a correction of about 1%. The resulting value · B exp = 3.65(10) · 10 −3 agrees well with the NLO result B excl = 3.571(1) · 10 −3 . For completeness, we also apply a similar procedure to the MEG result for the branching ratio of µ → e ννγ. The cuts given in (5.1) are applied in the rest frame of the muon. As can be seen from the result in the last column of Table 1, in this case the NLO effects are marginal and it is sufficient to use a LO description to remove the effect of the angular cuts in (5.1).
Obviously, this is only a simplistic and by far not complete simulation of the full analysis. While the cut on E * γ has the biggest impact, the results for the factor actually depend quite significantly on all the details of the cuts. In particular, in the presence of a second photon it is important to precisely specify how the cuts are applied. This can also be seen from the rather large difference between B excl and B incl for τ → e ννγ. We do not claim that this is the conclusive resolution to the apparent 3.5 σ deviation for the measured branching ratio of τ → e ννγ. However, we do claim that a proper inclusion of NLO effects is mandatory for such a measurement, in particular if stringent cuts on the decay products are applied.
Finally, we consider the branching ratio for µ → e ννγ with ω 0 = 10 MeV and an additional cut on the angle between the electron and the photon, θ eγ > 30 • . For this BR, a preliminary result B exp = 4.365(43) · 10 −3 is available [3] from the Pibeta collaboration. The error in (4.8) is the numerical error from the Monte Carlo integration only. Following the arguments given above, a conservative estimate of the theory error due to higherorder correction is 0.2%. There seems to be a tension between theory and experiment. However, the experimental result is only preliminary and the SM value quoted by the collaboration disagrees significantly with our result. Hence, further clarification is required before a conclusion can be drawn.

Distributions
While branching ratios with only an ω 0 cut are useful standard quantities, differential distributions with arbitrary cuts offer a more direct comparison between experiment and theory. Furthermore, they are important for precise background estimation. As we have seen in the previous section, not taking into account the cuts systematically at NLO can have surprisingly large effects. In this section we will provide some examples of such distributions that can be created with our Monte Carlo program. We will restrict ourselves to the radiative decay of the muons, but obviously similar results can be obtained for radiative decays of the taus.
A very promising place to study the process µ → e ννγ is the MEG experiment [42]. Because it was designed to search for LFV muon decays its detector geometry creates a very restrictive phase space. We mimic the MEG detector by imposing the following cuts: | cos (p e , z)| ≡ | cos θ e | < 0.5 , |φ e | < π 3 . (5.1c) The muon polarization is set to be 85% along the z-axis such that P µ = −0.85 z.
In addition, we require one and only one photon, meaning that an event with additional real radiation with energy larger than the detector resolution of roughly 2 MeV will be discarded if it hits the detector, i.e. also follows (5.1b).
This approximation is justified even though it would reject a pair of collinear photons because, in contrast to QCD, there is no mechanism preferring collinear photons. The cuts enforced by the finite size of the electron detector (5.1c) would automatically be satisfied indirectly through the photon cuts (5.1b), if they were emitted back-to-back.  In order to analyse the loss of events that is due to this approximation we consider the distributions dB d(cos θ e ) and dB dφ e , (5.2) without implementing the cuts (5.1c). These distributions of the branching ratio B are shown in Figure 2. Integrating them at NLO in the interval specified by (5.1c) we find 1 B Hence the θ e and φ e cuts result in a loss of about 4% and 7%, respectively. Due to the non-vanishing polarization, the cos θ e distribution is not symmetric w.r.t. zero. We also see that the corrections are negative and amount to roughly 5-10%. As we will see, this is a generic feature for the distributions, except for special corners of phase space. For illustrative purposes, we now consider some distributions with the full cuts, i.e. including (5.1c). We start with the missing energy distribution, dB/dE /, as it is of special interest to the MEG experiment. In the region with very little invisible energy the radiative decay of the muon is an irreducible background to the lepton-flavour violating decay µ → eγ that is searched for by MEG. We define the invisible energy as The result is shown in the top-left panel of Figure 3. For the bulk of the distribution, the corrections are of the order of −5%, but in the tail they increase substantially. We also note that the distribution itself falls rapidly towards zero for E / → 20 MeV, due to the kinematic constraints.
In the top-right panel of Figure 3 we show the distribution of the cosine of the angle between electron and photon, cos θ eγ . Due to the energy cuts, the neutrino energy is  restricted to be small and the electron and the photon are nearly back-to-back. Indeed, more than 99% of the events lie between −1 ≤ cos θ γe ≤ −0.9. Again, the corrections are −5% for the bulk of events. Because the energy of the electron and the photon can be measured by the MEG experiment they, too, are of particular interest. It is customary to show the energy fractions instead of the energy itself. Due to the cuts (5.1a) we have x 0.85 and y 0.76. Our NLO predictions for the spectrum dB/dx and dB/dy are shown in the lower panels of Figure 3.
Additionally to single differential distribution we have also computed the double differential distribution d 2 B/dE γ dE e ∼ d 2 B/dxdy to demonstrate the impact of the energy cuts. For this plot only the geometric MEG cuts (5.1b) and (5.1c) were used. Note that even though E γ starts at 0 in Figure 4, this is unphysical as a cut on the photon energy is required to ensure infrared safety. We take E γ ≥ 2 MeV. To be precise, apart from an elec-  Figure 4. The double differential distribution d 2 B/dE γ dE e . The K factor is shown through the colouring. A hard cut of E γ ≥ 2 MeV was imposed on the visible photon. tron with (5.1c) we require exactly one photon with E γ ≥ 2 MeV and (5.1b). Additional photons must satisfy (5.1d). The distribution has the expected sharp fall for increasing x ∼ E e and y ∼ E γ . The corrections are particularly important in the region of large x, i.e. for large electron energies. This is also expected from the bottom-left panel of Figure 3.
As a final illustration of the flexibility of our code, we now turn to the Mu3e experiment [43,44]. In the normal running mode Mu3e cannot detect photons. However, it is in principle possible to detect them should they convert into electron positron pairs of sufficient energy. In order to mimic this situation we consider the following cuts: The cuts on cos θ are a simplistic way to include the geometry of the detector. We require one and only one photon in the detector, meaning that an event with real radiation of more than the threshold of 20 MeV will be discarded if | cos θ γ | < 0.8. As for the MEG results, we set the polarization P µ = −0.85 z. As a demonstration we show the photon energy distribution dB/dy in Figure 5. We refrain from showing the LO curve, as it is basically indistinguishable from the NLO curve. In accordance with the fact that the cuts (5.5) are not very restrictive, the size of the corrections is rather modest. However, the K-factor, shown in green, is far from universal. The kink in the K factor at y 0.8 is due to the cut on E e . This shows once more that cuts can have a large impact on the shape and size of NLO corrections.

Conclusion
In this article we have presented a general purpose parton-level Monte Carlo program for the NLO corrections within the Fermi theory of the radiative decay of leptons. This generalizes an earlier NLO computation [23] of the branching ratio for these decays in that now arbitrary cuts can be implemented. By today's standards this is not a very complicated process and it is somewhat surprising that such a Monte Carlo program had not been presented earlier. Possible reasons might be that until recently, the experimental result were not very precise and the (pure QED) corrections could have been expected to be very modest, i.e. at the order of 1%. However, in many important circumstances, the corrections are much larger, as pointed out also in [23,24]. In regions of phase space corresponding to relatively stringent cuts, the corrections can easily reach 10%. This is the case for background studies to µ → eγ searches as well as for recent measurements of the branching ratio of radiative tau decays.
With increasing precision of the measurements [3,8] a fully-differential treatment at NLO is mandatory to obtain reliable predictions. Indeed, we have shown that it seems very plausible that the 3.5 standard deviation discrepancy between the BaBar measurement of B(τ → e ννγ) and the NLO SM result is related to not using a full NLO calculation when estimating the efficiency. With the program presented here, there is no longer any reason not to use a full NLO calculation and we hope this helps in making more precise comparisons between theory and experiment for radiative lepton decays in the future.