Universal structure of radiative QED amplitudes at one loop

We present two novel results about the universal structure of radiative QED amplitudes in the soft and in the collinear limit. On the one hand, we extend the well-known Low-Burnett-Kroll theorem to the one-loop level and give the explicit relation between the radiative and non-radiative amplitude at subleading power in the soft limit. On the other hand, we consider a factorisation formula at leading power in the limit where the emitted photon becomes collinear to a light fermion and provide the corresponding one-loop splitting function. In addition to being interesting in their own right these findings are particularly relevant in the context of fully-differential higher-order QED calculations. One of the main challenges in this regard is the numerical stability of radiative contributions in the soft and collinear regions. The results presented here allow for a stabilisation of real-virtual amplitudes in these delicate phase-space regions by switching to the corresponding approximation without the need of explicit computations.


Introduction
The limit of scattering amplitudes where an external massless gauge boson becomes soft or collinear to a fermion has important applications from a conceptual quantum field theory point of view, but also for concrete applications in higher-order calculations. In this paper we consider QED with massive fermions and using a diagrammatic approach provide new explicit results of the limiting behaviour of one-loop radiative amplitudes with practical applications in mind.
In order to meet the experimental precision the calculation of next-to-leading order (NLO) or even next-to-next-to-leading order (NNLO) QED corrections has become mandatory. A lot of work has therefore been put into the computation of the corresponding loop amplitudes where the presence of non-vanishing fermion masses results in significant complexity. Furthermore, the need to account for non-trivial detector geometries and acceptances has triggered the development of numerous Monte-Carlo codes that allow for the calculation of fully differential observables.
Most of these calculations face the difficulty of a strong scale hierarchy due to small fermion masses acting as regulators of collinear singularities. It is therefore one of the main challenges of fully differential higher-order QED calculations to achieve numerically stable implementations of the amplitudes. Particularly delicate in this regard are radiative loop amplitudes evaluated in the phase-space region where the emitted photon is soft. These numerical instabilities are further exacerbated if the photon becomes collinear to a light on-shell fermion. The approximation of the corresponding expression using the soft expansion up to subleading power yields an elegant solution to this problem. This next-to-soft stabilisation has facilitated the first calculation of fully differential NNLO QED corrections to Bhabha [1] and Møller [2] scattering. Even though the soft expansion can straightforwardly be computed with the method of regions [3] the corresponding calculation is cumbersome. The currently limited knowledge of universal soft and collinear behaviour in QED makes therefore a systematic study highly desirable. Hence, in this article we present two universal results for radiative amplitudes at one loop, the next-to-soft limit and the collinear limit.
It is a well-known fact that in QED radiative amplitudes exhibit universal factorisation at leading power in the soft limit given by the eikonal approximation to any order in perturbation theory [4]. Furthermore, it has been shown a long time ago by Low, Burnett, and Kroll [5,6] that also the subleading term is related to the non-radiative amplitude at tree level via a differential operator. 1 This so-called LBK theorem was later extended to massless particles [10] where a universal radiative jet function was introduced to take into account collinear effects. More recently, the massless version of the theorem has attracted some attention in the context of resummation of next-to-leading power threshold logarithms. To this end the theorem was extended to also include loop corrections in the framework of diagrammatic factorisation [11][12][13] as well as in soft-collinear effective theory (SCET) [14][15][16]. In this paper, however, we are interested in the case of QED where all fermion masses and all other scales are considered to be much larger than the energy of the emitted photon. In this case, these recent loop-level extensions are not applicable since the underlying effective theory is heavy-quark effective theory (HQET) and not SCET. In particular, there is no radiative jet function in this scenario due to the absence of any collinear scale. This leaves hard and soft modes as the only relevant degrees of freedom.
The collinear limit of radiative amplitudes has been extensively investigated in the context of QCD with massless quarks where it gives rise to infrared singularities. The factorisation into a process-independent splitting function multiplying the non-radiative amplitude has therefore been known for some time now [17,18]. While the splitting functions correspond to the Altarelli-Parisi kernels at tree level, this is no longer true if loop corrections are taken into account. The two-loop corrections to the QCD splitting functions have been calculated in [19,20]. Unfortunately, much less is known in the case of QED where collinear divergences are regularised by finite fermion masses. Taking these masses to be small, we expect a similar factorising structure. This is due to the applicability of SCET to the case of small but non-vanishing fermion masses. Nevertheless, the corresponding splitting function is currently only known at tree level where it coincides with the QCD version up to an additional mass term [21][22][23][24].
It is therefore the goal of this paper to present the one-loop generalisation of both, the LBK theorem as well as the small-mass collinear factorisation formula in QED. The paper is organised as follows: We start by giving the most important conventions and definitions in Section 2. Section 3 then discusses the soft limit at subleading power while we present our findings about the leading power collinear limit in Section 4. In both cases we first give a brief summary of the tree-level derivations before presenting the one-loop extension. The main results of this paper are given in (3.31), (4.22) and (4.23). In Section 5 we then apply and validate these formulas in the process e − e + → e − e + γγ. Finally, we comment on possible applications and future developments in Section 6.

Notation and conventions
We denote the amplitude for a process with n final state particles by A n and the corresponding QED l-loop correction by A (l) n . Analogously, we use M n and M (l) n for the unpolarised squared amplitude and refer to it as matrix element. The correponding quantities for the radiative process with one additional photon in the final state are given by n → n + 1. The perturbative expansion of the amplitudes and matrix elements is done in terms of the fermion charge Q where Q = −e for an incoming particle or an outgoing antiparticle and Q = +e otherwise. The symbol Γ is used in various places to generically parametrise part of an amplitude.
In what follows we always denote the loop momentum by ℓ. Furthermore, we consistently take k to represent the momentum of the emitted photon in a radiative process. It is the soft and collinear limit of this particle that is studied in the following sections. Other on-shell momenta are denoted by p i , i.e. p 2 i = m 2 i with m i the particle mass. The corresponding velocity is then given by In the situation where a particle may also be off shell we use q i instead.
To make the power counting in the soft expansion transparent we introduce the soft bookkeeping parameter ξ. The soft limit is then governed by the scaling k ∼ ξ ≪ m i ∼ ξ 0 which yields the expansion with the leading and subleading power contributions M sLP n+1 and M sNLP n+1 . Similarly, we introduce the book-keeping parameter λ for the power counting in the collinear limit. The corresponding expansion is only meaningful if the collinear fermion is light. We therefore take the mass to scale as p 2 = m 2 ∼ λ 2 . This in turn gives the behaviour k · p ∼ λ 2 due to At the level of the matrix element this gives rise to the expansion M n+1 In the collinear limit the amplitudes scale as A n+1 ∼ λ −1 and not as ∼ λ −2 as one would naively expect.
The parameters ξ and λ are just used to facilitate the book keeping of the corresponding expansions. In all our equations their numerical value is ξ = λ = 1. Furthermore, we emphasize that the different treatment of the fermion masses in the two limits has crucial implications. The small mass in the collinear limit introduces a collinear scale that is absent in the soft limit. From a more formal effective field theory perspective, the collinear case is therefore governed by SCET while the soft behaviour can be described in terms of HQET.
It is well known [4] that the leading-power term in the soft expansion is given by to all orders in perturbation theory with the eikonal factor summing over all external legs. It is the purpose of this paper to present similar universal results for the subleading term M sNLP n+1 (Section 3) as well as the leading collinear contribution M cLP n+1 (Section 4).

Soft factorisation at subleading power
It has been shown a long time ago by Low, Burnett, and Kroll that the soft expansion of radiative tree-level amplitudes can be related up to subleading power to the non-radiative process by means of a differential operator [5,6]. In this section we will generalise this LBK theorem to the one-loop level. We first start in Section 3.1 by giving a short review of the tree-level derivation. The oneloop extension is then discussed in detail in the following Section 3.2 with the main result given in (3.31).

The LBK theorem
Following [6] we split the radiative tree-level amplitude A (0) n+1 into contributions due to external and internal emission In addition to the set of on-shell momenta {p} = {p 1 , ... , p i , ... , p n } we define the sets of momenta {p} i = {p 1 , ... , p i − k, ... , p n } that are adapted to emission from line i. Taking all particles apart from the emitted photon to be incoming (but ignoring the complex conjugation of the polarisation vector ǫ) allows us to write the soft expansion of A ext n+1 as Since {p} satisfies the radiative momentum conservation i p i = k this is not a strict expansion in ξ. Following [25] we can make the above split gauge invariant (up to subleading power) via the modification . The leading contributions in ξ vanish due to i Q i = 0 and the subleading contributions cancel between the two terms of the last expression in (3.3a). Because the full amplitude is gauge invariant we also have k · A II n+1 ∼ O(ξ 2 ). Following the reasoning in [25], we observe that the leading O(ξ 0 ) term in A II n+1 must be independent of k due to the lack of 1/k poles in A int n+1 . Combined with the Ward identity this directly implies that A II n+1 ∼ O(ξ). As a consequence, the soft expansion of the total amplitude can be written as with the LBK operator Squaring the amplitude, summing over spins and polarisations, and using the identity then yields for the matrix element This shows that not only the leading term in the soft expansion is related to the non-radiative matrix element but that this is also true at subleading power at tree level. However, the nonradiative matrix element in (3.6) is evaluated with a set of momenta {p} that does not satisfy momentum conservation. This is unproblematic for the tree-level matrix elements considered here. If, on the other hand, loop corrections are taken into account (Section 3.2) this significantly complicates the evaluation of the corresponding integrals. In this case a different formulation of the LBK theorem is helpful. To this end, we reabsorb the first term of the LBK operator into the matrix element to undo the expansion and write Since {p} i satisfies momentum conservation we can now express the non-radiative matrix element in terms of invariants The corresponding expansion in k can then be written as where the sum L is over the set of independent invariants {s} = s({p}, {m 2 }) expressed in terms of the momenta {p} and the on-shell masses {m 2 }. Similarly, we can write Inserting (3.9) and (3.10) into (3.7), all derivatives with respect to the masses cancel and we obtain the simple formulation of the LBK theorem in terms of invariants with the modified LBK operator The advantage of (3.11) over (3.6) is that conventional one-loop techniques can be applied in this case. We emphasise that the choice of {s} = s({p}, {m 2 }) is ambiguous since the momenta {p} do not satisfy momentum conservation. This is however not an issue as long as the same definition is used in the calculation of the derivatives ∂s L /∂p µ i . The above formula can therefore be conveniently used to analytically compute the soft limit of tree-level matrix elements up to subleading power. An alternative approach that is particularly suitable for the numerical evaluation of the LBK theorem was recently presented in [26].
The above formula assumes all particles apart from the photon to be incoming. In the case of outgoing particles the corresponding momentum p has to be replaced with −p. In particular, this also implies ∂/∂p µ → −∂/∂p µ .

One-loop generalisation of the LBK theorem
The derivation of the previous section cannot be applied one-to-one in the presence of loop corrections due to contributions from regions where the loop momentum is small (ℓ ∼ ξ). In this soft region the power counting used to derive the LBK theorem is no longer valid. In particular, individual contributions from internal emissions already contribute at leading power. Fortunately, the method of regions [3] can be used to disentangle the soft contribution from the hard one. For the latter the LBK formula (3.6) still holds. The soft contribution will be evaluated in a generic way in the following sections. The collinear regions, on the other hand, all vanish due to the absence of collinear scales in the soft limit as defined in Section 2. The combination of hard and soft contributions, that will be given in Section 3.2.4, therefore generalises the LBK theorem to one loop.

General considerations regarding the soft contribution
In order to systematically analyse possible origins of the soft contribution we classify the one-loop integrals as illustrated in Figure 1 where the momenta q i can be off shell (internal) or on shell (external). The circle symbolises the one-loop integral associated with the one-particle irreducible (1PI) part of a particular Feynman diagram. The first class, L 1 , includes (m + 1)-point integrals from diagrams where the photon is directly attached to this 1PI part. Class L 2 includes m-point integrals from diagrams where the photon is attached to a leg that directly connects to the 1PI part with momentum q 1 − k. As we will see, the treatment of these integrals depends on whether the momentum of the adjacent leg q m is on shell or off shell. Finally, for integrals of the type L 3 the photon is attached indirectly to the m-point 1PI part such that the momentum flowing into the loop integral is q 1 − k +q with a non-zeroq. For integrals to have a non-vanishing soft contribution the momentum routing has to be chosen such that the loop momentum ℓ is aligned with a photon propagator. All other choices lead only to linear propagators in the soft momentum expansion and therefore vanish as a consequence of the residue theorem. There can thus be at most as many soft regions as the number of photons in the loop. However, most of them yield scaleless integrals and vanish in dimensional regularisation. This is in particular the case for all possible routings of L 3 . The presence of the momentumq allows to set ℓ = 0 for the soft contribution in all propagators except for the photon propagator with momentum ℓ. Hence, loop integrals of the form L 3 do not contribute to the soft region. For the second class, on the other hand, there is one non-vanishing soft contribution indicated by the momentum routing in Figure 1b if the corresponding internal propagator is given by a photon and if in addition q 1 is on shell. In the case where q m is off shell the soft expansion starts at O(ξ) and is given by .

(3.12)
For on-shell q m the leading integral reads instead which already contributes at O(ξ 0 ). Finally, the first class of loop integrals L 1 gives rise to up to two non-vanishing soft contributions given by the two momentum routings ℓ andl in Figure 1a if the corresponding propagators are photons. The integral for routing ℓ is only non-zero if q 1 is on shell and it starts to contribute at O(ξ −1 ) if q m is on shell and at O(ξ 0 ) otherwise. The analogous statements hold for thel momentum routing. The above reasoning allows to represent every possible soft contribution according to the three pairs of diagrams shown in Figure 2 where the external legs are now all on shell. Every R int, ext {e,a} corresponds to an amplitude with a specific choice of the momentum routing where the soft contribution does not vanish. The labels for emission and absorption {e, a} take on the values 1 . . . n or Γ. The superscript int, ext indicates whether the photon k is attached internally or externally. In the former (latter) case we are dealing with integrals of the type L 1 (L 2 ). As mentioned in connection with L 1 , it is possible that one amplitude contributes to two soft representations. Taking e = i and assuming p i to be incoming, we can write the corresponding expressions generically as All terms related to the emission from leg e = i and the soft photon propagator are given explicitly in (3.15). The vertex and fermion propagator related to the absorption is common to R ext {e,a} and R int {e,a} and is included in Γ µ {e,a} . This implies the scalings Γ µ {i,j} ∼ Γ µ {i,i} ∼ ξ −1 and Γ µ {i,Γ} ∼ ξ 0 . We then write the expansion in the soft region of the sum of the diagram pairs as with the leading and subleading power terms denoted by S LP {e,a} and S NLP {e,a} , respectively. Based on the previously discussed power counting of the integrals L 1 and L 2 we can deduce that S LP {i,Γ} = S LP,ext {i,Γ} = S LP,int {i,Γ} = 0.

Vanishing of the soft contribution at leading power
The leading-power soft contribution to (3.15) is given by The propagators can be brought to the same form with the partial fraction identity where the second term in the curly brackets can be neglected up to scaleless integrals. We then see immediately that We therefore arrive at the known result that the eikonal approximation in QED does not receive genuine loop corrections. Furthermore, this also shows that S NLP {i,Γ} = 0 since it effectively corresponds to a leading-power contribution. The third class of soft contributions, R soft {i,Γ} , can therefore be omitted in the following discussion.

Soft contribution at subleading power
At subleading power there are contributions in (3.15) from either the higher-order expansion of propagators (denominator) or from numerator terms proportional to ℓ or k. We therefore write For the denominator type the leading-power cancellation of Section 3.2.2 occurs if propagators other than [ℓ 2 + 2ℓ · (p i − k) − 2k · p i ] or [ℓ 2 + 2ℓ · p i ] are expanded. Furthermore, expansion in ℓ 2 of these two propagators results only in linear propagators. Consequently, we only have to consider the expansion in ℓ · k of the propagator [ℓ 2 + 2ℓ · (p i − k) − 2k · p i ] in (3.15). Using partial fraction then yields the simple contribution The numerator type can be written as

22c)
where we have used the replacement ℓ · p i → k · p i in the numerator which holds up to scaleless integrals.
To make further progress at this point we need to treat S NLP {i,j} and S NLP {i,i} separately in order to specify the form of Γ µ {e,a} . In the case where the photon is reabsorbed by the emitting leg, i.e.
where A (0) n = Γ (0) u(p i ) corresponds to the non-radiative tree-level amplitude. In this case, we further have the simple Passarino-Veltman decomposition where we can again replace ℓ · p i with k · p i . It is then straightforward to see that Hence, diagrams where the loop corrects only the emitting leg do not contribute at subleading power.
In the case of S NLP {i,j} , on the other hand, we find for an incoming particle p j that This in turn results after the tensor decomposition in the only non-vanishing subleading power contribution of the form where we have defined the function The analytic results for the integrals , (3.29c) can be found in Appendix A. The causal +iδ prescription is given explicitly in the above integrals. The result (3.29) is also valid for incoming antiparticles with the overall sign difference parametrised by the fermion charges Q i and Q j . The total soft contribution can thus be obtained by summing the above expression over all external charged fermions, i.e.
The corresponding expression for the matrix element can be obtained by interfering with the eikonal approximation of the tree-level amplitude. The resulting formula is given in the following section.

One-loop extension of the LBK theorem
Based on the previous discussion we find that the one-loop correction to a generic radiative process in QED in the limit where the emitted photon becomes soft satisfies the expansion This is the generalisation of the LBK theorem at one loop. We emphasize that the above result assumes all particles to be incoming. For outgoing particles one can simply replace the corresponding momentum p with −p. Furthermore, the LBK operatorD i and the function S(p i , p j , k) are defined in (3.11b) and (3.29b), respectively. A conceptual illustration of the factorisation formula (3.31) is shown in Figure 3. Contributions with hard and soft origin are depicted in green and orange, respectively. The first two diagrams on the r.h.s correspond to the hard sector given by (3.31b). The factorisation of (3.31c) into a universal soft function -connecting three external legs simultaneously -and the non-radiative matrix element is illustrated in the third diagram. Based on this, a naive extrapolation to higher orders in perturbation theory is possible by interpreting Figure 3 as an all-order statement. First of all, this would imply that the LBK operator yields the hard contribution also beyond one loop. More interestingly, however, it would significantly constrain the mixed hard-soft structure. At two loops, for example, the hard-soft region would be fixed through objects that already enter in (3.31). In particular, it would correspond to (3.31c) with M

Collinear factorisation at leading power
In this section we consider the leading-power behaviour of radiative amplitudes in the limit where the emitted photon (k) becomes collinear to a light on-shell fermion (p). This configuration is governed by the scale hierarchy where {s} denotes all other scales in the process. As mentioned in the introduction it was shown a long time ago that the leading-power contribution to the tree-level amplitude factorises in this limit into a process-independent splitting function multiplying the non-radiative amplitude with radiative kinematics [21][22][23][24]. It is the goal of this section to present the one-loop generalisation of this formula. We first start in Section 4.1 by reproducing the tree-level derivation from [21] and then discuss the one-loop extension in Section 4.2. The final factorisation formula is then given in (4.22) for initial-state radiation (ISR) and in (4.23) for final-state radiation (FSR).

Collinear factorisation at tree level
Contrary to the soft limit, care has to be taken in the collinear case when treating the gauge dependence of the emitted photon. Only axial gauge, where the sum over photon polarisations is given by does not mix up the power counting of individual diagrams. At tree level a convenient choice for the gauge vector r isk ≡ (E k , − k). At leading power in the collinear limit we therefore only need to consider diagrams where the photon is emitted from the collinear fermion leg. Restricting the discussion for the moment to ISR, we have We then write the fermion propagator in terms of quasi-real spinors [21] with energy It is then straightforward to derive the factorised result for the matrix element Figure 4: Interference terms that contribute at leading power in the limit where the emitted photon becomes collinear to the initial-state electron.
where the tree-level splitting function for ISR in its standard form [24] reads The analogous derivation for FSR yields Alternatively, J FSR (z, m) can be derived from J ISR (x, m) via the crossing relation p → −p, i.e. by replacing k · p → −k · p and x → z −1 .

Collinear factorisation at one loop
The quasi-real electron method from the previous section does not work anymore if loop corrections are taken into account due to non-factorisable diagrams where the photon is emitted from a loop. However, the method of regions can be applied in this case to disentangle universal contributions from collinear degrees of freedom from the process-dependent hard part. It is therefore possible to determine the splitting function based on a specific process. To this end we consider muon-electron scattering, i.e. e(p 1 )µ(q 1 ) → e(p 2 )µ(q 2 )γ(k) (4.10) and calculate the small-mass collinear limit of the one-loop corrections to the electron line. The scale hierarchy for ISR then reads Working at leading power and in axial gauge, we only need to take the four interference terms shown in Figure 4 into account. In particular, we have used the convenient choice r = p 2 for the gauge vector. This choice is allowed since p 2 2 = m 2 is small. Note that r = p 1 would not be permissible in this case since the small numerator k · p 1 in (4.2) would mix up the power counting.
The momenta p i = E i (1, n i β i ) of the energetic electrons can be decomposed into large and small components via the set of light-cone basis vectors {n i = (1, n i )/ √ 2,n i = (1, − n i )/ √ 2} which allows us to write any momentum as In the case of the energetic particle p i this takes the form of the desired decomposition where the scaling can be deduced from (4.14) Applying the light-cone decomposition to the external momenta in our process we find Based on the achieved disentanglement of scales it is possible to identify the contributing momentum regions. In order to do so it is helpful to use the formulation of the method of regions in the alpha parameter representation [27] automatised in the public code asy.m [28]. The following four regions are then found to contribute to the individual interference terms M i : hard: ℓ ∼ (1, 1, 1) 1 ∼ (1, 1, 1) 2 (4.16a) The terms that correct the incoming electron line, i.e. M 1 and M 2 , get only contributions from the p 1 -collinear region. Furthermore, at leading power the hard region only contributes to the factorisable diagram M 3 . Since we can apply the quasi-real electron method in this case it follows immediately that In addition to the hard region, all other three scalings contribute to M 3 . In the case of M 4 , on the other hand, only the p 1 -collinear and p 2 -ultra-collinear regions are present at leading power. As can be expected, the unphysical ultra-collinear region cancels between M 3 and M 4 . We are then left with the two collinear contributions that turn out to factorise according to Apart from the interference terms M i we also need to take into account mass and wave function renormalisation. All counterterms connected to the heavy particles (muon) enter in (4.17) in the renormalisation of the non-radiative massless one-loop matrix element M (1) n . The counterterms for the emitting electron, on the other hand, renormalise the one-loop splitting function J (1) ISR , while the ones for the other light particle (outgoing electron) contributes to Z (1) . The renormalised results for J (1) ISR and Z (1) are given in Appendix B. The factor Z (1) corresponds to the one-loop massification constant of [29]. Massification is a method to efficiently determine the leading mass effects of an amplitude solely based on the massless result. All mass terms that are not polynomially suppressed are recovered in this way. It is the universality of the aforementioned massification constant that makes such a reconstruction possible. The p 2 -collinear contribution therefore takes the leading-order mass effects of the outgoing electron into account. The one-loop splitting function J (1) ISR contains both the corresponding mass terms as well as leading-power corrections due to the collinear emission.
A non-trivial check for the above result is the behaviour in the soft limit where E coll corresponds to the eikonal factor in the collinear limit. We therefore get the expected form of the matrix element in the collinear-soft limit given by (4.20) As already mentioned previously the collinear contributions are expected to be process independent. Thus, one-loop diagrams for muon-electron scattering other than those shown in Figure 4 should not lead to such contributions. We have explicitly checked that this is the case due to a cancellation between diagram pairs that are related (up to a sign) through the crossing q 1 ↔ −q 2 . The only additional contribution is therefore the hard one originating from factorisable diagrams that trivially exhibit the factorising structure of (4.17). If, on the other hand, we take the muon to be light as well, i.e. M 2 ∼ m 2 ∼ λ 2 , there are two additional collinear contributions with exactly the same structure as for the outgoing electron consistent with the expectation based on massification. We are thus lead to the main result of this section that at one loop can be written through the factorisation formula where we have defined the all order quantities In (4.22) the product is over all external fermion lines j = i with a small mass m 2 j ∼ λ 2 . Furthermore, the same calculation with only minor modifications can also be applied to the case of FSR yielding the analogous formula A schematic illustration of these factorisation formulas is given in Figure 5. Furthermore, the corresponding expressions for J ISR , J FSR , and Z can be found in Appendix B. As expected, we find that the ISR and FSR splitting functions are related via crossing symmetry.
It is useful to compare our findings to the corresponding factorisation formula for massless fermions that can be extracted from the QCD results of [17]. Suppressing the separation into ISR and FSR, the massless collinear limit can be written as where the only process-independent contributionJ comes from the collinear fermion (p i ) and y ∈ {x, z}. The corresponding expressions at tree level and at one loop are given in Appendix B.
For massive particles, on the other hand, every light fermion contributes an additional factor in the factorisation formula taking into account the corresponding small-mass effects. This results in the more complex collinear structure of (4.22) and (4.23) than one would have naively expected based on the known QCD formula. Nevertheless, it turns out that there is a relation between the massive and massless splitting functions. In particular, we find where the massive splitting function reduces in the massless limit to the massless one plus singular corrections from massification. It is conceivable that the same relation will also hold beyond one loop. In this case, however, there will be non-vanishing soft contributions from massification [29]. In addition to being an interesting result in its own right, this represents a strong check for the validity of the results presented in this section and in Appendix B.

Validation
To demonstrate the correctness and applicability of equations (3.31), (4.22) and (4.23) we consider the soft and collinear limits in the process at one loop where k can become soft or collinear. This 2 → 4 process is a highly non-trivial test of our formalism since the full one-loop matrix element is rather involved and contains hexagon functions. We will compare our approximations to OpenLoops [30,31] running in quadruple precision mode [32]. This, while obviously slower, is remarkably stable and produces reliable results deep into the soft and collinear limits. The process (5.1) could also be considered to be the real-real-virtual contribution to the N 3 LO corrections to Bhabha scattering. Hence, implementing this matrix element in a way that remains sufficiently stable for collinear and soft emission would be essential for any future N 3 LO calculation.
In the following discussion we use a centre-of-mass energy of √ s = 10.583 GeV, tailored to the beam energy of the Belle II experiment.

Soft limit
Let us begin by considering the limit where one of the two photons becomes soft while the other photon remains hard, i.e. k → 0. Looking at (3.31), we have We emphasise again that this choice is not unique. It is therefore crucial to use the same definition both in the evaluation of the non-radiative matrix element as well as for the derivatives ∂/∂p i,µ in the LBK operator (3.11b). Since already the one-loop matrix element for ee → eeγ is rather complicated, we have implemented the corresponding derivatives numerically to a very high precision in Mathematica. Combining this with the soft contribution from (3.31c) then yields the complete subleading power approximation. The corresponding ξ −2 and ξ −1 terms can then be compared to OpenLoops as a function of the 'softness' 2E k / √ s. The result is shown in Figure 6 down to values of 10 −10 . It is clearly visible that including the ξ −1 (NLP) terms significantly improves the precision of the approximation. This behaviour clearly validates our one-loop generalisation of the LBK theorem presented in (3.31).

Collinear limit
Next, let us consider the collinear limit. Once the massless one-loop matrix element for the nonradiative process e − (p 1 )e + (p 2 ) → e − (p 3 )e + (p 4 )γ(p 5 ) is known, the application of the factorisation formulas (4.22) and (4.23) is rather straightforward. As an example we consider the case of the photon k becoming collinear to p 1 (ISR) or the case of it becoming collinear to p 3 (FSR). The cases of p 2 and p 4 are completely analogous.
The massified approximation [29] for the matrix element reads which is valid for the bulk of the phase space, i.e. assuming k is neither soft nor collinear. Note that the masses are only given indices so that the different Z can be better disentangled once k becomes collinear. Of course all m i are equal.
In the ISR limit we replace Z(m 1 ) with J ISR , reducing the number of particles in the massless matrix element M n+1 In complete analogy the FSR limit is given by M n+1 All that is left to do before we can compare to OpenLoops is to multiply out the terms in (5.5) and (5.6). The result of this comparison is shown in Figure 7 for ISR and FSR as a function of the 'collinearity' 1 − cos ∢(p i , k). To understand the observed convergence behaviour it is important to realise that the expansion is not performed in the collinearity but in k · p i , m 2 i → 0. The approximation thus only improves while k · p i gets smaller. At the point, however, where k · p i approximately satisfies (2.2), the limit saturates since m i is kept constant. This explains the kink at 10 −7 . We can therefore conclude that Figure 7 represents a strong validation of our factorisation formulas given in (4.22) and (4.23).

Conclusion and outlook
In this paper we have presented two novel findings about the universal structure of radiative QED amplitudes in the soft and in the collinear limit. We have extended the well-known LBK theorem to the one-loop level, relating the radiative with the non-radiative amplitude at subleading power in the soft limit. The additional loop effects are taken into account by supplementing the LBK theorem with a soft function that we have evaluated in a universal way. In addition, we have derived a factorisation formula at one loop that describes the leading-power collinear limit in the presence of small but non-vanishing fermion masses. Contrary to the analogous result for massless QCD we also get contributions from non-collinear light external fermions. These additional terms in the factorisation formula take into account the corresponding leading small-mass effects.
The approximation of real-virtual amplitudes with the soft limit at subleading power can be used to achieve a stable and efficient implementation of this contribution. In the case of Bhabha and Møller scattering this next-to-soft stabilisation enabled the first fully differential NNLO calculation. Even though the computation of the soft expansion was straightforward it turned out to be cumbersome. The extension of the LBK theorem presented in this paper significantly facilitates the application of the next-to-soft stabilisation method to other processes. This is particularly relevant in light of the MUonE experiment [33][34][35] where a theory prediction at the level of 10ppm is needed to achieve the targeted experimental precision. The minimal requirement for this is a fixed-order NNLO QED Monte Carlo matched to a parton shower. This has triggered a wide theory effort [36] where many partial results have been calculated in the past years. Two Monte Carlo codes are available that include the dominant electron-line corrections at NNLO [37,38]. In addition, the subset of the NNLO corrections with closed (and open) fermion loops are also known [39][40][41]. Further, electroweak effects [42] and possible contaminations from physics beyond the Standard Model [43][44][45] have been studied. Very recently a crucial step towards the full set of NNLO QED corrections was accomplished with the calculation of the two-loop amplitude with a non-zero muon mass [46]. After massification [29] of the massless electron this result can be included in a Monte Carlo code. The remaining bottleneck is therefore a numerically stable implementation of the real-virtual contribution. Next-to-soft stabilisation combined with the extension of the LBK theorem presented in this paper represents an elegant solution to this problem.
One could in principle also follow a similar approach in the case of the collinear limit. Contrary to the soft limit, however, the reliability of the approximation is limited by the scale hierarchy between the light fermion mass and the typical energy scale of the considered process. Leadingcollinear stabilisation could therefore only be used reliably in the case of high energies. On the other hand, also for low-energy processes the small electron mass leads to large peaks in the radiative amplitudes in collinear regions complicating the numerical integration over the phase space. The collinear factorisation formula presented in this paper could thus be used as basis for a subtraction scheme for these collinear pseudo-singularities similar to the NLO formalism developed in [24].
In principle, the same strategies applied in this paper can also be used to extend the nextto-soft and leading-collinear factorisation formulas beyond one loop. However, a more formal understanding of the presented results in terms of an effective field theory could facilitate this task. For the LBK theorem (3.31) this would entail a definition of the soft function (3.31c) in the framework of HQET. In the case of the collinear splitting function, on the other hand, a construction in the context of SCET would be needed.
Finally, we remark that this paper makes it possible to calculate contributions from radiative amplitudes based on the corresponding massless result. In order to do so one would use massification [29] for the bulk of the phase space and otherwise switch to the next-to-soft or leading-collinear approximations. While this is not strictly necessary at one loop it will be indispensable at two loop. The extension of the presented results for radiative two-loop amplitudes is therefore planned for the future. This in turn would represent a major step towards the fully differential calculation of the dominant N 3 LO corrections for muon-electron scattering.

A Soft integrals
In (3.29c) and (3.29d) we have defined the two integrals necessary to construct the soft contribution. These integrals are universal and are given here in d = 4 − 2ǫ dimensions with µ denoting the scale of dimensional regularisation. With s ij = 2p i · p j , s iγ = 2k · p i , and we have where v ij = 1 − 4m 2 i m 2 j /s 2 ij = (1 − χ)/(1 + χ). The Mathematica package HypExp [47,48] was used to expand the hypergeometric function in I 2 in terms of the harmonic polylogarithms (HPLs) of Remiddi and Vermaseren [49]. We then applied the related Mathematica code HPL [50,51] to simplify the resulting expression. Since 0 < v ij < 1 and thus 0 < χ < 1 all HPLs are manifestly real.

B Splitting functions
In the following we give the explicit expressions for all the quantities that enter the collinear factorisation formulas (4.22) (ISR) and (4.23) (FSR). The results are presented in a form that can be used in three major flavours of dimensional regularisation: the four-dimensional helicity scheme (fdh), 't Hooft-Veltman scheme (hv), and conventional dimensional regularisation (cdr) (see [52] and references therein for the definitions of these schemes). To this end, we keep the dimensionality of ǫ scalars, n ǫ , explicit in the poles but set it to zero in the finite parts. The regularisation-scheme dependence is therefore manifest as terms ∝ n ǫ . The corresponding results in hv and cdr can be obtained by setting n ǫ = 0. Inserting n ǫ = 2ǫ, on the other hand, retrieves the expressions in fdh. Furthermore, in the case of hv and fdh ǫ has to be set to zero in the tree-level splitting function.
We then define the invariant s kp = 2k · p where the photon momentum k is collinear to an initial-or final-state fermion p. The initial-state collinear splitting function can conveniently be written in terms of + H 2 (u) + H 1,1 (u) (−2u 2 x + 2u 2 + 2ux − 4u + 2). can be obtained from J ISR via the crossing relation p → −p. In particular, this implies x → z −1 and u → v −1 with The corresponding analytic continuation is unambiguously defined via the prescription s kp → s kp + iδ or equivalently u → u − iδ. We then find The imaginary part is given explicitly leaving all of the HPLs real for the physical region where 0 < z, u < 1. The massless version of the FSR splitting function entering (4.25) can be extracted from the spin-summed result of equations (II.10) and (II.11) in [17] by taking the QED limit. The corresponding expressions in the fdh scheme read