A subtraction scheme for massive QED

We present an extension of the FKS subtraction scheme beyond next-to-leading order to deal with soft singularities in fully differential calculations within QED with mas- sive fermions. After a detailed discussion of the next-to-next-to-leading order case, we show how to extend the scheme to even higher orders in perturbation theory. As an application we discuss the computation of the next-to-next-to-leading order QED corrections to the muon decay and present differential results with full electron mass dependence.


Introduction
One of the crucial ingredients needed for perturbative calculations in QED and QCD is a method to perform the usually divergent phase-space integration. Since we are often interested in matching experimental procedures as closely as possible, it is essential to be able to compute fully differential cross sections with several cuts applied. This renders the phase-space integration too complicated for analytic evaluation.

JHEP01(2020)085
One possibility to perform the divergent phase-space integrations numerically is through universal infrared subtraction schemes. These methods allow the calculation of the real corrections at next-to-leading order (NLO), at least in principle, for any process in QED or QCD. This is accomplished by the construction of a counterterm that will point-wise subtract the singularities in the integrand such that a numerical integration over the phase space is possible. The counterterm has to be in a form that allows analytic integration. By exploiting the universal structure of soft and collinear singularities, these counterterms can be constructed in a process-independent way.
Two of the most widely used truly universal schemes at NLO are FKS [1,2] and the dipole formalism [3,4]. FKS treats soft and collinear singularities separately by constructing different local counterterms for each. The dipole formalism does away with this distinction by using Lorentz invariant dipole terms that remove both singularities simultaneously.
Recent years have seen a huge effort towards the development of universal schemes for next-to-next-to-leading order (NNLO) calculations (see e.g. [5][6][7][8][9][10][11][12]). Typically, these schemes were designed with QCD calculations in mind. This is both their greatest strength and their greatest weakness: by being applicable to non-Abelian theories, they all sacrifice simplicity to some extent. This makes them awkward to use for calculations in massive QED because they effectively treat collinear singularities that are not present.
Hence, we present the subtraction scheme FKS 2 that, while limited to massive QED, does not suffer from this problem and is very simple to implement. FKS 2 is a natural extension of the FKS scheme to NNLO for double-soft singularities. It can also readily be extended to even higher orders in perturbation theory.
Even though the electromagnetic coupling α is much smaller than the strong coupling, two-loop calculations in QED can be important in cases where a very high precision is required. As an example we mention Bhabha scattering (see [13] and references therein), which is used for a luminosity measurement at lepton colliders. Another potential application is muon-electron scattering which can be used for an alternative determination of the hadronic contribution to the running of α [14]. Recently, also NNLO corrections due to emission from the electron line for electron-proton scattering have been computed [15].
In addition to fixed-order contributions it is often also required to include multiple soft and/or collinear emissions of photons. Typically this is done by combining a parton-shower approach with fixed-order computations [16][17][18][19][20][21][22]. This allows to resum logarithmically enhanced terms. A recent example where resummation is combined with fixed-order NNLO contributions in QED can be found in [23]. While we will not address resummation directly, it is important to keep it in mind when constructing a subtraction scheme. This paper is structured as follows: in section 2 we briefly review the FKS scheme at NLO to introduce our notation and familiarise the reader with the necessary concepts. Next, in section 3 we present FKS 2 , the extension of the scheme to NNLO. Further extensions beyond NNLO are discussed in section 4, while referring details of the N 3 LO case to appendix B. In addition, we comment on some generic properties of the scheme in section 5. Next, we use section 6 to demonstrate the validity of FKS 2 by calculating as an example the NNLO QED corrections to the muon decay in a fully differential way. Finally, we conclude in section 7.

JHEP01(2020)085 2 Notation and concepts
Beyond leading order, a physical cross section is computed as a sum of several separately divergent parts. As a concrete example we consider a NNLO contribution to a n-particle process, which can be written as The double-virtual corrections are obtained by integrating M (2) n over the Born phase space dΦ n . Here M (2) n contains all terms of the n-particle (renormalised) matrix element squared with two additional powers of the coupling α. This includes the interference term of the two-loop amplitude with the tree-level amplitude as well as the one-loop amplitude squared. Similarly, the real-virtual contribution is obtained by integration of M (1) n+1 , the interference of the (renormalised) (n + 1)-particle one-loop amplitude with the corresponding treelevel amplitude, over the (n + 1)-particle phase space dΦ n+1 . Finally, for the doublereal contribution the tree-level matrix element with two additional particles, M (0) n+2 , is integrated over the corresponding phase space. In (2.1) we implicitly assume the presence of the flux factor (or the analogous factor for a decay rate) as well as a measurement function that defines the observable in terms of the particle momenta. The measurement function has to respect infrared safety, i.e. the observable it defines must not depend on whether or not one or more additional soft photons are present as arguments of this function.
In order to have sufficient flexibility in defining the observable, these phase-space integrals have to be done numerically. However, the presence of infrared singularities make a direct integration impossible. In dimensional regularisation with d = 4 − 2 this would lead to 1/ poles. Instead, a suitable subtraction has to be made such that the numerical integration is carried out only with expressions that neither contain implicit soft singularities from real emissions nor explicit 1/ singularities from loop integrations. Since we are dealing with QED processes with massive fermions there are no collinear singularities. Put differently, the collinear poles 1/ are replaced by ln(m) terms, where m is a fermion mass. This considerably simplifies the subtraction procedure and in what follows we present a scheme that is tailored to this situation.
The structure of soft singularities in QED has been studied a long time ago by Yennie, Frautschi and Suura (YFS) [24] to all orders in α. The key feature is that after splitting the amplitude (squared) into a contribution containing the soft singularity and a contribution free of soft singularities, the former exponentiate. Thus we can write where M ( )f n are free from infrared poles and all singularities are contained in the eikonal E. This does not yet completely defineÊ as there is some freedom to include finite terms. The precise definition we will use will be given in (2.11).

JHEP01(2020)085
The simple structure of (2.2) can be exploited to resum leading logarithmic corrections and even combine this with fixed-order computations. Indeed, there is a long history of using the YFS approach to construct Monte Carlo algorithms [20][21][22] to include QED effects in scattering processes. We will focus on a fixed-order approach and use the YFS formalism to extend the FKS subtraction scheme [1,2] to deal with soft singularities in QED beyond NLO.

FKS for soft singularities at NLO
In this subsection we will briefly summarise the necessary aspects of the FKS scheme at NLO. Because we only treat soft singularities, FKS is dramatically simplified. The NLO correction to a cross section is split into virtual and real parts The real corrections are obtained by integrating the tree-level matrix element M (0) n+1 over the phase space dΦ n+1 . To simplify the discussion we assume that in the tree-level process described by M (0) n no final-state photons are present. Hence, in M (0) n+1 only the particle (photon) with label n + 1 can potentially become soft. If there are additional photons (i.e. photons in the tree-level process) the combinatorics becomes slightly more involved, but the essential part of the discussion is not affected.
When computing a cross section in the centre-of-mass frame, we choose coordinates where the beam axis is in z direction. Further, we denote the (partonic) centre-of-mass energy by √ s. When computing a decay width we instead parametrise one of the outgoing particles in z direction and, if necessary, rotate the coordinate system afterwards.
Following [1] we parametrise the momentum of the additionally radiated particle n + 1 as 1 where e ⊥ is a d − 2 dimensional unit vector and the range of y 1 (the cosine of the angle) and ξ 1 (the scaled energy) are −1 ≤ y 1 ≤ 1 and 0 ≤ ξ 1 ≤ ξ max , respectively. The upper bound ξ max depends on the masses of the outgoing particles. Following [2] we find 1 Note that this parametrisation could also tackle initial-state collinear singularities because y corresponds to the angle between the photon and the incoming particles. However, a different parametrisation may be sensible (and is allowed here) to better account for pseudo-collinear singularities from light particles (cf. section 5.3). What is important in the following is that the scaled energy ξ1 is chosen as a variable in the parametrisation to ensure a consistent implementation of the distributions defined in (2.8).

JHEP01(2020)085
Further kinematic constraints are assumed to be implemented through the measurement function. We write the single-particle phase-space measure for particle n + 1 as where the d − 3 angular integrations dΩ and trivial factors are collected in dΥ 1 . Denoting by dΦ n,1 the remainder of the (n + 1)-parton phase space, i.e. dΦ n+1 = dΦ n,1 dφ 1 , we write the real part of the NLO differential cross section as To isolate the soft singularities in the phase-space integration we use the identity in terms of distributions. Here we have introduced an unphysical free parameter ξ c that can be chosen arbitrarily [1,2] as long as The dependence of ξ c has to drop out exactly since no approximation was made. Therefore, any fixed value could be chosen. However, keeping it variable is useful to test the implementation of the scheme.
Using (2.8) we split the real cross section into a hard and a soft part 2 In dσ (1) s we can now (trivially) perform the ξ 1 integration. To do this systematically, we define for photons the general soft limit S i of the i-th particle

JHEP01(2020)085
where M (0) m−1 is the matrix element for the process without particle i. The eikonal factor is assembled from self-and mixed-eikonals and sign jk = (−1) n jk +1 , where n jk is the number of incoming particles or outgoing antiparticles among the particles i and j. Further, we define the integrated eikonal completing the definition in (2.2). After dΥ 1 and dξ 1 integration (under which dΦ n,1 → dΦ n ) we obtain This part now contains explicit 1/ poles that cancel against poles in the virtual cross section. The second term of the real corrections, dσ (1) h given in (2.9c), is finite and can be integrated numerically after setting d = 4. Combining the real and virtual corrections, the NLO correction is given by We have used M  .13). According to (2.2) the explicit 1/ poles cancel between the two terms in the integrand of (2.13b) and the phase-space integration in (2.13c) is also manifestly finite.
We note that S i is invariant under rotations, but not Lorentz invariant, because it contains the explicit energy E i . Hence, also E i andÊ are only invariant under rotations but not under general Lorentz transformations. The integrated eikonalÊ jk has been computed in [2], dropping terms of O( ). As we will see this is sufficient even beyond NLO. The expression is given in section A, using our conventions.

FKS 2 : NNLO extension
In the following, we discuss the extension of FKS to NNLO, while still limiting ourselves to massive QED. To simplify the discussion in this section, we assume that all (suitably renormalised) matrix elements are known to sufficient order in the coupling and expansion in . In section 5.1 we will state what precisely is needed for a NNLO computation.

Real-virtual correction
The treatment of the real-virtual contribution proceeds along the lines of normal FKS because it is a (n + 1)-particle contribution. Again we assume that there is only one external particle, with label n + 1, that can potentially become soft. We use (2.8) with another unphysical cut-parameter ξ c A to split the realvirtual cross section into a soft and a hard part For dσ (2) s the analogy to the NLO case is particularly strong because there is no genuine one-loop eikonal contribution [25,26], i.e. the soft limit of the real-virtual matrix element is with the same E n+1 as in (2.10). Therefore, compared to (2.12) the definition of the soft part remains essentially unchanged However, dσ n . Unfortunately, dσ (2) h is not yet finite as it contains an explicit 1/ pole from the loop integration. To deal with this pole we use eikonal subtraction, i.e. we split the real-virtual matrix element according to (1)f n+1 is free from poles. This is again the YFS split, mentioned in (2.2). In (3.4) we have introduced yet another initially independent cut-parameter ξ c B .
With the help of (3.4) we can now write where c A indicates that the subtraction should be performed with the cut parameter ξ c A . The finite piece

JHEP01(2020)085
can be integrated numerically with = 0. Integrating the divergent piece, dσ (2) d , over the complete phase space we obtain dσ (2) where in I the first argument refers to the cut-parameter of the ξ integration and the second to the argument ofÊ. This process-and observable-dependent function is not finite and generally very tedious to compute. Even for the simplest cases it gives rise to complicated analytic expressions including for example Appell's F i functions. However, as we will see it is possible to cancel its contribution exactly with the double-real emission.
To summarise, the real-virtual corrections are given by where the expressions for dσ (2) s , dσ (2) f and dσ (2) d can be read off from (3.3), (3.5b) and (3.5c), respectively. We point out that dσ (2) rv is independent of both ξ c A and ξ c B .

Double-real correction
For the double-real contribution we have to consider M (0) n+2 , the matrix element for the process with two additional photons (with labels n + 1 and n + 2) w.r.t. the tree-level process. We extend the parametrisation (2.6) accordingly to Writing the phase space as dΦ n+2 = dΦ n,2 dφ 1 dφ 2 , the double-real contribution becomes where we have used analogous definitions as in (2.6) and (2.7). The only difference between dΦ n,1 and dΦ n,2 is in the argument of the δ function that ensures momentum conservation. Note that the factor 1/2! is due to the symmetry of identical particles.
For the mixed contributions dσ (2) hs and dσ (2) sh we use Considering first dσ (2) hs , we perform the ξ 2 integration (under which dΦ n,2 → dΦ n,1 ) and use (2.11) to do the dΥ 2 integration to obtain dσ (2) hs (ξ c 1 , ξ c 2 ) = dΥ 1 dΦ n,1 Thus, we find again the integral I of (3.5c). Finally, we turn to the double-soft contribution dσ (2) ss . Since the ξ integrals in dσ (2) ss factorise. Therefore, we can do the ξ 1 , Υ 1 integrations independently from the ξ 2 , Υ 2 integrations and obtain dσ (2) ss (ξ c 1 , ξ c 2 ) It is clear that the simplicity of the infrared structure of QED with massive fermions is crucial for reducing the complexity of the procedure described in the steps above.

Combination
At this stage we have introduced four different cutting parameters ξ c A and ξ c B as well as ξ c 1 and ξ c 2 . All of these are unphysical, arbitrary parameters that can take any value 0 < ξ c i ≤ ξ max . In total we have to deal with seven different contributions. Two of them, dσ (2) s and dσ (2) ss , are very simple as they just depend on the eikonal. Another two contributions dσ (2) f and dσ (2) hh can be calculated numerically with = 0. The sum of the three remaining auxiliary contributions dσ (2) d , as well as dσ (2) sh and dσ (2) hs , only depend on the function I defined above Note that, due to the sign difference and the symmetry factor, dσ (2) aux vanishes if we choose (3.14) This cancellation will not be affected by the measurement function. Thus, in what follows we will make the choice (3.14), avoiding the computation of the potentially difficult I function. 3 We can now collect the non-vanishing contributions, sorted by remaining integrations (3.15d) The three terms of the integrand of σ (2) n are separately divergent. However, in the sum the 1/ poles cancel. The other parts, σ (2) n+1 and σ (2) n+2 , are finite by construction. Hence, we can set d = 4 everywhere (except in the individual pieces of the integrand of σ (2) n ) and obtain

JHEP01(2020)085
This is the generalisation of (2.13) to NNLO. In the integrand of (3.16a) the build up of the exponentiated singular part e αÊ is recognisable. For M ( )f n to be finite,Ê has to contain the soft 1/ pole. However, any choice of the finite part is possible in principle. We have chosen to define the finite matrix elements through eikonal subtraction, (3.4). This ensures that the auxiliary contributions cancel and the remaining parts σ First steps towards extending universal schemes beyond NNLO have been made in QCD [27]. The simplicity of FKS 2 suggests that this paradigm is a promising starting point for further extension to N 3 LO in massive QED, provided that all matrix elements are known.
At N 3 LO, we have four terms which are separately divergent. In order to reorganise these four terms into individually finite terms, we repeatedly use ( all terms are separately finite and, as discussed in detail in appendix B, given by Once more we have used the fact that for tree-level amplitudes M As always, the ξ c dependence cancels between the various parts such that dσ (3) is independent of this unphysical parameter.

FKS : extension to N LO
The pattern that has emerged in the previous cases leads to the following extension to an arbitrary order in perturbation theory: The eikonal subtracted matrix elements

Comments on and properties of FKS
With the scheme now established, let us discuss a few non-trivial properties that are helpful during implementation and testing.

Regularisation-scheme and scale dependence
In QED calculations it is advantageous to calculate the matrix elements M ( ) n in the on-shell scheme for α (and the masses). This way the only µ dependence is in a global prefactor µ 2 induced through the integral measure. The same holds for the integrated eikonal. Hence, for the finite matrix elements M ( )f n there is no µ dependence after setting d = 4. A similar argument can be made for the regularisation-scheme dependence. So far we have implicitly assumed that the computation is performed in conventional dimensional regularisation. However, it is often more convenient to use other dimensional schemes, where e.g. external particles are treated in four dimensions [28]. As discussed in [29,30], after proper renormalisation and infrared subtraction the matrix elements are regularisationscheme independent for d = 4. As there are no collinear singularities, the eikonal-subtracted M ( )f n are scheme independent. Moreover, the integrated eikonal is scheme independent. To be concrete, we list the input that is required for a computation of a physical cross section at NNLO in QED. The important point is that once the final expressions for a NNLO cross section, (3.16), or beyond, (4.3), are obtained, we can set d = 4 everywhere.
• The two-loop matrix element M (2) n is known with non-vanishing masses up to O( 0 ). In general this is a bottleneck because the necessary master integrals are only known for a very select class of processes, not to mention the algebraic complexity. However, it is possible to approximate dσ (2) using 'massification' of M

JHEP01(2020)085
• The renormalised one-loop matrix element M (1) n of the n-particle process is known including O( 1 ) terms. This is usually the case for NNLO calculation as it is needed for the sub-renormalisation M (2) n ⊃ δZ × M (1) n as well as the one-loop amplitude squared, which is part of M (2) n . Once these pieces are assembled to M (2)f n , the O( ) terms can be dropped.
• The renormalised real-virtual matrix element M (1)f n+1 is assembled, the O( ) terms can be dropped.

Massification
The bottleneck in the computation of cross sections for massive QED at N LO is the availability of the matrix element M n (m = 0) into log(m) terms that will cancel against corresponding 'singularities' of the real corrections. In addition, it retains in dσ ( ) n the finite log(m) terms that are present in differential distributions. However, terms m log(m) that vanish in the limit m → 0 will be neglected.
It should be noted that a similar procedure in dσ

Phase-space parametrisation
A further issue in connection with small lepton masses is related to the phase-space parametrisation. The phase space has to be constructed in any way that allows the distributions to be well implemented, i.e. ξ i should be an integration variable of the numerical integrator. In addition, for small m there are potentially numerical problems due to pseudo-collinear singularities. In fact, these regions produce precisely the log(m) terms that correspond to the collinear 'singularities' of the real part. These log(m) terms will cancel the virtual collinear 'singularities' mentioned above. Hence, for small m there is a numerically delicate cancellation. For the simple observables presented in section 6 this JHEP01(2020)085 problem could be solved through a dedicated tuning of the phase-space parametrisation. For more complicated processes another solution might need to be implemented in the future. The idea is to subtract such regions from the integrand and add them back in integrated form [35]. This amounts to actually extend the subtraction scheme to collinear singularities. It also offers the possibility to use massification for dσ ( ) n+j . To this end, M ( −j)f n+j (m = 0) is massified. The resulting collinear singularities due to phase-space integration with massless external fermions are treated as follows: first they are subtracted at integrand level; second their integrated contribution is added back in massified form.

Muon decay
As an example for the subtraction scheme presented in the previous section we will discuss the muon decay at NNLO in the Fermi theory of weak interactions. 4 NLO corrections to the muon decay have been known for many decades [36,37]. Using the optical theorem, the NNLO QED corrections to the decay width were calculated around the turn of the millennium assuming vanishing electron masses [38]. Over the course of the next decade, the electron energy spectrum, which is not infrared finite in the limit m e → 0, was calculated. At first, only its logarithms were known analytically [39,40]. A few years later, the full spectrum was calculated with a numerical loop integration [41] and the original calculation of [38] was extended to include mass effects [42]. It was only recently that the form factors necessary for a fully differential calculation were published [31,43]. We have included muon and electron loops but neither tau nor hadronic contributions [44,45]. We treat the electromagnetic coupling α in the on-shell scheme, except in table 1 where, in order to compare to [38] we need the MS couplingᾱ ≡ᾱ(µ = M ).

Calculation
The momenta of the muon and electron are written as Furthermore, β is the velocity of the electron in the muon rest frame. Apart from the form factors needed for dσ (2) n , we also need matrix elements for dσ (2) n+1 and dσ (2) n+2 . 5 We have generated the diagrams for M (1) n+1 and M (0) n+2 using QGraf [46] 4 To compare with literature we also have to consider the case of an {e + (p5), e − (p6)} pair in the final state. 5 Note that, to remain consistent with the discussion above, we will denote the decay by σ instead of Γ. f contribution with a fit (blue line). The orange line corresponds to dσ In green we show the sum dσ (2) h , (3.5a). The right panel is zoomed into the region of interest showing 1σ confidence bounds for the fit. All plots are normalised such that the average dσ (2) h (ξ c A ) = 1. and calculated them using Package-X [47]. The numerical integration of dσ (2) f and dσ (2) hh was performed in Fortran using vegas [48]. Most loop integrals in M (1) n+1 were included explicitly, while the more complicated triangle-and box-functions were evaluated using the COLLIER library [49].

ξ c dependence
Due to the simplicity of the process, it is actually possible to explicitly compute the contribution of the integral I(ξ c A , ξ c B ) of (3.5c) to the total decay rate, and check that the dependence of all four ξ c parameters vanishes for the physical result. To verify this we perform weighted two-dimensional fits of the form c 11 log ξ c i log ξ c j +c 10 log ξ c i +c 01 log ξ c j +c 00 for the numerical data of dσ (2) f and dσ (2) hh and check whether the ξ c dependences vanish within the numerical error of this fit.
The ξ c B (in)dependence (according to (3.5a)) of the combination dσ (2) h = dσ (2) f − I is shown in figure 1 for two different values of ξ c A . To numerically evaluate these expressions in the plots, we drop the 1/ poles consistently in all intermediate expressions. In a next step, dσ (2) h is then combined with dσ   Figure 2. ξ c A (in)dependence of the ( 0 part of the) real-virtual contribution. The combined hard contributions dσ (2) h (blue, dots and fit) and the soft contribution dσ (2) s (orange) sum up to the total real-virtual contribution dσ (2) rv (green), (3.6). The right panel is zoomed into the region of interest, showing 1σ confidence bounds to the fit. All plots are normalised such that dσ  Figure 3. The ξ c2 (in)dependence of the ( 0 coefficient of the) double-real corrections for two ξ c1 . The dσ (2) hh -data (blue, data and fit), dσ (2) ss (yellow) and 1 2 (I(ξ c1 , ξ c2 ) + I(ξ c1 , ξ c2 )) contribution (green) sum up to the total n + 2 particle contributions (red). The right panel magnifies the region of interest including 1σ bounds. Note that the swap c 1 ↔ c 2 is trivial. All plots are normalised such that dσ (2) rr (ξ c2 ) = 1.   Table 1. The different contributions to σ (2) . Note that σ (2) e also includes the process µ → ννe ee. See text for interpretation. The couplingᾱ(µ = M ) is renormalised in the MS scheme. ∆ rel denotes the relative difference of our results to the massless result [38]. JHEP01(2020)085

Results for the decay rate
The first quantity we consider is the full decay width where we have pulled out factors of the MS couplingᾱ/π. We compute σ (2) using the massified form factors as well as the form factor with full m dependence [31,43]. We will label these two results 'massified' and 'massive', respectively. In the case of the massified result, we expand all three parts of the integrand contributing to σ (2) n , see (3.15b) and (3.16a). Of course, the exact mass dependence of dσ (1) s and dσ (2) ss is usually much easier to obtain than for M (2) n . However, the complete cancellation of singularities requires a consistent expansion in z of all contributions at the n-particle level.
Because the full decay rate does not contain terms log z ∼ log m the limit m → 0 exists and we can compare our massified and massive results with the result for a massless electron [38]. We note that in this particular case (contrary to distributions, where log m terms exist), the massified result is not expected to be superior to the massless computation.
Following [38], we split the result into three parts: photonic corrections σ (1) γ and σ (2) γ , corrections due to an electron pair (real or virtual) σ (2) e , and corrections due to a muon pair (virtual) σ (2) µ . These parts have been defined and their analytic results in the massless case given in equations (2.11), (2.13) and (2.15) of [38]. The individual results for the NNLO corrections are shown in table 1, where the Monte Carlo error is smaller than the significant digits. Note that [38] had to include the 'open-lepton production' µ → ννe ee into their calculation of σ (2) e to guarantee finiteness. We have included this process as well [50] since it contributes to σ (2) e (two-trace contribution) and σ (2) γ (one-trace contribution). 6 The results of table 1 merit a few comments: • The good agreement for the purely photonic contributions σ (2) γ between the massive and massless result is due to the absence of terms log m and m log m as discussed by [38].
• The massified results differs by about 3% from the massive (and massless) result for σ (2) γ . This is due to the mismatch between the real corrections, that were calculated with the full electron mass dependence, and the massified two-loop amplitude that only includes logarithmically enhanced mass effects.
• The massified results agrees perfectly with [38] for the σ (2) µ part because the contribution comes purely from one two-loop diagram that is free of any soft or collinear logarithms and hence effectively massless.
• The massive and massless results for σ production. In the pure electron trace m must not be neglected to lead to finite expressions. However, in the other trace the electron mass can be set to zero. Our value of 3.16 was calculated with full electron mass dependence. If we were to set m → 0 in the this trace, we would obtain 3.23 in much better agreement with [38].
Note that in any case the 'massive' result should be considered the reference. Our results agree with [42]. For the pure mass effects of the photonic part, this agreement is only at the 20% level. This is due to large numerical cancellations between σ (2) n , σ (2) n+1 , and σ (2) n+2 which make the extraction of a few-percent effect on the NNLO corrections numerically challenging. In fact, an efficient numerical evaluation of the integrals with full mass dependence [43] has only recently been implemented [51].

The electron energy spectrum
In order to validate our computation, we consider the NNLO corrections to the normalised electron energy spectrum x e = 2E/M and compare them to results available in the literature. If two (negatively charged) electrons are present in the final state, we include both of them in the x e distribution. The leading and sub-leading logarithmic contributions for this observable were calculated in [39,40]. Because this corresponds to a strict expansion in z, we expect good agreement for large x e as noticed in [31]. In figure 4 we compare the two results and see that the differences are compatible with the constant (logarithm-free) terms missing in [39,40]. These terms were computed numerically and shown in a plot for x e > 0.3 in [41]. If we include these constant terms of [41], we obtain perfect agreement with our result, using the massive form factors. Note that the difference between massified and massive result in figure 4 is at the percent level and only becomes visible around the zero crossing at x e ≈ 0.21 and x e ≈ 0.88, never changing the overall picture. The on-shell coupling (α/π) 2 is omitted in the results shown in the figure 4.
With a fully differential Monte Carlo code, we can compute arbitrary distributions, including cuts. As an example, we consider again the normalised electron energy spectrum but impose a cut on photon emission. Concretely, we restrict the total energy of all photons within a cone of angle θ ≡ (p e , p γ ) = 37 • (i.e. a cone with | cos θ| > 0.8) around the electron to be less than 10 MeV.
The results are shown in figure 5. Comparing the normalised NNLO result (blue histogram) to the normalised LO result (orange line) in the top panel reveals that only for large x e the corrections to the shape are relevant. This is driven by the NLO corrections. They are large at both ends of the x e spectrum, as shown by the NLO K (1) factor -10 0 10 20 Figure 4. The NNLO corrections to the electron energy spectrum, omitting a factor (α/π) 2 . The logarithmic contributions of [39,40] (green dashed) agree reasonably well with our massive result (blue histogram) for large x e . Adding the constant terms of [41] (orange) we obtain good agreement.

Conclusion
We have presented a subtraction scheme tailored to the case of QED with massive fermions. This allows to perform the phase-space integration and obtain predictions for arbitrary physical cross sections if the corresponding matrix elements are known. While we have primarily NNLO calculations in mind, the extension of the scheme beyond NNLO is also discussed. After describing its implementation, we have commented on its properties. We have noted that, while intermediary results may look complicated, most of this complexity drops out in the final result. We have verified our scheme by calculating the muon decay at NNLO and compared with known results. While our scheme leads to very simple expressions, there are still several avenues for further developments. Since the analytic computation of matrix elements with massive fermions is very challenging, one possibility is to investigate the option of computing the subtracted matrix elements numerically. After all, these are finite expressions. If it is possible to perform the subtraction and UV renormalisation at the integrand level, a direct numerical evaluation in four dimensions should be feasible.
Even if the matrix elements are available, the numerical integration can be challenging. First of all, the matrix elements need to be implemented in a stable way in all corners of phase space. The presence of small fermion masses also results in pseudo-collinear singularities that can lead to numerical instabilities. An option to deal with these is to treat them as singularities, i.e. subtract them and add them back in integrated form. Obviously, this tarnishes the simplicity of the subtraction scheme. But it offers a possibility to take into account small masses in an expanded form for the double-real and real-virtual corrections. Such an implementation can be seen as a generalisation of massification for all parts of a NNLO cross-section computation.

JHEP01(2020)085
The scheme we have presented relies on eikonal subtraction and is closely related to the YFS exponentiation of the soft singularities. The YFS picture has been used to construct Monte Carlo codes that describe multiple emission of soft photons. Thus, our scheme for fixed-order computations lends itself to be combined with a YFS Monte Carlo program. Such a generalisation will allow to resum large logarithms and combine this with fixed order results. Ultimately, this is required to obtain very precise predictions for fully differential QED observables.

JHEP01(2020)085
For the case j = k this expression simplifies tô for the self-eikonals.
B Details for FKS 3

B.1 Real-virtual-virtual contribution
Let us begin with the real-virtual-virtual part that we split again into a hard and soft contribution as in (3.2). Using that even at the two-loop level the soft contribution in analogy to (2.12) and (3.3) is given by The hard contribution is now dσ (3) h (ξ c ) = dΥ 1 dΦ n,1 dξ where dσ (3) f is finite and the divergent part dσ  Above we have defined two functions I (1) and J that are potentially tedious to compute. However, as we will see they cancel in the final result, similar to the function I at NNLO.
Obviously dσ   Here we have defined a third auxiliary function K that will cancel in the final result.

B.3 Triple-real contributions
The evaluation of the triple-real contributions proceeds along the lines of the FKS 2 doublereal part, albeit with more (individually ξ c dependent) terms dσ (

B.4 Combination
Combining all contributions at N 3 LO we need to evaluate (4.1). Collecting the terms with an n-parton phase space we get dσ

JHEP01(2020)085
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.