Three-loop corrections to the muon and heavy quark decay rates

Recently, $O(\alpha^3)$ corrections to the muon decay rate and $O(\alpha_s^3)$ to heavy quark decays have been determined by Fael, Sch\"onwald, and Steinhauser. This is the first such perturbative improvement of these important quantities in more than two decades. We reveal and explain a symmetry pattern in these new corrections, and confirm the most technically difficult parts of their evaluation.


I. INTRODUCTION
Muon decay provides one of the pillars of the Standard Model, the Fermi constant G F [1]. Its determination relies on precise measurements of the muon lifetime [2,3] and mass [4], and on theoretical evaluation of radiative corrections. These corrections arise primarily from quantum electrodynamic (QED) interactions involving the muon and the daughter electron, and are expressed as a power series in the fine structure constant α ≃ 1/137. Very recently, a new term in this series has been calculated [5] using an expansion in the difference of the muon and electron masses. Here we confirm the most demanding parts of the expansion published in [5]. We also explain a pattern governing the first two terms of the expansion.
Determination of the coefficients of increasing powers of α requires monumental efforts. A new result appears very rarely, approximately once every one or two generations of theorists. Each new result is witness to a new technology available in perturbative quantum field theory. First order corrections, calculated in 1955 [6], were in fact the first loop effects calculated for a decay process. They served as a model for subsequent studies of quantum chromodynamics (QCD) processes in heavy quark decays [7,8].
It took more than 40 years before the second-order coefficient became known [9,10]. That progress was made possible by the technique of recurrence relations based on integration-by-parts [11,12] and on symbolic manipulations with computers [13]. Now, another two decades later, the coefficient of α 3 has been determined [5]. Several theoretical developments underly this advancement. Laporta algorithm [14] allows to express multi-loop integrals in terms of a relatively small number of master integrals. This algorithm can be implemented in powerful symbolic algebra software capable of parallelization [15,16]. Finally, despite the electron being about 207 times lighter that the muon, the calculation is done as an expansion around the situation where the electron and the muon masses, m e and m µ , are equal [17][18][19].
Historically, radiative corrections to the muon decay were calculated before corrections to heavy quark decays, because QCD was developed only later. Computational technology is the same in both QED and perturbative QCD, and once the QED corrections are known, it is relatively easy to evaluate additional diagrams involving multi-gluon vertices. Ref. [5] presents corrections both for the muon and for the b-quark decays. Wherever our discussion below is relevant for both types of processes, we refer to the decaying particle's mass as M , and denote the mass of the produced charged particle by m. We use δ = 1 − m/M to denote the expansion parameter: the relative difference of masses.
In Section II we present our calculation of five terms of the expansion in δ for a subset of corrections, including the most demanding part with all quanta (photons or gluons) coupling to incoming or outgoing fermions, rather than, for example, to virtual loops. In Section III we demonstrate how the first two terms in the expansion in δ can be found from a form factor, first determined in 2004 in Ref. [17]. Section IV contains our conclusions.

II. EXPANSION IN THE MASS DIFFERENCE
As we shall demonstrate in Section III, the first two terms of the expansion of the decay rate in the mass difference may be obtained without a new calculation. Unfortunately, the remainder of the expansion requires a challenging calculation. For this publication, we have evaluated the first five terms of the expansion, δ 5,...,9 of the three most difficult contributions: those proportional to C 3 F , C 2 F n b , and C F n 2 b where n b labels loops containing a fermion of mass M . (We use the standard notation for the SU(N ) factors: C F = (N 2 − 1)/(2N ), C A = N , T F = 1/2 with N = 3 for QCD and C F = 1, C A = 0, T F = 1 for QED [20].) Our method, originally developed for corrections of O α 2 or O α 2 s in Ref. [19], is the same as the one employed in Ref. [5] for O α 3 s (α s denotes the strong coupling constant). However, we used a different software implementation.
In this section we briefly describe our method using the example of the muon decay. It was found in Ref. [19] that an expansion in terms of δ is useful for computing high order corrections. In particular, the calculation is much less involved than an expansion around m/M = 0, and converges very well all the way to the physical value of m c /m b or m e /m µ , even though m e /m µ = 0.005 is very far from the equal mass limit.
Our crucial tool is the optical theorem. It allows us to write contributions to the decay rate of the b-quark or muon as self energy diagrams. To generate the diagrams, we use DiaGen/IdSolver [21] modified to work with the types of propagators and expansions that appear in the present calculation. This code also produces the required integration by parts (IBP) identities and carries out the reduction to master integrals. The reduction of the large number of necessary integrals required the use of the WestGrid cluster which is part of the high performance computing infrastructure Compute Canada.
The decay is induced by an effective Fermi interaction with the W boson contracted to a point, as shown in Fig. 1. Corrections due to the W propagator including electroweak effects are known up to two loops [22]. Here we are interested in three loop QED corrections in the limit of a heavy W . They are described by five loop self energy diagrams (the extra two loops are responsible for the tree-level decay in Fig. 1(b)). We show examples of such diagrams in Fig. 2. The integral over the neutrino loop is carried out first by integrating over k 1 shown in Fig. 1 Examples of third-order radiative corrections to the muon decay. Wavy lines represent photons, contributing to both virtual and real radiation. In our approach, the imaginary part of such five loop diagrams is calculated.
be done at any QED order because the neutrinos do not interact with each other in the approximation considered here. The number of loops is thus reduced by one, at the price of introducing a propagator with non-integer power. This complicates the IBP reduction. For the remaining four loop integrals, we use an asymptotic expansion in δ, which we now describe.

A. Asymptotic Expansion
General properties of the asymptotic expansion of the loop integrals can be explained with the help of Fig. 1. Since we have already integrated over the neutrino loop, we only need to consider k 2 when carrying out the expansion. The following integral remains where we use dimensional regularization with d = 4 − 2ǫ. Without QED corrections, there are two regions: soft with components k 2 of the order of M − m, thus much smaller than M ; and hard with k 2 ∼ M . In the hard region, the term M 2 δ(2 − δ) can be treated as small. After expanding in it, the integral takes the form This is an on-shell self-energy integral with no imaginary part and so it does not contribute to the decay. This is the crucial observation which holds in any order in QED: if the momentum flowing through the neutrinos is large, the resulting integral has no imaginary part. This reduces the number of regions that must be considered by half.
At O α 3 this means that we only need to consider eight regions instead of the general 16. Moreover, among the non-contributing regions is the computationally very difficult one with all five hard loops.
In the soft region, k 2 ∼ δ · M and we simplify the integral in Eq. (1) by noticing that some terms are suppressed by an extra power of δ, Thus we encounter integrals with so called eikonal propagators [23,24] whose general form is where ∆ does not depend on the integration variable k. We have included another eikonal propagator that appears in higher-order loop diagrams. As it turns out, in all cases ∆ takes the form of a constant or another eikonal propagator. This makes it possible to apply this integral iteratively when computing the required master integrals.

B. Master Integrals
Expansion in relevant momentum regions of all diagrams leads to a large number of eikonal and on-shell integrals. We encounter at most three-loop on-shell integrals because the four-loop on-shell integrals contribute only to the real part and not to the imaginary part that we need.
Fortunately the three-loop on-shell master integrals have been computed in Ref. [25] to the required order in ǫ for the present calculation. We note that these integrals were originally computed in Ref. [26] in the context of the fermion anomalous magnetic moment, (g − 2), but it was not until [25] that enough terms in the ǫ expansion where known to make the present calculation possible.
The soft integrals appear up to four loops. After reducing to master integrals, there are one one-loop, one two-loop, three three-loop and 13 four-loop master integrals that must be computed. The general one-loop eikonal integral in Eq. (4) is computed in [23]. The result has the property As discussed in Section II A, the form of ∆ in this result makes it possible to iteratively apply the one loop computation to all two and three-loop and almost all four-loop eikonal type master integrals that are required. The remaining four-loop eikonal integrals can be computed analytically using a single Mellin-Barnes parameter. In all cases, the Mellin-Barnes integral can be evaluated by summing over the poles of the Gamma functions. Thus all eikonal integrals can be written explicitly in terms of Gamma and hypergeometric functions which are then expanded in ǫ to the needed order. For the expansion of hypergeometric functions, we use the Mathematica package HypExp [27]. Detailed examples of evaluation of eikonal integrals can be found in Ref. [28,29]. The final ingredient are the three-loop renormalization constants. The mass dependent constants are known in an expansion around the limit m/M = 0 [30]. We require the opposite expansion and thus had to compute the wave-function and mass renormalization terms. The calculation was done in an identical way to [30] except, of course, expanding around m/M = 1.
Interestingly, authors of [30] were able to present their solution only in terms of one-and two-dimensional integrals which must be computed numerically. In the limit m/M = 1 however, we obtain fully analytic results. [After this calculation was completed, we learned about Ref. [31] where analytic results in the limits m/M = 0, 1, ∞ are presented.] Using the approach described in this section, we have been able to confirm terms in the first five orders, δ 5 , . . . , δ 9 , of Ref. [5], in the most technically difficult groups: the ones that contain three photons, all coupling to the main muon-electron line, one example of which is shown in Fig. 2(a); and diagrams containing one or two virtual loops with mass M , as shown in Fig. 2(b,c).
In principle, an extension of our calculation to diagrams containing virtual m loops and massless loops, and threegluon vertices is possible. Instead, in the following Section III, we show how one can reproduce the first two terms of all types from already known form factors.

A. Hidden symmetry under M ↔ m
The starting point of the expansion in the difference of masses is the equal mass limit. There, the three loop correction has been known for a long time [17]. That result can be used to obtain the first two terms in the mass difference expansion, as explained in Ref. [19].
Without radiative corrections, the effect of the electron mass on the muon decay rate is known exactly, 192π 3 is the decay rate in the limit of a massless electron and ρ = m/M is the ratio of electron and muon masses. Near the limit of equal electron and muon masses a useful parameter is δ = (M − m)/M = 1 − ρ, introduced in Section II. The first two terms in the expansion of the tree-level decay rate around δ = 0 are We see that the phase space suppresses the decay rate by five powers of δ.
The following discussion applies equally to the muon decay and to the semileptonic b → c and b → u quark decays. Quark decays require a more extensive calculation because QCD corrections have a richer gauge group structure. We therefore use the b → c decay as an example, use the language of QCD, and denote δ = (m b − m c )/m b . Corrections to the muon decay can be deduced from a subset of results obtained for quarks.
We want to argue that when radiative corrections are included, the first two orders in δ, displayed in Eq. (7), are affected only by virtual corrections and not by the real radiation. Real radiation emission is suppressed by two powers of the daughter fermion velocity in the rest frame of the decaying fermion. That velocity, v = 1 − , is at most about δ. So, the real radiation affects the decay rate only in the third order in the expansion, δ 7 .
It has also been demonstrated that the form factors arising from virtual corrections do not have a linear term in the δ expansion [19], as long as the coupling constant α s is renormalized in the symmetric point µ 2 = m c m b . This is because of the symmetry of the virtual gluon diagrams under the interchange m b ↔ m c . We have therefore used Eq. (7) to reproduce the first two terms of the expansion in [5].
The available information about the form factors in [17] is limited to the axial form factor. Fortunately, in the limit of equal masses, the vector form factor does not receive radiative corrections. Corrections are suppressed by the square of momentum transfer q 2 = (p b − p c ) 2 which contributes again only in the third order, δ 7 .
Since the results of [17] assume the exact equal mass limit, they do not differentiate between the two heavy fermions. In Ref. [5] effects of b and c loops are labeled with factors n b and n c , to indicate that they are proportional to the number of those quarks in virtual loops (in reality of course there is n b = n c = 1). We thus replace n b and n c in the formulas of [5] by N H /2, so that the number of heavy fermions N H in [17] corresponds to n b + n c .
Finally, when we apply Eq. (7), we obtain corrections in terms of the coupling constant normalized at the geometric mean of the decaying and daughter fermion. On the other hand, results of Ref. [5] use the mass of the decaying fermion as the renormalization scale. To remedy this, we run the coupling constant in the formulas resulting from Eq. (7) using [32], where n f = n l + N H . In the above formulas, δ originates from the logarithm ln m b mc = − ln(1 − δ) = δ + O δ 2 . After this adjustment of the coupling constant in the axial form factor, we reproduce all coefficients of δ 5 and δ 6 in Ref. [5].
B. Example of the pattern in terms δ 5 and δ 6 As an example of the procedure outlined above, we consider the first two terms in the color part proportional to where a 4 = Li 4 (1/2) ≃ 0.517 [33] and ζ n is the Riemann zeta function [34]. Following the discussion around Eq. (7), the bulk of these terms should have the form d 5 δ 5 (1−3δ/2). Consider the deviation from that pattern, D = d 6 +3d 5 /2, a much simpler expression than either d 5 or d 6 , but not zero. We now show how this difference is removed by changing the renormalization scale. To this end, consider one-and two-loop form factor [18,[35][36][37], as summarized in [17], η (1) where we have neglected other color structures in η (2) A . Change of the scale of α s in the square of the form factor produces an additional contribution α 3 where the factor 48/5 is related to the normalization of the results in [5]. The result of changing the scale from the symmetric point m b m c to m 2 b , given by Eq. (18), provides precisely the deviation from the pattern d 5 δ 5 (1 − 3δ/2) we found in Eq. (13).

IV. CONCLUSIONS
The third order correction to the muon and the heavy quark decay rates found in Ref. [5] is an important milestone in perturbative quantum field theory. Here we have confirmed the first five terms of the three technically most difficult structures involving no light vacuum polarization loops and no three-gluon couplings. We have also revealed a pattern in the first two terms of the expansion of all flavor and color structures, hidden when the symmetry between the initial and the final fermion is broken by renormalization. This pattern becomes visible when a symmetric scale M m is used to renormalize the coupling.
We close by congratulating the authors of Ref. [5] on completing their heroic calculation.