Fully differential NLO predictions for the rare muon decay

Using the automation program GoSam, fully differential NLO corrections were obtained for the rare decay of the muon $\mu\to e\nu\bar\nu ee$. This process is an important Standard Model background to searches of the Mu3e collaboration for lepton-flavour violation, as it becomes indistinguishable from the signal $\mu\to 3e$ if the neutrinos carry little energy. With our NLO program we are able to compute the branching ratio as well as custom-tailored observables for the experiment. With minor modifications, related decays of the tau can also be computed.


Introduction
The rare decay of the muon µ → eννee is not one of those Standard Model processes that has been in the limelight in the past decades. The branching ratio B = (3.4 ± 0.4) × 10 −5 has been measured more than thirty years ago by the SINDRUM collaboration [1,2] in the context of searches for the lepton-flavour violating process µ → 3e. For the corresponding decay τ → eννee a measurement B = (2.8 ± 1.5) × 10 −5 is available from CLEO [3], but for other leptonic five-body decays of the τ only upper limits exist.
On the theory side, several tree-level calculations of the branching ratios for the various rare decays have been made long ago [4][5][6]. These calculations might have been sufficient for the previous experimental situation. However, with the upcoming Mu3e [7][8][9] experiment that is dedicated to the search for the lepton-flavour violating decay µ → 3e at the Paul Scherrer Institute (PSI) an improved theoretical description of the rare muon decay is highly desirable for several reasons. First, the rare muon decay is an interesting process in its own right as it can be measured very precisely by Mu3e. Second, in the limit where the neutrinos have very little energy, it is a background to the decay µ → 3e.
In order to avoid infrared singularities and to get a reliable prediction in all corners of phase space, it is necessary to keep finite electron mass terms. Furthermore, a realistic background study benefits from a fully differential description of the rare muon decay. This requires the computation of the corresponding matrix elements and their implementation in a Monte Carlo program.
The squared matrix elements at tree level for the unpolarized rare muon decay and a study of the energy spectrum have been published less than ten years ago [10]. For a more realistic investigation of the background to µ → 3e, the polarization of the muon has to be taken into account. Even though the tree-level calculation of the matrix elements in the Fermi theory is trivial, these matrix elements have been made available to the Mu3e collaboration only very recently [11]. A tree-level study of the rare muon and tau decay has been published at a later stage [12].
Typically the next-to-leading (NLO) order QED corrections are expected to be at the percent level. However, in extreme regions of phase space, as for example at the endpoint of the lepton energy spectrum, corrections can be enhanced due to large logarithms. As seen in a recent NLO calculation [13] of the radiative muon and tau decays → ννγ, QED corrections can easily be as large as 10%. Thus, a fully differential NLO Monte Carlo for the rare decays of a polarized muon is an important ingredient to fully exploit the future data to be taken by the Mu3e collaboration. In this article we present such a calculation. We use GoSam [14] to obtain the matrix elements and perform the phase-space integration using FKS subtraction [15,16].
The paper is organized as follows: in section 2 we will discuss our methodology. Since this is a standard NLO calculation we focus on changes to GoSam that were necessary to compute the required one-loop matrix elements. In section 3 we present some branching ratios that were independently confirmed by [17]. Finally, in section 4 we present some distributions and conclude in section 5.

Methodology
We perform the computation in the Fierz-rearranged effective 4-Fermi interaction, i.e. with the Lagrangian where P L = (1 − γ 5 )/2 is the usual left-handed projector. Even though this is an effective theory, it can be shown that the Fermi constant G F does not get renormalized and all QED corrections are finite after the usual QED renormalization [18]. The QED Lagrangian L QED includes electron ψ e and muon fields ψ µ (and tau fields for the rare tau decay) but no quark fields. The effect of the latter is very small for muon decays. As mentioned in the introduction we keep m e = 0.

The rare decay
At leading order the rare decay µ → eννee is given through four diagrams. At NLO there are about 40 one-loop diagrams and 20 diagrams involving a real emission. These matrix elements are generated using the automation tool GoSam [14] and reduced at run time using Ninja [19][20][21] or golem95 [22,23]. Scalar integrals were computed using OneLOop [24,25]. The arising soft singularities from the real emission diagrams are treated using FKS subtraction [15,16]. In the absence of collinear singularities, the FKS method is particularly simple as it generally treats soft and collinear singularities separately.
Finally, a numerical integration of the full phase space is carried out using VEGAS [26]. This allows for the production of any differential observables with arbitrary cuts. Instead of a general purpose phase-space generator that creates momenta recursively, we use a tailored phase-space generator. It is designed such that pseudo singularities, which cause numerical instabilities due to the smallness of m e , are aligned with the variables of the VEGAS integrator to optimally utilize the VEGAS adaption.
For the renormalization of the lepton masses we always choose the on-shell scheme. The QED coupling can be renormalized either in the on-shell scheme α = α os or in the MS schemeᾱ = α MS (µ = m µ ). The results can of course easily be converted usinḡ Since GoSam returns the NLO matrix elements in the four-dimensional helicity scheme (FDH) [27], care has to be taken to use this scheme also for the external wave-function renormalization and the real corrections [28].
Regarding the γ 5 that is present in (2.1) we note that the radiative corrections to the axial vector contribution are related to those of the vector contribution by setting m e → −m e [18]. This is due to the fact that the simultaneous transformation ψ e → γ 5 ψ e and m e → −m e leaves L QED invariant but exchanges the vector and axial-vector currents. Using an anticommuting γ 5 , as done in GoSam [29], respects this relation and hence avoids problems with γ 5 in this particular case.

Changes in GoSam
The GoSam package is designed to compute one-loop amplitudes for processes involving an arbitrary number of particles in QED and QCD. To use the 4-Fermi interaction (2.1), the corresponding Feynman rules have to be incorporated into GoSam's model file.
GoSam applies the spinor helicity formalism using light-cone decomposition. For a lepton of mass m the momentum q ρ is decomposed with the help of two light-like vectors ρ and n ρ as Then the spinor for the massive lepton is written in terms of the usual massless spinors as Unfortunately, the reference vector n ρ relative to which the polarization is defined is always chosen by GoSam to be an internal vector, typically the momentum of a massless particle in the process. While this greatly simplifies the amplitudes that are evaluated, it is inconvenient to compute polarized observables this way. An external vector as the reference vector is needed. This can be changed by adapting the GoSam process template such that the light-cone decomposition of massive fermions is done with respect to the new external light-like reference instead of an internal light-like vector. Additionally, the actual GoSam code needs to be adapted to allow for massive lepton counterterms.

Branching ratios
In this section we present some results for the branching ratio where Γ µ = 1/τ µ with τ µ = 2.19698 × 10 −6 s is the experimentally measured width of the muon. We have used the standard values for the various inputs [2]: In addition to the full branching ratio we also consider the branching ratio with a constraint on the invisible energy where E 1 ≥ E 2 are the energies of the positrons and E 3 is the energy of the electron, respectively.
In table 1 we summarize the branching ratios in the on-shell scheme for various cuts on E /. The error indicated in the table is the numerical error of the Monte Carlo integration only. For ever more stringent cuts on E / it becomes increasingly more difficult to obtain stable numerical results with our general purpose Monte Carlo program. As can be seen, the corrections are moderate if no cut is applied, but become substantial for stringent cuts on E /. Results for the branching ratio with such stringent cuts on E / are particularly relevant for a background estimate of the Mu3e experiment. Furthermore, these results allow for a comparison with [10] at leading order and with [17] at NLO. We find perfect agreement in all cases, assuming that [10] used Γ to normalize the branching ratio. Of course, we can also compute the leptonic five-body branching ratios of the τ . As an example we give
value for Γ τ = 1/τ τ with τ τ = 2.903×10 −13 s and m τ = 1776.8 MeV. As for the muon case, in the on-shell scheme the corrections are very small in the absence of additional cuts. In the MS-scheme, the NLO corrections are about −3.5%. We note that the muonic and hadronic photon vacuum polarization contributions both give a very small effect O(10 −9 ) to B τ e . The latter have been estimated using the Fortran code hadr5n12 [30][31][32]. The result given in (3.4) is to be compared with the experimental result from CLEO B τ e = (2.8±1.5)×10 −5 [3]. For other leptonic decays of the τ only upper limits are available from experiment. They can be computed easily with our code even if processes with two different lepton flavours in the final state require trivial modifications.

Distributions
In addition to branching ratios, differential distributions with arbitrary cuts can also be computed with our Monte Carlo program. To provide some examples, the following distributions will be considered: • dB/dE /: The invisible energy distribution is important to correctly estimate the background from the rare decay to the LFV decay µ + → e + e + e − .
• dB/dx i where x i ≡ 2E i /m µ : The momentum fraction distributions of the three charged final-state leptons can be used to experimentally discriminate between different LFV models.
• dB/d cos θ i where θ i ≡ (q i ,ẑ) is the angle between the z-axis and the momentum q i of the various final-state leptons: The angular distributions of the charged leptons can be used to study the V − A structure of the Fermi interaction, or in the case of LFV, the Lorentz structure of the effective operator.
For illustration, these calculations were carried out using the polarization s = −0.85ẑ and imposing E i > 10 MeV and | cos θ i | < 0.8 .

Invisible energy spectrum
One of the most important distributions, dB/dE /, is shown in figure 2 with numerical errors indicated by the error bars. This distribution falls sharply in the region E / → 0 that is particularly important for the Mu3e experiment. Thus, using a standard configuration (shown in blue in figure 2) it is challenging to get enough statistics in the low energy tail to obtain predictions with reasonable statistical uncertainties. To improve the tail, we focused on the low energy region by imposing an additional cut on E / of 20 MeV (shown in orange), 10 MeV (shown in green) and 5 MeV (shown in red) and combined the results. As can be seen from figure 2, the NLO corrections are negative except for a small region of maximal E /. In the low-energy tail, the corrections exceed −10%. Hence, there are substantially fewer background events to µ → 3e from the rare decay than expected from tree-level simulations. The cuts E i > 10 MeV are the reason for the sharp fall of the distribution at E / → m µ − 30 MeV. The kink in the distribution is at about m µ /2, shifted to somewhat lower values due to the effects of the non-vanishing electron mass. In fact, due to the additional real radiation of a photon, the NLO corrections amount to shifting the distribution dB/dE / to higher energies.

Momentum fraction and angular distribution
As an example of a distribution where the polarization of the muon has an effect, we consider the angular distribution of the charged leptons. We consider three cases: no cut on the missing energy, a cut E / ≤ 20 MeV and finally a cut E / ≤ 10 MeV, focusing on events in the tail of the E / distribution. In the left panels of figure 3 we show the normalized distributions dB/d cos θ at NLO for the more energetic positron (blue), the less energetic positron (orange) and the electron (green). In the right panel of figure 3 we show the normalized energy fraction distributions. The K-factor for the positrons is shown in the subpanels. The electron's K-factor is virtually indistinguishable to the one of the soft positron.
As observed already in section 3, the size of the corrections tends to increase if cuts are applied. With the cuts of (4.1) the corrections for the full branching ratio amount to −1.7% (−3.4%) in the on-shell scheme (MS scheme). If in addition we require E / ≤ 20 MeV the corrections increase to −7.3% (−9.1%) and for E / ≤ 10 MeV they are −10.2% (−11.9%). At NLO, the difference between results in the on-shell and MS scheme is less than 0.5%.
With the chosen polarization of the muon, if no cut on E / is applied, all final state leptons prefer to be emitted backwards, as is the case for the normal muon decay. However, with a cut on E / the hard positron behaves distinctly different to the other two leptons. Loosely speaking, this can be understood by noting that the more energetic positron is the primary positron, whereas the soft positron (and the electron) are produced from the conversion of the internal photon. As the cut on E / becomes more restrictive, this effect becomes more pronounced. The soft positron behaves ever more similar to the electron. This can also be seen in the energy distribution where the hard positron behaves drastically different from the soft positron and the electron, as shown in the right panel of figure 3. Furthermore, the size of the NLO corrections increases with more restrictive cuts, in analogy to what can be seen in table 1 and figure 2.
We stress that the distributions shown in figures 2 and 3 just serve as an illustration and that we can compute any distribution with arbitrary cuts at NLO. In regions of phase space, where there are not many events, the numerical error is large and a dedicated run is needed to improve the numerical precision.

Conclusion
We have presented a Monte Carlo code for the rare polarized muon decay µ → eννee. This code allows to compute branching ratios and distributions with arbitrary cuts at NLO in the Fermi theory. Our results for the branching ratio have been compared with Ref. [17] and full agreement has been found.
If no stringent cuts are applied to the final state, the size of the corrections in the on-shell scheme is very modest, i.e. at the order of 1%. Whilst it is notoriously difficult to give a reliable theoretical error of higher-order calculations, this is a good indication that the NLO calculation in this case provides a theoretical prediction with an error well below 1%. Corrections beyond NLO in the Fermi theory are very unlikely to be larger. Effects beyond the Fermi theory are suppressed by the large electroweak mass scale M ew and are of the order (m µ /M ew ) 2 ∼ 10 −6 .
In connection with background studies for the Mu3e experiment, it is often important to consider stringent cuts. In particular, the region of phase space where the invisible energy E / is very small is of importance. With such cuts the corrections can easily reach 10% or even more. The fact that the corrections are large and negative is favourable to µ → 3e searches. However, the size of the corrections indicates that in this case corrections beyond NLO in the Fermi theory are likely in the region of 1%.
From a phenomenological point of view five-body leptonic decays of the muon are much more important than the corresponding decays of the tau. Nevertheless, the latter can also be obtained at NLO with minor modifications of our code.