Locally finite two-loop amplitudes for off-shell multi-photon production in electron-positron annihilation

We study the singularity structure of two-loop QED amplitudes for the production of multiple off-shell photons in massless electron-positron annihilation and develop counterterms that remove their infrared and ultraviolet divergences point by point in the loop integrand. The remainders of the subtraction are integrable in four dimensions and can be computed in the future with numerical integration. The counterterms capture the divergences of the amplitudes and factorize in terms of the Born amplitude and the finite remainder of the one-loop amplitude. They consist of simple one- and two-loop integrals with at most three external momenta and can be integrated analytically in a simple manner with established methods. We uncover novel aspects of fully local IR factorization, where vertex and self-energy subdiagrams must be modified by new symmetrizations over loop momenta, in order to expose their tree-like tensor structures and hence factorization of IR singularities prior to loop integration. This work is a first step towards isolating locally the hard contributions of generic gauge theory amplitudes and rendering them integrable in exactly four dimensions with numerical methods.


Introduction
Perturbative Quantum Chromodynamics (QCD) has enjoyed rapid progress in the past two decades. Most remarkably, we have witnessed an "next-to-leading order (NLO) revolution" [1][2][3][4][5][6]: all important hard-scattering LHC processes can be computed, in an automated way, at next-to-leading-order in the strong-coupling expansion [7,8]. The automation of NLO computations has marked the culmination of years of challenging research which revealed elegant structures in gauge-theory amplitudes, loop integration and phase-space integration over real partonic radiation. Impressive progress on the same fronts (amplitudes, loop and phase-space integration) has been achieved beyond NLO. Many 2 → 2 processes are by now computable at next-to-next-to-leading-order (NNLO) [9][10][11][12][13][14][15], while hadron collider processes with simpler kinematics can be calculated at N 3 LO [16][17][18][19][20][21][22]. A first cross-section for a 2 → 3 process has been computed at NNLO in the approximation of a large number of colors [23]. Here we shall explore a factorization-based approach to that may be useful in the computation of multiloop amplitudes with multiple external lines, and give a first application at two loops in massless quantum electrodynamics. Currently, theoretical simulations and experimental measurements are in most cases of comparable accuracy. However, it is expected that with the collection of 15 times more data by the end of the LHC physics program, in the next two decades, and further advances in data analysis methods, the precision of many measurements will surpass the accuracy of theoretical predictions [24]. Theory will therefore be the limiting factor in the precise extraction of fundamental parameters, in testing the Standard Model of particle physics and in revealing potential signals of new phenomena. It will become important in the future to perform the many more NNLO and N 3 LO cross-section computations necessary to match the projected experimental uncertainties.
In their current form, however, the aforementioned methods are seriously challenged by the core of multi-loop amplitude computations needed for the purposes of the future precision physics program at the LHC. Their computational cost increases fast with complexity, measured either by the number of loops or the number of scattered particles. Next generation problems, such as three-loop amplitudes for 2 → 2 processes 1 beyond the massless case (for which we note recent spectacular breakthroughs [84][85][86]) and two-loop amplitudes for 2 → 3 processes with two or more mass parameters, appear to have a prohibitive algebraic complexity. Beyond the algebraic obstructions, on the analytic front, the "basis" of special functions (elliptic polylogarithms, etc.) which emerge routinely in such computations is not fully known. There is a realistic chance that the computational complexity cannot be tamed by following these well established and, up to now, very successful paths. New approaches must also be pursued. In this spirit, new ideas aiming at fundamental aspects of perturbative computations are emerging [87][88][89][90][91][92][93][94][95][96][97][98][99].
We believe that numerical integration directly in momentum space can provide a powerful solution to the computational complexity of multi-loop amplitudes. In momentum space, Feynman diagrams take the form of simple rational functions before loop integration. The evaluation of the multi-loop amplitude integrands is inexpensive, especially when modern recurrence and other generation methods [7,8,[100][101][102][103][104] are employed. The idea of integrating scattering amplitudes directly in momentum-space is appealing because at a loop order N L the number of integrations (≤ 4N L ) is independent of the number of scattered particles. This is very helpful in the effort to automate the evaluation of amplitudes for generic processes, since the number of integrations in other approaches depends strongly on the number of external particles. This is, for example, one main reason that while analytic approaches have been very successful in two-loop 2 → 2 amplitudes, the computation of 2 → 3 two-loop amplitudes is much more complicated [23,[105][106][107][108][109][110][111][112][113][114][115][116].
However, it is also true that direct momentum-space integration cannot be applied in a straightforward manner either. Feynman amplitudes are divergent in d = 4 space-time dimensions due to ultraviolet (UV) and infrared (IR) singularities. While UV singularities are canceled by renormalization, IR singularities cancel when all components of physical cross-sections (including phase-space integration over real radiation and convolution with parton densities/fragmentation functions) are assembled. In order to make this cancellation manifest and to obtain convergent representations of the physical remainder, the singularities need to be explicitly extracted from the integrand of multi-loop amplitudes. Similar extractions are of course required for phase-space integrands [13,15,19,22,.
Rendering amplitudes integrable with removal of IR and UV divergences followed by appropriate contour deformations has been pursued at one loop, in foundational works of Refs. [146][147][148][149][150][151][152][153]. At higher loop orders, the structure of IR and UV divergences is severely entangled and designing subtractions which remove simultaneously all singularities is challenging. In a recent work [154], it was shown that such singularities can be successfully removed in practice at higher loops by an amplitude-level "nested subtraction" approach [155][156][157][158]. There has also been related work on expressing amplitude integrands in a basis of integrals which makes the divergence structures manifest, known as the "local numerator" approach. This method has been applied in N = 4 super-Yang-Mills theory [159][160][161][162][163] and to less supersymmetric theories, including QCD [110,[164][165][166].
In this article, we demonstrate for the first time that two-loop amplitudes for the processes e − + e + → γ * 's can be freed locally, i.e., at the integrand level, from all infrared and ultraviolet divergences by a limited number of counterterm subtractions. The process of annihilation into off-shell photons is chosen for its relative simplicity, and because its IR structure is closely related to annihilation into massive electroweak bosons. The counterterms remove all types of soft, collinear and ultraviolet singularities. Infrared counterterms are summed into simple universal factors which are easy to integrate analytically, regardless of the number of final-state photons considered. Our approach is inspired by the picture of infrared factorization [155-157, 167-176, 176-183] which is established for amplitudes after integration of the loop momenta. Unlike Ref. [154], we will not attempt to construct counterterms for every individual diagram, but will construct counterterms for the sum of all diagrams, using a suitable choice of loop momentum routings to combine all diagrams into a single integrand. The benefit of summing over diagrams is that IR singularities become factorized, that is, proportional to lower-loop amplitudes, and the needed counterterms are very simple even when the number of external legs is very large.
The origin of gauge theory infrared divergences in terms of factorized and universal jet and soft functions has been demonstrated by a number of methods [156,183]. In this paper, we use a consequence of this factorization: processes with the same set of jet functions and soft infrared functions can differ only in their short-distance, "hard" functions. In this way, a known process can be used as a guide to isolate the hard functions in a more complex process. We shall use below the color-singlet annihilation form factor at two loops (here, in QED) as the known function, which we shall use to determine the much more complex annihilation cross section for n off-shell photons. The construction of counterterms below will follow the procedures implemented in proofs of factorization, and will rely on Ward identities that decouple unphysically polarized vectors from physical processes.
The efficacy of such a subtraction method depends on the logarithmic nature of infrared divergences in gauge theories. The power counting procedure that gives this important result has been known for a long time, and has been described in detail in the literature [184][185][186]. Key steps in these arguments depend on essential features of self energy and three-point vertex subdiagrams in massless renormalized perturbation theory. First, a renormalized fermion self energy diagram takes the form / p π F (p 2 /µ 2 ), with µ the renormalization scale. This is the inverse propagator times a function, π F , of the invariant mass that is no more than logarithmically divergent. In the same way, the renormalized vector self energy takes an analogous form, with a potentially logarithmically divergent function, π V , times the unique transverse tensor, (g µν q 2 − q µ q ν )π V (q 2 /µ 2 ). Finally, after integration and renormalization, a three-point subdiagram can depend only on its external momenta, and must obey the relevant Ward identities of the theory.
As discussed below, the tensor properties of two-and three-point functions that we have just described are a challenge to any local subtraction method, because they are only guaranteed to hold after the internal loop integrals have been carried out. At fixed values of internal loop momenta, self energy subdiagrams can produce power divergences in infrared limits. In addition, we will see that the internal loop momenta of vertex subdiagrams generically produce polarizations of their external photons that are in the directions of those loop momenta. We shall refer to these below as "loop polarizations". Loop polarizations, of course, are not necessarily in the directions of external momenta, and, as we shall see, interfere with the factorization, and hence universality, of infrared enhancements. All of these complications begin at two loops, and much of this article will be involved with showing how such obstacles can be overcome. Specifically, we will modify the loop integrand of certain diagrams by adding terms that vanish upon integration, to restore manifest logarithmic singularities that factorize into lower-loop expressions before integration is carried out.
After setting the notation in the following section, we discuss the underlying factorization properties of the class of processes under consideration, and point out the possibility of realizing this factorization at the level of the integrands. In Sec. 3 we construct a set of infrared and ultraviolet counterterms for the one-loop amplitude of multiple off-shell photon production in electron-positron annihilation. The infrared counterterms are manifestly proportional to the tree-level amplitude, thanks to Ward identities that reveal the factorized nature of the singularities. We continue to the construction of IR and UV counterterms for the two-loop amplitude, first for diagrams with fermion loops (Sec. 4), and then for diagrams with loops formed by additional photons (Sec. 5). In these constructions, we will modify certain one-loop self energy and vertex subgraphs to make factorization properties of IR divergences manifest before integration. UV subdivergence counterterms for one-loop subgraphs will be carefully constructed to obey Ward identities even away from the strict UV limit, to prevent spoiling the IR factorization properties of the remaining loop. We will present in Sec. 6 numerical checks to demonstrate that all IR and UV singularities are removed by our subtraction terms. The subtraction terms can themselves be analytically integrated by standard methods, and the detailed calculation will be presented in a forthcoming publication. We conclude with a summary of the procedures developed in this paper, and a brief discussion of the outlook for further progress.

Expansions and factorization for integrands
We consider the QED Lagrangian in the Feynman gauge, with N f massless lepton fields which carry an electric charge q f , normalized in units of the positron charge e, In addition to soft singularities, the theory also exhibits collinear singularities due to the masslessness of the fermionic fields. We will focus on processes with multiple photons produced in electron-positron annihilation, e − (p 1 ) + e + (p 2 ) → γ * (q 1 ) + . . . + γ * (q n ) , (2.2) where in parenthesis we denote the momenta of the external electron, positron and photons. We will consider here the case of off-shell photons, q 2 i = 0, focusing on the treatment of initial-state singularities. Although off-shell, these amplitudes are gauge invariant, and share many of the properties of amplitudes for the production of massive vector bosons. By restricting ourselves to off-shell photons only, we can concentrate on the IR divergences associated with initial states. Building on the formalism we develop below, we intend in future work to treat amplitudes with onshell photons and charged particles in the final state. 2

Expansions
The infrared structure of a probability amplitude for a process in this class, Eq. (2.2), is independent of the number of photons. We will therefore treat the final state in complete generality. For the initial state, we define the Mandelstam variable s ≡ (p 1 + p 2 ) 2 .
As an expansion series in the coupling constant, a generic n-photon production amplitude is given by where M (0) denotes the tree-level process and M (k) is the k-loop amplitude. At two loops, the n-photon production amplitude acquires an explicit dependence on the electric charges of all fermionic flavors due to closed fermion loops (as in Fig. 8  below). The two-loop amplitude can then be further decomposed as where q f e is the charge associated to the fermion species f . The term M (2) contains diagrams with no fermion loops. The amplitude coefficients M (2) c originate from Feynman diagrams with a closed fermion-loop, and c photons coupled to the loop. At two loops, as many as two of these photons may be virtual. Since we are working in QED, due to Furry's theorem, these coefficients vanish for odd values of c.
The one-loop, two-loop, etc., orders of the n-photon production amplitude are integrals over loop momenta We will adopt the above notation also for counterterms and finite amplitudes after subtraction; in all cases, M , with various superscripts and subscripts, will be used to denote integrated quantities while M will be used to denote unintegrated quantities. For the tree amplitude, there is no loop integration, so M (0) and M (0) will be used interchangeably. The tree-level and one-loop amplitudes, M (0) and M (1) , respectively, will be written down from diagrams in the Feynman gauge; the same is true for a starting expression of the two-loop M (2) , but a modified version of amplitude M (2) will be constructed, to make IR factorization properties more manifest without changing the integrated amplitude.
To denote counterterms, we will use the generic notation to represent a local (i.e. integrand-level) approximation of some amplitude M in a certain IR or UV limit r. For example, we will write R soft M (1) to denote the soft limit of the one-loop amplitude. We stress the inherent freedom and ambiguity left in this notation: R r M only needs to match the divergent part of M in the limit r, while the finite part is unspecified. The specific choice of counterterms will be a main topic of this paper. Now it's natural to use the notation to denote the result of integrating Eq. (2.6) over loop momenta. Due to the electron-positron initial state, the amplitude takes the form where M is a matrix in spinor space, and where the dependence of the amplitude on electron and positron spinors is shown explicitly. We will adopt a similar notation for various amplitude components, whereby M (X) ( M (X) ) denotes the integrated (unintegrated) component X of the amplitude with external spinors removed. For future reference, it will be useful to define the Dirac projectors, and The first of these projectors, also commonly used in SCET [187], will serve in our definition of infrared counterterms at one and two loops.

Factorization for integrals and integrands
As noted above, for the processes under study, all IR singularities are associated with virtual lines that are either collinear to the incoming electron or positron, or have vanishing momenta and are attached to one of the external lines of the pair or to fermion lines that are collinear to the incoming lines. All such IR behavior factorizes from the "hard scattering", in this case the production of off-shell photons. All such amplitudes enjoy the same pattern of factorization as the (Sudakov) annihilation form factor for massless fermions in QED. Let us briefly review the reasons for this result. Infrared singularities are generated by the loops of jet functions for the electron and for the positron, and by loops in a soft function, which describes the coupling of soft photons to those jets, and their scattering through virtual massless fermion loops. In fact, for our discussion, we will not need to construct these functions explicitly. We need only know that they factor from a hard-scattering, or short-distance function. The hard function, in turn, depends only on the momenta of outgoing photons and the incoming lines, and is completely insensitive to collinear evolution before the hard scattering, as well as to all soft photons. The factorization is therefore in terms of a product, rather than a convolution. Product factorization is characteristic of all wide-angle or large-momentum transfer scattering amplitudes in gauge theories [156,171].
There are many ways to express the product factorization described above. For us, it will be convenient to give it as a product in Dirac space, and in terms of integrands, following the notation of Eq. (2.5). (We suppress the subscript "full" in Eq. (2.5), but the following discussion applies independently to each of the terms in the expansion in the order of fermion loops, Eq. (2.4).) For the amplitude that describes the annihilation process e − (p 1 ) + e + (p 2 ) → γ * (q 1 ) + · · · + γ * (q n ), with all photons off-shell, we will represent the factorized amplitude as where in the first form we define the corresponding (annihilation) form factor as a matrix element. In the second expression, the matrix element is written in terms of a function F , which represents the all-order sum of form factor diagrams. Each diagram has a single, local, electron-positron vertex represented by P 1 T ({q j }) P 1 , at which the total momentum of the initial state, p 1 + p 2 , flows out of the diagram. The factor T is itself a Dirac matrix, representing a hard-scattering function that contains all dependence on the final state momenta {q j }. It appears between the projectors defined in Eq. (2.9), whose purpose we discuss at the close of this subsection. Like F , the matrix P 1 T P 1 will be computed from sums of loop diagrams. Its integrals, however, will be completely independent of the loop integrals of F , as described below.
The content of Eq. (2.12) is that all IR singularities arise from form factor integration, and are independent of the final state momenta. Correspondingly, the matrix P 1 T ({q j }) P 1 is IR finite in four dimensions. This relation holds to all orders in perturbation theory [156] for renormalized amplitudes. Equation (2.12) is also a sum over all orders in this fully factorized expression. To a given power of the coupling, say e 2a , a ≥ 1 in the notation of Eq. (2.3), M (a) ({q j }) is a function of a loop momenta, while F and T appear at all loop orders that add up to a, as (2.13) In this product factorization, the notation implies that the integrals over loops in the form factor F are completely decoupled from those in matrix T , whose integrands are represented by F and T , respectively. The form factor integrand for F (2) is illustrated Figure 1: Diagrams that contribute to F (2) for arbitrary order in the local vertex T .
diagrammatically in Fig. 1. Also included implicitly are the UV counterterms that are part of the perturbative expansion in these expressions for renormalized amplitudes.
Because of the separation of integrations, we can think of Eq. (2.12) as an effective theory expression for this class of amplitudes, where hard scales, of the order of external invariants, have been integrated out within T . From this point of view as well, the matrix P 1 T P 1 is a local vertex. The Dirac structure of this vertex depends on the nature of the hard scattering. For our case, with photon production only, T will be a Dirac vector, dressed by an IR finite coefficient function. For other electroweak processes, potentially including heavy lepton flavors in the hard scattering, axial vector, tensor and scalar matrices will all occur in general. In both this picture and the factorization formula, it is necessary to regulate UV divergences associated with the separation of scales, and to renormalize the operator represented here by the vertex.
We can now exploit the perturbative expansion of the factorized amplitude. To begin, in the sum over orders in Eq. (2.13), we have already separated the case b = a, for which the form factor at zeroth order is just the local vertex, evaluated between the external Dirac spinorsv(p 2 ) and u(p 1 ). We then find, using the properties of the projectors in Eq. (2.11), that at the level of integrands as well as integrals, we may identify (2.14) In the final definitions we recognize from Eq. (2.13) a term from which all IR subdivergences have been subtracted. The finiteness of this term at all loop orders is equivalent to the statement of factorization for this amplitude. Solving for M (a) finite in (2.13), we can interpret the additional terms as subtractions, what we will refer to as "infrared counterterms" below, which regulate all IR singularities in the integrand of the amplitude at this order.
Given (2.14), we can now identify, following Eq. (2.5) for amplitudes without external spinors, that This result enables us to solve iteratively for the finite remainder at order a + 1 in terms of those at lower order, At the lowest orders, we can immediately write down expressions that we will encounter again below, The advantage of this procedure is that all dependence on the multiplicity and masses of the produced photons is isolated in the functions M finite , which are infrared and ultraviolet finite after the integrals have been carried out. For a complex process, however, with high multiplicity in the final state, these integrals may be difficult to carry out analytically.
So far, our discussion is still entirely at the level of integrated amplitudes. As explained above, our goal is to formulate the M (a) finite so that they may be evaluated numerically.
With this in mind, we turn to the question of how to derive analogs of Eqs. (2.18) and (2.19) at the level of the integrands for these quantities, in forms that are at the same time locally integrable at singularities and manifestly convergent at infinity. To emphasize these two aspects, we write the analogs of (2.18) and (2.19) as In these expressions and below, M (a) will denote a sum of a-loop integrands without UV counterterms. Here the subscript "UV-finite" represents the relevant set of diagrams with ultraviolet counterterms subtracted. Schematically, we may write at a loops, for both M (a) and F (a) , where R i is an operator that inserts a counterterm appropriate for UV region i.
Terms not shown involve nestings of UV limits and UV limits where more than a single loop momentum diverges. We will use this notation, adapted to the specific cases being treated, below. In the same spirit, the form factor terms F in Eqs. (2.21) and (2.22), explicitly subtract the singular infrared behavior of the integrands of the original diagrams. We emphasize that these relations do not generally hold for amplitude and form factor integrands written down directly from the QED Lagrangian. Beyond one loop, it will be necessary to modify both the original amplitudes and the IR subtraction terms. These changes are necessary because factorization properties are expressions of the symmetries, especially gauge symmetry, of the theory, and such symmetries are not always manifest locally. We explain the challenges involved in finding a construction that obeys these relations in the following sections. We close this section with a brief discussion of counterterms, followed by an explanation of our use of the projection matrices, P 1 .
In our treatment, the ultraviolet counterterms in Eq. (2.23) will be written, not directly as a Laurent series in = 2 − d/2, but as a set of explicit integrals, whose integrands, once combined with the integrands of the original diagrams, produce convergent integrals overall. The counterterm integrals, which can be carried out independently as well, will depend on arbitrary masses, which serve to regularize them in the infrared, and in principle to match them to (massless) QED counterterms in any renormalization scheme. The set of counterterms for each diagram that contributes to M will include a subset that is in one-to-one correspondence to the original counterterms of the theory for these diagrams. We will find, however, that additional counterterms are necessary to make all integrals converge, even though these counterterms may add up to zero in the final integral. An example will be found below in light-by-light scattering diagrams.
Many of the same UV counterterms that appear in M will appear in the form factor integrals of F. We will also encounter UV counterterms associated with the renormalization of the vertex we have denoted by P 1 T P 1 .
In our discussion, we shall not consider contributions to the amplitude with the self-energies of the external lines, which we simply truncate and multiply by spinors or photon polarizations. Correspondingly, we do not carry out a multiplicative renormalization for our amplitudes. This process, however, is completely standard and unaltered by our construction. We will implement our construction explicitly for the one-and two-loop amplitudes of the full class of processes in Eq. (2.2).
Finally, before going on to the treatment of the one-loop amplitudes, Eq. (2.21), we will address the motivation for the use of the projectors P 1 . The example of M (1) finite shows, at the first nontrivial order, the essential role of these projections. Consider the factor F (1) , which is illustrated in Fig. 2. Because the integral in F (1) is fully decoupled from M (0) , the fermion propagators that are contracted with the local vertex in Fig. 2 are generically off-shell. This would imply that, without the projectors, the tree-level hard scattering function M (0) would lose gauge invariance. At tree level, this means that it would not be transverse to unphysical polarizations of the external photons.
To see the loss of gauge invariance explicitly, we consider a generic contribution to M (0) , with n external photons attached in some fixed order Π, with fixed polarizations, * (q i ). We isolate one of these, say * ν (q e ), which we will replace by its momentum q ν e . We will refer to such a polarization as "scalar" or "longitudinal", terms which are equivalent for massless vectors. For massive neutral scalars, q ν e is strictly speaking the scalar polarization.
Let us denote by M (0,Π) ν the sum of all insertions of the scalar-polarized photon of momentum q e into the fermion line from which the final-state photons emerge, keeping the order of the other photons fixed. This sum reduces to only two terms, determined by the QED Ward identity. The result is easily derived in complete generality by repeated use of the lowest order identity, which will be revisited in the next section. Here, we simply give the result, If this expression is inserted in Eq. (2.21) without the projectors P 1 , both terms will vanish only if both virtual fermion lines in Fig. 2 are on-shell. As a result, there will be gauge dependence in all regions except when the virtual photon carries zero momentum. In particular, in collinear regions for the photon's momentum, gauge dependence is induced. In contrast, with the projectors in place, this gaugedependence vanishes identically because / p 2 1 = / p 2 2 = 0. At the same time, in all singular regions, the projectors act as the identity for on-shell lines.
In summary, the result of our procedure will be expressions for the integrands of amplitudes, as in Eqs. (2.21) and (2.22), in which all infrared singularities are subtracted as form factor integrands (the F (L) in Eq. (2.13)), up to two loops. All process-dependent information is isolated in each M (a) finite , which can in principle be evaluated numerically. We anticipate that it will be possible to extend this procedure beyond two loops, and to other theories for a wide class of related processes.
In the following section, we turn to this formalism at one loop.

The one-loop amplitude
The one-loop amplitude M (1) for the process e − (p 1 ) + e + (p 2 ) → γ * (q 1 ) + . . . γ * (q n ) is a d-dimensional integral (with d = 4 − 2 in dimensional regularization) over the components of the loop momentum, as in Eq. (2.5). M (1) is singular in d = 4 spacetime dimensions, due to IR and UV divergences. In this section, we will describe how to remove these singularities systematically from the integrand, through the introduction of simple counterterms. Our approach will be by direct construction of the necessary subtractions. As we proceed, the general formula for the fullyintegrable and process-dependent result, M finite in Eq. (2.21), will emerge from the analysis. We emphasize that, even at one loop, we will encounter many individual diagrams that do not have the structure of Eq. (2.21).
Our aim is to identify a set of counterterms that are local in momentum space, and which add up to a single function M (1) singular (l) of the loop momentum that approximates the full one-loop integrand M (1) in all divergent limits. The integration will then be re-organized into two integrals, (3.1) The first integral on the right-hand side, which we will identify with M finite in Eq. (2.18), will be convergent in four space-time dimensions and amenable, in principle, to direct numerical integration. We will not discuss its practical computation here, but note that relevant methods have been developed in Ref. [188]. The second integral will contain all divergences in the dimensional regularization parameter . We will detail the construction of its integrand and its evaluation within dimensional regularization. We emphasize that the finiteness of the first term in Eq. (3.1) is a property of the sum over all one-loop diagrams, and the construction of counterterms will reflect properties of the sums of their integrands.
Generic IR and UV counterterms at one loop have already been derived in Refs. [149,150,189]. The counterterms that we present in this section are equivalent to the ones in the literature, in approximating all singular limits of one-loop amplitudes. Their derivation here serves to describe generic features of our methodology and as an introduction to the construction of the more complicated two-loop counterterms that will be presented for the first time, in later sections. Our derivation of local UV counterterms at one loop is designed, in particular, to respect Ward identities. The full set of counterterms will become important at two loops in constructing UV-subdivergence counterterms that do not spoil IR-factorization properties, which depend crucially on these Ward identities. As suggested in Sec. 2.2, we will see that the IR and UV counterterms necessary for a process with n photons match those of the e + e − annihilation form factor. In this way, we use our knowledge of the simplest process under consideration, e + + e − → γ * , to organize singular contributions to a much larger set of processes.

IR divergences of the one-loop amplitude
In this subsection we will present the subtraction of IR divergences from M (1) , reflected at the integrand level by the terms F (1) in Eq. (2.21). We will study their UV structure in the following subsection.
We may decompose the one-loop amplitude into four classes, each of which is characterized by the following diagrams: • the virtual photon is adjacent to both incoming fermions,

3)
• the virtual photon is adjacent to the incoming electron but not to the incoming positron, where the ellipses denote additional photon insertions into the "hard" subdiagram, consisting of all lines that are off-shell for l||p 1 , as well as permutations among the external photons, • the virtual photon is adjacent to the incoming positron but not to the incoming electron,

5)
• and, finally, the virtual photon is not adjacent to either of the incoming fermions We assign a loop momentum label l to all the diagrams, and the loop momentum flow is directed towards the vertex closest to the incoming electron. To compute diagrams analytically, it is often convenient to choose different momentum routings in different diagrams, as permitted by the shift-invariance of loop integrals. However, as will become apparent soon, such loop momentum label assignments can obscure important cancellations of collinear singularities in the sum of diagrams. A judicious assignment of the loop momentum for all diagrams is necessary in order to render these cancellations local in momentum space, and our momentum assignment achieves this purpose. In addition, the collinear limits will factorize in terms of the tree-amplitude, again locally at the level of the integrand. IR singularities are associated with configurations where the loop momentum l either vanishes in all four components, or is collinear to one of the two incoming lines. For the diagrams of M (1) , the complete list of such singularities is We refer to these points or lines as "pinch surfaces" [157,178]. To identify these regions, we recall that diagrammatic integrals are defined in the complex plane for each loop momentum component, and carried out along a contour that is specified by the "i " prescription for propagator denominators. In general, these contours can be deformed away from poles where propagator denominators vanish, unless poles coalesce from opposite sides of the contour, producing a "pinch". This happens from two poles of a single denominator, whenever the corresponding line carries zero momentum in all four components, or from pairs of denominators, when they share the light-like momentum of an external particle. The latter are collinear pinches, the former are soft. Note that when the virtual photon line of diagram M (1,A) carries zero momentum, both adjacent internal fermion lines are automatically on-shell, carrying their respective external momenta, so that in this case a soft pinch forces three lines to be on shell at once.
The set of collinear and soft singularities identified above do not exhaust the full set of pinch surfaces in the amplitudes that we treat here, but they are the only ones that produce infrared singularities. That is, they are the only pinch surfaces that are not integrable in four dimensions. The simplest examples of integrable pinch surfaces are points in loop momentum space where a fermion line carries zero momentum. All of our diagrams include such points within their loop momenta, but they are always integrable, because of the extra factor of momentum in the numerator of fermion propagators. We will not reproduce the power counting necessary to show these results here, but refer the reader to, for example, Refs. [157,178]. A recent and very general discussion has been given in Ref. [190].
An essential feature of the list of divergent pinch surfaces given above is that it is the same for any number of final-state photons in our basic set of processes, Eq. (2.2). As we will show, as indicated in Eq. (2.21), that this universality will be inherited by the IR counterterms necessary to factorize infrared behavior. Indeed, the full infrared behavior for every such amplitude is already present in the case n = 1 in Eq. (2.2), the annihilation form factor e + + e − → γ * .
We now seek a function M IR that approximates the one-loop integrand M (1) in all of the above IR-divergent limits for the general process with an arbitrary number of photons in the final state. The construction of such an approximation is complicated somewhat by overlaps of soft and collinear divergences. For example, the region l ∼ −x 1 p 1 , which yields a collinear singularity, overlaps for x 1 = 0 with the l ∼ 0 region, which yields a soft singularity. The entanglement of IR divergences takes more complicated patterns at higher loop orders. However, it has been shown that approximations of the integrands of arbitrary loop integrals can be constructed systematically with the method of nested subtractions [155][156][157]186]. In this approach, the diverse singular regions are ordered according to their extent within the integration volume. Then, the singularities are removed recursively, subtracting first the ones corresponding to the smallest volumes. Example applications of the nested subtractions method for one-and two-loop Feynman integrals can be found in Ref. [154]. In this article, we perform an explicit application of the method for physical amplitudes.
Following the nested subtraction approach, we first need to find an approximation of the integrand in the soft limit, which receives contributions only from the diagrams of the class M (1,A) . The sum of their integrands can be written compactly in the form where the numerator N (1) reads In the soft limit l µ → 0, we neglect l inside M (0) , and approximate M (1) with In this definition of the soft counterterm, we have left the denominators of the outermost electron and positron lines intact, in particular not dropping terms quadratic in l compared to p i · l, with i = 1, 2. In addition, we keep the exact l-dependence in the numerator except inside M (0) , where we set l = 0. We now turn our attention to the collinear region l||p 1 . Naively, we would have to analyze each diagram individually, but the situation is actually much simpler. In covariant gauges, collinear divergences appear when vector particles connect jet subdiagrams to the hard scattering. In individual diagrams, these contributions are non-factoring in their momentum dependence. Such collinear vector lines, however, always carry scalar polarizations (equivalent to longitudinal polarizations for massless particles), contracted with the hard scattering subdiagram. The Ward identities of the theory (in this case, QED) ensure that the sum of all such contributions factorizes at fixed momenta for the collinear photons [186]. This, by now standard, result enables us to introduce local counterterms for the sum of diagrams, avoiding the potentially complicated non-factoring dependence of each diagram individually.
Only the diagrams in M (1,A) and M (1,B 1 ) contribute to the collinear divergence in the region l p 1 . We write the sum of their integrands in the form (3.10) In this relation, M (0,A+B 1 ) µ (p 1 + l, p 2 , l; q 1 , . . . , q n ) represents the sum of tree diagrams from the set A + B 1 including all vertices at which the final state photons and an additional external photon of momentum l attach, but excluding the diagram where l is directly attached to p 1 . Momentum l flows out of each tree diagram and into the vertex adjacent to the electron line that is external to the full one-loop diagram, to be consistent with Eq. (3.7). The polarization vector of the exchanged photon (momentum l) as well as the electron and positron spinors are removed from M (0,A+B 1 ) µ . The momenta of the additional photon and the incoming electron are generically off-shell, but they approach the mass shell in the collinear limit l µ ≈ x 1 p µ 1 . Using η 1 as a reference vector, any q µ in the direction of p µ 1 , i.e., q µ = x q p µ 1 can be approximately rewritten as 11) or alternatively, Either of these approximations becomes accurate in the diagrams above when the photon of momentum l becomes collinear to p 1 , with q µ a linear combination of p µ 1 and l µ . While Eq. (3.11), with a denominator linear in l, suffices in order to demonstrate the factorization of collinear singularities (see, for example, Ref. [156]), in this article we use the equivalent version of Eq. (3.12) with a quadratic denominator.
Similarly, in the collinear limit l ≈ −x 1 p 1 , the metric tensor can be decomposed in terms of l µ and the generic vector η µ 1 as Staying always within the collinear approximation, where l 2 → 0, we can equivalently approximate the metric with (3.14) Using anti-commutation relations for Dirac matrices and the Dirac equation satisfied by the external spinor u(p 1 ), in the collinear limit l µ ≈ −x 1 p 1 , the rightmost factors in Eq. (3.10) can be re-written as, We now use these considerations to show how collinear singularities factorize on a point-by-point basis in our one-loop diagrams. By decomposing the metric tensor η µν in Eq. (3.10) according to Eq. (3.14), we notice that only one term gives a non-vanishing contribution in the collinear limit. Thus, in the limit l ≈ −xp 1 we can make the replacement η µν → 2η ν 1 l µ /[(l +η 1 ) 2 −η 2 1 ] without changing the behavior of M (1,A+B 1 ) in the region l||p 1 . This allows us to write We represent the right-hand side of Eq. (3.16) graphically by The triangle arrow at the lower end of photon lines denotes the following polarization approximation in the photon propagator: where η 1 is chosen to have a large rapidity separation from p 1 (i.e. it is not collinear to p 1 ). Here we are effectively considering a photon line whose polarization is proportional to its own momentum, a scalar polarization. Here and below we choose the triangle arrow notation to define a quadratic, rather than linear, denominator in the approximate polarization tensor. For the opposite collinear limit, we use a gray triangle to denote an analogous approximation, with a reference vector η 2 chosen to have a large rapidity separation from p 2 , The expression in Eq. (3.16) is now ready to be simplified with a partial fractioning identity, which is a manifestation of the QED Ward identity, This identity describes the effect of attaching a photon with a scalar polarization on a fermion line. We visualize the identity in Fig. 3. We also note the following identities for insertions of scalar photons in external fermion lines: which are depicted pictorially in Fig. 4. We now see that repeated applications of where M (0) is the lowest-order tree diagram, whose external vertices do not include the virtual photon. The cancellations are illustrated, as an example, for graphs in two-photon production in Fig. 5. With the action of the Ward identity, non- Figure 5: The application of Ward identities to factorize the l p 1 collinear divergence in the case of two-photon production. This is a graphical representation of Eq. (3.23).
factorizing collinear singularities cancel, and the one-loop collinear limit is expressed in terms of the lowest order, on-shell tree amplitude. The general case of an arbitrary number of photon emissions is proven analogously, and Eq. (3.23) holds generally. Interestingly, if we choose η 1 to be equal to −p 2 in Eq. (3.23), then the three denominators become identical to those of the soft-photon approximation given by the right-hand side of Eq. (3.9). After some simple Dirac algebra, it can be shown that the numerator of Eq. (3.23) also coincides with that of Eq. (3.9) in the collinear limit l p 1 . An analogous observation can be made for the opposite collinear limit l p 2 . Therefore, the right-hand side of Eq. (3.9), which is simply a form factor whose vertex is the tree amplitude (with l set to 0), may be used as a "global" IR counterterm that simultaneously cancels soft and collinear singularities of the one-loop amplitude. This is an expression of the universality of IR behavior in the complete set of amplitudes under consideration. As discussed in Sec. 2.2, we will place the tree amplitude between the Dirac projectors defined in Eq. (2.9). We arrive at the following global IR counterterm, which is equivalent to the soft photon approximation in Eq. (3.9) because P 1 acts as the identity in the soft, and collinear, regions. This is the usual one-loop form factor, which we will denote by F (1) , but with the QED vertex replaced by the truncated tree-level spinor matrix sandwiched between a pair of projectors, just as in Eq. (2.12). Therefore, the one-loop IR counterterm is equivalent to where we have recalled the corresponding diagrammatic notation of Fig. 2. The IR subtraction with the form factor counterterm of Eq. (3.25) suffices to remove locally infrared singularities for one-loop amplitudes where the final-state photons are on-shell as well. Even though diagrams such as in M (1,C) have collinear pinches when two fermion lines share the momentum of a light-like outgoing photon, numerator factors conspire to make them IR finite in physical amplitudes. In contrast, certain two-loop diagrams with on-shell photons are only finite after integration, in their original forms. An extension of our formalism at two loops for on-shell photons will be presented in future work.
We comment here on how the P 1 projector, as defined by Eq. (2.9), in Eq. (3.24) yields a tree-amplitude factor explicitly, in a regularization scheme where helicities are treated in four dimensions. Using the spin sum for either u(p 1 ) orv(p 2 ), the projector P 1 has the following two alternative representations, Without loss of generality, let us consider the case in which the incoming electron is left-handed, and the incoming positron is right handed. We replace u(p 1 ) by u L (p 1 ) and replacev(p 2 ) byv R (p 2 ) in Eq. (3.24). (Due to the chirality-preserving nature of interactions in massless QED, the amplitude vanishes when the incoming electron and positron have the same chirality.) These spinors satisfy, The IR counterterm R IR M (1) in Eq. (3.24) is rewritten as Since an odd number of gamma matrices preserves chirality, the above expression reduces to The square bracket in this equation is precisely the tree amplitude contracted with the external spinorsv R (p 2 ) and u L (p 1 ), so the counterterm of Eq. (3.24) is proportional to the tree amplitude even before integration. In summary, in this subsection we have shown that the following remainder of the one-loop amplitude leads to an IR-finite integration, with R IR M (1) given by the explicit subtraction of Eq. (3.31), or equivalently, in form factor notation by Eq. (3.25). Although this procedure provides an integral that is free of infrared singularities, to provide a numerically computable expression, we must also subtract UV-divergent behavior at the level of the integrand. We now turn to this procedure.

Ward identity-preserving Ultraviolet counterterms
The IR-finite one-loop integrand, Eq. (3.32), retains its non-convergent behavior in the UV limit. It is therefore not yet suitable for numerical evaluation. In order to remove the integral's UV-singularities through a local subtraction, we need to find an approximating function R UV M IR-finite that matches the singular behavior of the integrand in the UV-limit, with R IR M (1) given in Eq. (3.24). It is at this point that we will encounter UV divergences associated with the factorized amplitude, Eq. (2.12). At one loop in the process under study, UV divergences occur only in triangle diagrams and fermion self-energies, which we will denote by Γ (1),ν eeγ and Π (1) e , respectively. Following our convention for the routing of the loop momenta, these Green's functions have the integrands where p is in general off-shell and q can be the momentum of a real or (starting at two loops) virtual photon. The Ward identity for these Green's functions takes the form For the vertex diagram, we choose, as suggested above, a UV counterterm that is defined by its integrand, given in this case by which picks out the only term of the numerator that results in a logarithmic UV divergence, while regulating denominators with a common mass, M , which will play the role of a renormalization scale. We stress that in this paper, a "counterterm" is always defined at the integrand level, although the graphical notation above draws the counterterm as a local vertex. Our counterterms include all local counterterms in the conventional sense once, but only after, loop integration is carried out. For the fermion self energy, we must cancel both linear and logarithmic UV divergences, where the latter are linear in the external momentum. Again regulating denominators with a common mass, the resulting counterterm that we choose is In each case, the label (1) indicates that this is a one loop local counterterm. Notice that the above counterterms satisfy the Ward identity as well, for any values of l whether or not they are large compared to the scales of external momenta. Therefore, the renormalized remainders of the above Green's functions will automatically satisfy the Ward identity as well. Though not essential for the one-loop case, this will turn out to be a particularly useful property at two loops for obtaining factorized counterterms of collinear limits in the presence of UV subdivergences. Note that the Ward identity for the UV counterterms Eq. (3.39) is by no means guaranteed by the Ward identity for the original diagrams Eq. (3.36), since we always have the freedom of adjusting non-divergent contributions to the UV counterterms. For example, an alternative UV counterterm for the self energy diagram, by Nagy and Soper [189], is This alternative counterterm perfectly matches the UV-divergent behavior of Eq. (3.34) at both linearly and logarithmically divergent orders, but its finite parts would not preserve the Ward identities in conjunction with our UV vertex counterterm, Eq. (3.37). We have chosen non-divergent terms so that the UV counterterm for the self energy, Eq. (3.38), is aligned with our UV counterterm for the vertex, Eq. (3.37), for both finite and divergent parts, preserving the Ward identity locally. 3 We have now specified the local UV-counterterms of the one-loop amplitude, Finally, we need to find a UV approximation for the infrared counterterm, R IR M (1) in Eq. (3.24). We choose In the above, we have drawn a square vertex with an extra cross to denote the UV limit of the triangle diagram in Eq. (3.25) with a form factor vertex P 1 M (0) P 1 , analogous to how the diagram in Eq. (3.37) denotes the UV limit of the vertex diagram in Eq. (3.34), with a label "(1)" indicating that it is a 1-loop counterterm. Again the projector P 1 is always omitted in the diagrammatic notation. In terms of the relation, Eq. (2.21) that defines M (1) finite , this counterterm will produce a contribution to F We have now succeeded in finding an approximation of the one-loop amplitude integrand in all singular limits. We can then decompose the integrand as, (3.43) The first term on the right-hand side of Eq. (3.43), M finite (l), is integrable in four space-time dimensions. All singularities are now contained in Expressions for the counterterms on the right-hand side of Eq. (3.44) have been given by Eqs.
where in the second relation for F finite we have used Eq. (3.25). We note that the single IR subtraction, F (1) , removed soft and collinear divergences from a large class of diagrams, thus simplifying the application of the method. At one loop, our construction was carried out with unmodified integrands, specified directly by the QED Lagrangian. At two loops, we will find again that the set of IR subtractions is much simpler than the full set of diagrams. To achieve this simplification, however, it will be necessary to develop modified, but equivalent integrands.

Integration of singular counterterms
Before going on to the two-loop application of out method, we will complete the discussion of the one-loop case by integrating our one-loop counterterms in dimensional regularization. This is the lowest-order example of how the factorization procedure gives a set of closed expression with all singularities of the amplitude, both IR and UV, isolating all dependence on the final state in an integrable function, specified by Eq. (3.43).
Let us consider first the UV counterterms for the one-loop vertex and propagator, which appear in Eq. (3.41). The vertex counterterm integrates to, (3.46) or, diagrammatically, For the integral of the electron-propagator counterterm we find As expected, every vertex or propagator counterterm in Eq. (3.41) is proportional to the corresponding tree diagram. The proportionality factor is equal and opposite for the two types of counterterms, which is a direct consequence of the Ward identity (3.36). Therefore, in the integral, vertex and propagator counterterms cancel pairwise. Since there is one more vertex than propagators in any such graph, we deduce a very simple result for the integral of the UV counterterms. They are proportional to the tree-amplitude, We will return to the role of the Ward identity in the renormalized amplitude below. Equally simple is the integration of the UV counterterm of the form factor vertex, the integral of Eq. (3.42), Meanwhile, the integration of the IR counterterm yields Combining this result with the ultraviolet counterterm of Eq. (3.51) and the UV-IR combination counterterm of Eq. (3.52) to get the full singular contribution to the amplitude in Eq. (3.44), we find that the divergent part of the bare one-loop amplitude is 55) and the divergent part as → 0 is proportional to the tree-amplitude in agreement with previous expectations [170,171].

Two-loop diagrams with fermion loops
In the previous section, we exploited factorization to construct local infrared and ultraviolet counterterms for one-loop amplitudes. We are now ready to develop counterterms that achieve the same goal for two-loop amplitudes. The result will be a set of integrands that are locally integrable in exactly d = 4 dimensions. As mentioned above, to achieve integrability it will be necessary to develop a modified, but analytically equivalent, integrand for this set of diagrams with respect to expressions obtained by the direct application of Feynman rules. The finite remainders found from our modified integrands will satisfy Eq. (2.22), thus implementing factorization at the level of the integrand.
In this section, we begin the treatment of two-loop diagrams with the study of those diagrams with a fermion loop. As we shall see, the IR structure of these diagrams is closely related to that of the one-loop diagrams treated above. Nevertheless, we will encounter several obstacles to the construction of local infrared counterterms. In these diagrams, such issues will involve integrals over the fermion loops.
Once integrated, fermion loops enjoy features associated with the Ward identities of QED, such as the transversality of the vacuum polarization, and the UV-finiteness of light-by-light scattering. These are routinely invoked in treatments of the renormalization and IR structure of the full amplitude. Here, we seek to implement these features on a local basis in momentum space.
As we shall see, the integrands of diagrams with a one-loop photon vacuum polarization subgraph exhibit power singularities in addition to logarithmic singularities. They also give rise to anomalous polarizations ("loop polarizations" below) that spoil the local factorization in collinear limits that we found at one loop. These problems can be sidestepped at the cost of giving up manifest locality, if we integrate first over the vacuum polarization loop. This, however, is just what we want to avoid doing here. In addition, local collinear factorization is not manifest for individual Feynman diagrams with four-point or higher fermion loop subgraphs, which, while they do not contribute to collinear divergences after integration, exhibit collinear singularities locally.
We will solve these problems by deriving alternative forms for the amplitude integrands, for which the power-counting of singularities is canonical, and whose properties in collinear regions are locally consistent with factorization. For these modified amplitude integrands, we will be able to design local counterterms that cancel all infrared and ultraviolet singularities, point by point in the integration domain.

Diagrams with one-loop photon vacuum polarization subgraphs
Recalling the notation of Eq. (2.4), we start with M 2 , consisting of two-loop diagrams with a photon vacuum polarization subgraph, as in Fig. 6. The integrand  (3.6), by replacing the free photon propagator with its one-loop correction, We observe that in the single soft limit, l µ ∼ δ → 0, µ = 0, . . . 3, at fixed, non-lightlike k, the one-loop photon propagator integrand scales as which is more singular than at tree-level by two powers. This leads to a power-like singularity for the two-loop amplitude in the limit of soft l, and it is no longer possible to rely on the leading-power approximation to write down a factorized counterterm analogous to Eq. (3.9) of the one-loop case. Similarly enhanced are collinear singularities. For example, in the collinear limit l p 1 at generic values of the vacuum polarization loop momentum, the divergence in l is linear instead of logarithmic, simply because of the factor (l 2 ) −2 . At each such point, the polarization of the photon l may be proportional to k µ , which is arbitrary, and hence non-factorizing in general. Hence, at the local level, the two-loop integrand loses important features that we found at one loop. The soft and collinear non-factorizing and power-like singularities of the amplitude we have just identified, however, are only apparent at the local, fully unintegrated level. For example, after integrating over the fermion loop momentum k, a tree-like δ −2 scaling is restored for the propagator in the soft limit, where all components of l µ scale like δ. However, it is not necessary to integrate the fermion loop to reduce the degree of divergences in l to the logarithmic level. The vacuum polarization loop integral can be rewritten in a fully equivalent manner by, for example, Passarino-Veltman tensor reduction [191]. After this reduction, the integrand for the one-loop photon propagator reads, The result of the k loop integral for this expression is precisely equivalent to the original form, but in this expression the transversality of the photon is manifest. The above form can be replaced by an even simpler expression, in which logarithmic power counting is also regained. The same gauge-invariant amplitude will be found if the l µ l ν /l 2 term from the vacuum polarization is dropped in every diagram in which it appears. In practice, by repeated use of the same Ward identity as in Eq. (3.20), it is easy to show that when all diagrams are combined these terms produce a fully scaleless integral. Therefore, the longitudinal tensor gives a zero contribution to the amplitude integral and it can be dropped from the propagator. 4 After doing so, the tensorial form of the effective one-loop virtual propagator matches the one of the tree photon-propagator, (4.5) Using this result, it is simple to generate a suitable integrand for M 2 (k, l) from the corresponding one-loop amplitude integrand M (1) (l) as (4.6) This is our first example of a modified, alternative two-loop integrand, whose integral is equivalent to the original form, but whose infrared behavior is amenable to local subtraction and factorization as in Fig. 7. Infrared singularities in Eq. (4.6) are approximated by the counterterm where the "form factor" diagram R IR M (1) (l) = F (1) (l) is given by Eq. (3.24). The second identify defines the modified two-loop form factor subtraction for the contributions of this set of diagrams to the integrable two-loop integrand in Eq. (2.22). This modified form factor subtraction would emerge from exactly the same reasoning as above applied to the vacuum polarization in its diagrams. We now turn our attention to ultraviolet divergences, subtracting the UV subdivergences of M 2 as k → ∞ or l → ∞, as suggested in Eq. (2.23). Instead of giving separate but straightforward derivations of the UV counterterms, it is equivalent in this case simply to present the following full remainder after the subtraction, whose integrals are free of infrared and ultraviolet singularities, Notice that M 2, finite in Eq. (4.8) is not singular in the simultaneous k, l → ∞ limit, so it is not necessary to subtract an additional counterterm for the "double-UV" limit. In Eq. (4.9), the one-loop IR subtraction produces an integrand that vanishes like a power when l µ or l 2 vanishes, and thus controls logarithmic enhancements in momentum l µ from IR limits of the vacuum polarization integral over k. A straightforward expansion of the terms in Eq. (4.8), using (4.9), shows that it is precisely of the form anticipated in Eq. (2.22) in terms of our modified amplitude and form factor integrands.
In summary, we have decomposed the M

Other fermion-loop contributions
We will now treat the remaining classes of fermion-loop diagrams, with c > 2, where κ ≡ c − 2 final-state photons are emitted from the fermion loop, as in Fig. 8. The loop momentum l is always chosen to be the photon momentum flowing to a vertex adjacent to the p 1 electron, and the other loop momentum k is always chosen to be one of the fermion lines in the fermion box. 5 Due to Furry's theorem, amplitude  Some diagrams that contribute to M (2) κ+2 , κ > 0, illustrated by Fig. 9, develop collinear singularities when a virtual photon that is attached at one end to the incoming electron or positron and at the other to the fermion loop becomes parallel to the incoming fermion. The functional form of all such diagrams that contribute to the collinear singularity from l p 1 is (4.11) In the collinear limit l ≈ −xp 1 , the last few factors times the Dirac spinor u(p 1 ) become Since this is proportional to p 1µ , we can apply the approximation of Eqs. (3.12) and (3.18) with q replaced by l to rewrite Eq. (4.11) as The one-loop (κ+2)-photon amplitude M µν can be simplified further with the use of Ward identities. The two virtual photons that attach to these loops have momenta l and −l − κ i=1 q i flowing out of the loops. Of the two, the photon with momentum l is scalar-polarized. The action of the Ward identity is illustrated below for the case of two photons in the final state. We have, in the notation of Fig. 3 and Eq. (3.20), (4.14) We have used the shorthand notationl = −l − q 1 − q 2 , while "symmetric terms" denotes diagrams with other permutations of external photons q i . In the final equality, we find the difference of integrands for the same three-photon amplitude evaluated at two different values of the loop momentum, where M ν photons is defined in Eq. (4.14). This difference integrates to zero in dimensional regularization, which respects momentum-shift invariance. Therefore We again encounter a situation, however, in which integration over a loop momentum removes singularities which are present locally, and need to be removed for local integrability.
As for the vacuum polarization diagrams above, we will achieve the removal of local collinear singularities by a modification of the integrand, adding to it terms that vanish upon integration, as follows. Again, gauge invariance will make this possible.
Without changing the integrated value of the amplitude, we can exploit Eq. (4.16) and modify the integrand of Eq. (4.11) as, Here, we have modified the vertex adjacent to the incoming electron line, replacing γ µ byγ where η 1 is chosen so that it produces no new pinches from the extra denominator.
In the collinear limit, l = −xp 1 , we have An analogous modification will be made for diagrams which develop a singularity when a photon is collinear to the positron. For diagrams which can develop both types of collinear singularities, we perform these modifications simultaneously. We write, (4.20) The above expression integrates to the same value as withγ µ → γ µ , and is free of collinear singularities. We have chosen the reference vector for the p 2 -collinear divergence to be −η 2 , in accordance with Eq. (3.19). We note that in this case, there is no corresponding modification to a form factor subtraction term in Eq. (2.22). The modified contributions are now infrared-integrable on their own, without further subtraction.
In summary, we have modified the integrand of M k + q 1 + q 2 k k + q 1 k + q 1 + q 2 + q 3 q 2 q 1 q 4 q 3 Figure 10: The four-photon subdiagram with a fermion loop, which has unphysical UV divergences diagram by diagram.

UV counterterms for c = 4
In QED, the fermion loops in M c for c = 4 are UV finite after summing over diagrams that permute the photons connected to the loop, but only after the fermion loop integral has been carried out. The complete result must be finite, simply because there is no gauge-invariant local four-photon vertex. Nevertheless, individual diagrams for c = 4 (light-by-light scattering) are UV divergent in d = 4. To preserve local convergence, we must construct additional counterterms for these diagrams. Analogous counterterms were developed for QCD in Ref. [189]. Unlike QCD, however, for QED such counterterms have no counterpart in the renormalization of the theory.
The general form of a four-photon subdiagram is as illustrated in Fig. 10. In the UV limit, we write down the counterterm just as for the vacuum polarization in Eq. (4.8), by dropping all external momenta and IR-regulating denominators with a mass M , The counterterm in this expression corresponds to the specific subdiagram of Eq. (4.21). There are six such counterterms, corresponding to the number of non-cyclic permutations of indices, and hence to distinguishable attachments of the two virtual and two real photons to the loop. Our UV-convergent contribution to M 4, finite is just the sum of these six distinguishable diagrams, each minus its counterterm. Each Figure 11: The four-photon UV counterterm inserted into a one-loop diagram.
one of these combinations specifies a locally integrable two-loop integral, which is also UV convergent.
The corresponding contribution to M 4, singular is just the sum of the six counterterms, Eq. (4.22) plus permutations. Their sum can be combined into a single, local vertex embedded in the one-loop diagram in Fig. 11. Using tensor reduction, we obtain for this sum, which is necessarily symmetric in the indices µ 1 . . . µ 4 , an integrand that is proportional to = 2 − d/2, (4.23) The integral of loop momentum k has a pole, so that the sum of counterterms is finite in d = 4, (4.24) As noted above, this finite result does not correspond to a counterterm of QED renormalization. In fact, it is not by itself gauge invariant, simply because it does not vanish when index µ i is contracted with the corresponding photon's momentum. This is not a problem, however, because the integral of the diagram in Fig. 11 is locally IR finite and UV convergent. Figure 11 simply produces finite shifts in M 4, finite and M 4, singular , with the same magnitudes but with opposite signs. Its only effect is to enable us to express M 2 , the integrable remainder is found by subtracting singular terms that are closely related to one-loop IR and UV singularities. The remainder is given in convenient form by Eq. (4.8), which is equivalent to the subtraction specified in Eq. (2.22), in terms of the one-loop integrable remainder function. The sums of diagrams with c > 2, with final-state photons attached to the fermion loop, are finite, even though individual diagrams have collinear singularities, and in the case c = 4, ultraviolet divergences. For all such diagrams, modifications of the integrands themselves, as in Eq. (4.17), lead to expressions that are locally integrable on a diagram-by-diagram basis, and which taken together give the original integral. For the c = 4 light-by-light case, the introduction of a set of UV counterterms whose sum is finite provides a convergent integrand, while giving the same finite result as in the original form after integration. In the following section, we turn to the remaining two-loop diagrams, without fermion loops.

Two-loop photonic contributions
In this section, we will describe the construction of subtraction terms for the ultraviolet and infrared singularities in the M (2) integrand that consists of two-loop Feynman diagrams with no fermion loops. All of these diagrams will have two virtual photons. As for the fermion loop diagrams, the finite remainders found from our modified integrands will satisfy Eq. (2.22), again implementing factorization at the level of the integrand.
For the removal of infrared divergences, we will adopt the same form factor approach as we used for the one-loop amplitude, where the local integrand of a 2 → 1 form factor with an appropriately defined vertex served as the infrared counterterm for the amplitude of a generic 2 → n process. For this approach to succeed, collinear photons must factorize locally from the off-shell subdiagram in all singular limits, that is, prior to any integration.
As we shall see, the factorization of a collinear photon in certain two loop diagrams is less obvious than the one loop, precisely because of the presence of the extra loop momentum. The procedure we will describe below entails two steps: • first generating the two-loop integrand with an appropriate routing of loop momenta, • and then replacing some terms in the integrand with equivalent ones (i.e., which integrate to the same value) that are locally integrable in collinear singular limits.
The following three subsections detail and summarize these steps, which will bring the infrared singularities of this set of diagrams under control.
Once we have derived a locally integrable integrand for the two-loop amplitude, we will still need two sets of ultraviolet counterterms, just as we did at one loop. The elements of the first set are in one-to-one correspondence with the counterterms of the QED Lagrangian, but presented in integral form. They are formulated so that, when added to the full integrand, they give a fully convergent integral in all UV limits of the original amplitude. The second set of UV counterterms is required because form factor-based infrared counterterms induce additional UV divergences. We emphasize again that the latter counterterms are a familiar feature of many factorization techniques. Subsection 5.4 will deal with the construction of both sets of UV counterterms.
The outcome of the construction described below will be a reorganized expression for the full two-loop photonic contributions to the amplitude. It will be the sum of a set of "simple" but singular integrals, which can be carried out in dimensional regularization for any of the 2 → n processes in question, plus an integral that depends on all the details of the process, but which is numerically integrable in four dimensions. The number of IR and UV counterterms generated by this procedure will depend on the process, but will not be qualitatively different than for the original set of diagrams with their QED counterterms.

Generation of the two-loop integrand
We construct the two-loop integrand starting from the usual application of QED Feynman rules in the Feynman gauge. This construction is not unique, given the invariance of individual Feynman diagrams under shifts of the loop momenta. To combine all diagrams into a single integrand, we need to make an unambiguous choice of momentum labeling for every diagram. As at one loop, we aim to label loop momenta in a way which allows the application of Ward identities for a local cancellation of non-factorizing collinear singularities.
We label the momenta of the two virtual photons with l and k, both directed away from the vertex closest to the incoming positron. It is unimportant to specify which of the two photons carries momentum k and which carries momentum l, due to a symmetrization over momentum routings, which we will discuss shortly. This initial integrand is illustrated by, We will need to manipulate this form further, in order to produce a modified integrand, which we will denote as M (2) (k, l) , which will factorize in all singular collinear limits. Our manipulations will include symmetrizations of the loop momenta, and modifications of integrands that do not change integrated results. We detail these manipulations in the following subsections.

Global amplitude symmetrizations
As we will soon explain, we need to perform the four-fold symmetrization generated by where k ⊥ and l ⊥ are the components of the loop momenta k µ and l µ that are orthogonal to both the incoming electron momentum, p 1 , and the incoming positron momentum, p 2 . This symmetrization is applied globally, to all diagrams of the amplitude. The symmetrized integrand, , exhibits more complete factorization properties than the individual terms that make it up on the right-hand side. Here, the notation M To see the need for the k ↔ l symmetrization, let us examine, as an example, two-loop diagrams of the e − (p 1 )+e + (p 2 ) → γ * (q 1 )+γ * (q 2 ) amplitude that are singular in two different collinear limits. First, we consider a loop momentum k collinear to the positron momentum p 2 , with the second loop momentum l in the hard region. Analogous to the analysis of collinear singularities at one loop using Ward identities, a sum of diagrams which factorize in the collinear k p 2 singularity are, In the above, grey triangles denote the collinear approximation for photon propagators collinear to p 2 , defined in Eq. (3.19). The identity follows, as at one loop, by successive application of the lowest order Ward identity, Eq. (3.20). In Eq. (5.5), we have included only a typical class of diagrams where the momentum k is adjacent to the collinear leg p 2 and attaches to all propagators of the same one-loop subgraph. The dots denote omitted diagrams that can also be grouped into similar classes. In the collinear limit that we consider here, Ward identities (leading to factorization) can be applied to each class independently. Let's now look at a distinct limit, where the photon of momentum k becomes collinear to the incoming electron, k p 1 , In Eq. (5.6), we show the analogous set of Feynman diagrams that are needed as a sum for factorization, using white triangles to denote collinear approximations defined in Eq. (3.18). Once again, the collinear singularity occurs when the loop momentum k is adjacent to the external leg, which carries momentum p 1 this time. We now observe that the second Feynman diagram on the right-hand side of Eqs. (5.6) is actually also needed in Eq. (5.5), but with the momenta l and k exchanged. If we make a single assignment of loop momenta that is consistent with factorization of the k p 1 limit, then we will spoil the factorization of the k p 2 limit, and vice versa. In effect, the assignment of momenta can be thought of as an instruction for how to combine the integrands of the diagrams shown in these figures. With fixed assignments of line momenta as shown in the figures, the association of the points in momentum space where the Ward identity acts would not be automatic, but would require a change of variables. This conflict is resolved simply by symmetrizing over the two momentum assignments, and both collinear limits factorize for the combination M (c) (d) Figure 12: Diagrams requiring special treatment. We refer to the diagrams of the first line as type S, and of the second line as type V .
in Fig. 12, for the case when the self energy or vertex corrections are on the electron side. The case when they appear on the positron side is completely analogous and is omitted from the figure as well as much of the later discussions. We will see below the role of reflections in the transverse planes when we study collinear limits involving these diagrams.
We will refer to the diagrams with an electron or positron self energy correction outside the gray off-shell subdiagram, from which the final-state photons emerge, as type S, and the diagrams with a QED vertex correction outside the gray photon emission blob as type V. For concreteness, we will assign l to be the momentum of the virtual photon of the self energy or vertex sub-loop in type S and type V diagrams as in figure 12. This choice is simply a convention, since the integrand will be symmetrized according to Eq. (5.4) in the end. We call all the other diagrams regular diagrams, decomposing the two-loop integrand as a sum of contributions from these three classes, where M (2) sym denotes the full initial integrand with fixed momentum assignments and symmetrizations in Eq. (5.3) for each of the terms on the right. We will not alter the integrand of M (2) regular beyond these momentum symmetrizations, while M (2) type S and M (2) type V will require a somewhat more elaborate treatment.

Collinear regions in type S and type V diagrams
The diagrams in M (2) type S and M (2) type V of Fig. 12, until altered, do not exhibit factorization in all infrared singular limits locally. In this discussion, we demonstrate that local factorization can be achieved by overcoming three challenges, as summarized below and then elaborated upon.
1. Stronger than logarithmic divergences due to doubled fermion propagators, 2. apparent contributions of unphysical "loop polarizations" that are not compatible with factorization in collinear limits, 3. and, finally, scaleless integrals appearing as self energy corrections to external legs.
The first challenge, which we encountered for photon self-energies in the previous section, is apparent in diagrams of type S, where self energy corrections to a fermionic propagator adjacent to an external electron or positron are inserted. These diagrams possess single-soft k → 0 and collinear k p 1 singularities, whose strength is exacerbated by the presence of the self energy correction and the double occurrence of a denominator 6 . A straightforward power counting shows power divergences, instead of tamer logarithmic divergences. A naive subtraction of the infrared divergences would require counterterms for both the leading and next-to-leading singular terms, while the factorization of k p 1 collinear singularities with the application of Ward identities, schematically illustrated in Fig. 13, would apply only to the leading terms. Of course, if we perform the integral over l, the amputated self energy diagram gives a result proportional to / p 1 + / k which combines with the identical / p 1 + / k from the numerator of the adjacent fermion propagator to cancel the extra (p 1 +k) 2 denominator, but this would not achieve our goal of fully local factorization. Figure 13: Illustration of the desired factorization property of type S diagrams in the collinear k p 1 limit. Prior to integration, the propagator denominator (p 1 + k) 2 appears twice, giving rise to power divergences in the collinear limit. Appropriate modifications of the diagrams are needed to reduce the divergence to a logarithmic one, so that we can apply the collinear approximations, in particular the polarization approximations Eq. (3.18) with q = k indicated by the triangle arrow. The second challenge posed by "loop polarizations" emerges, for example, in the calculation of the collinear limit k → −xp 1 in type V diagrams. Let us focus on the vertex subgraph given by Eq. (3.34) multiplied with a common fermion propagator and a spinor factor, as they appear in type V diagrams, The tensor structure of the subgraph in the collinear limit is crucial for factorization. Analogous to the one-loop case in Eq. (3.15), the term with the vector p ν 1 in the numerator corresponds to a photon with a scalar polarization in the collinear limit. This term, when all Feynman diagrams of type V are summed and cancellations due to Ward identities are carried out, yields a factorized contribution in this limit. The l ν vector, on the other hand, corresponds to a photon with a random polarization, which we call "loop polarization", and it is problematic. Ward-identity cancellations are not present for these terms. Before going on, we note that both loop polarization terms in (5.8) have one fewer propagator than the original diagram, and appear in integrals characteristic of self energies, rather than the full vertex. This feature will be used to advantage below.
The third challenge concerns the factorization of k p 2 collinear singularities. Regular diagrams (i.e., not of type S or type V) factorize straightforwardly by Ward identities after using the collinear polarization approximation of Eq. (3.19) with q = k. But the sum of type S and type V diagrams, through Ward identities, becomes an expression which does not factorize into lower-loop amplitudes times an one-loop correction on the incoming positron (p 2 ) leg, as illustrated in Fig. 14. Of course, this problem goes away if we do not attempt to formulate factorization locally; after integrating over the sub-loop l, the self energy correction to the incoming electron line produces a factor of / p 1 , which annihilates the massless external spinor u(p 1 ), and the right-hand side of Fig. 14 vanishes. However, in order to be able to construct local counterterms for the amplitude, we need to exhibit factorization at the integrand level. It is therefore our aim to rewrite the type S and type V diagram expressions so that the sum in Fig. 14 vanishes locally in the limit k p 2 . Later we will construct the modified vertex integrand, Eq. (5.12) and modified self energy integrand, Eq. (5.18), which exactly satisfy the desired Ward identity, as illustrated in Fig. 15.

Strategy for factorizable vertex and self energy corrections
In the following, local factorization is made possible by specifying modified one-loop vertex and self energy integrands for subdiagrams that appear adjacent to the incoming lines. The modified integrands are designed to leave the integrated amplitudes unchanged in dimensional regularization. These modifications will apply to generic e − (p 1 )+e + (p 2 ) → γ * (q 1 )+. . .+γ * (q n ) amplitudes of an arbitrary photon multiplicity, n, in the final state, as well as to the two-loop 2 → 1 form factor. We will see that suitably modified integrands of the latter can serve as local IR counterterms of all other amplitudes in the class, providing a non-trivial two-loop generalization of the one-loop IR counterterm of Eq. (3.24).
Each integral in the type S and type V diagrams of Fig. 12 is characterized by an internal loop momentum for the vertex or self energy, which we will denote as l, and the momentum of the exchanged photon, which we denote by k. As noted above, the loop labels will be symmetrized, but the arguments here will not require that.
In these diagrams, the singular regions can be denoted by seven pairs (A, B). The first index, A, refers to the loop momentum k and takes the values A = 1, 2, S, H when k is collinear to p 1 , collinear to p 2 , soft or hard respectively. Here "hard" means the typical hard momentum transfer scale of the process, for example, the energy of final-state photons, which are taken to be all of the same order.  H).
To factorize the amplitude in the three single-collinear limits, where a collinear line attaches to a vertex carrying hard loop momentum, we will try to construct integrands that have four necessary and sufficient properties: • In region (1, H) (exchange photon k collinear to p 1 and loop momentum l hard and fixed) there should be no "loop polarizations" in singular terms. Leading singular behavior in this region should arise only from terms in which the photon with momentum k µ is scalar-polarized, i.e., with a polarization vector proportional to p µ 1 or k µ , on the index that is contracted with the hard scattering subdiagram. This polarization is captured by the polarization approximation introduced in Eq. (3.18).
• In the same region (1, H), power-like singularities coming from doubled propagators in type S diagrams, shown in Fig. 12(a)-(b), should be reduced to no stronger than logarithmic.
• In region (2, H) (k collinear to p 2 , l fixed), singular terms should obey the modified Ward identity for the scalar polarized photon of momentum k when it attaches to the type V vertex diagram on the p 1 line or is adjacent to the type S self energy diagram, as in Fig. 15.
• The Ward identity for a scalar-polarized photon carrying loop momentum l in the region (H, 1) (k hard and l collinear to p 1 ), illustrated in Fig. 16, is responsible for factorizing the l p 1 singularity into a product of a singular collinear function and a one-loop amplitude with loop momentum k. This identity is respected by the original Feynman diagram expressions and must remain valid after any modifications that we perform on type V diagrams. In practice this means that the modification of the type V diagram, which appears on the left-hand side of Fig. 16, should not change its singular behavior in the l p 1 limit. Figure 16: Result of the Ward identity in the l collinear to p 1 limit.

Modification of type V integrands
For type V diagrams in Fig. 12(c)-(d), i.e., QED vertex diagrams on an external fermion leg, with k the exchanged photon and l the loop momentum, the l sub-loop integrand can be written as Note that we have included the adjacent ( / p 1 + / k)/(p 1 + k) 2 fermion propagator in the expression for the vertex. The subscript in V 0 indicates that we will later modify this factor of the integrand.
After some Dirac algebra in d = 4 − 2 dimensions, using identities such as This factor includes all the dependence on loop momentum, l µ . The overall denominator (p 1 + k) 2 produces a pinch for k collinear to p 1 in combination with the denominator k 2 of the propagator for the photon emerging from the one-loop vertex. Similarly, a pinch singularity is produced for l collinear to p 1 due to the denominators l 2 and (l + p 1 ) 2 . Crucially, we observe that these two collinear singularities are disentangled. In particular, the terms in the first two lines of Eq. (5.11) that are singular for both k p 1 and l p 1 cancel in the l p 1 limit. Similarly, the final term of Eq. (5.11), proportional to / p 1 + / k, is suppressed in the k p 1 limit and is singular only when l p 1 . This natural separation of the (1, H) and (H, 1) collinear regions allows us to modify terms of the integrand which contribute to the k p 1 limit without affecting the singular behavior in the l p 1 region.
Based on this, we propose to replace V µ 0 in Eq. (5.11) by an equivalent form, in which loop polarizations proportional to l µ are eliminated from terms that are singular when the exchanged photon momentum k becomes collinear to p 1 . This is done by averaging the full expression on the first line of Eq. (5.11), containing the factor (p µ 1 + l µ )/[l 2 (p 1 + l) 2 ] , over a change of variables from l to −p 1 − l, and by averaging over an analogous change of variables in the 1/[(p 1 + l) 2 (p 1 + k + l) 2 ] term of the second line, exchanging p 1 + l with −(p 1 + k + l). Finally, anticipating our treatment of region (2, H), we explicitly write down the average k T ↔ −k T , l T ↔ −l T , as required by the global symmetrization of Eq. (5.3) applied to both regular diagrams and type S/type V diagrams, which we denote by a subscript T . Note that in either limits k p 1 or k p 2 , k T is close to zero and the symmetrization is essentially just l T ↔ −l T .
The resulting expression, only slightly more complicated than (5.11) is On the right, we can readily confirm the absence of factorization-breaking loop polarizations for singular contributions in the region (1, H). Explicit factors of l µ are suppressed in this region by the factor of k / acting on the Dirac spinor, or by an explicit factor of k 2 . We emphasize that although V µ mod in (5.12) is not equal to V µ 0 , its integral over l gives exactly the same result, so that it can be used in place of V µ 0 in the calculation of the amplitude.

Modification of type S integrands
The type S self energy integrand is treated in a manner analogous to the type V vertex. We seek a functional form that can be combined with the expressions for the vertex in the Ward identity relevant to the region (2, H), Fig. 15, when k is collinear to p 2 and l is hard. Here we shall simply present the form that will provide the correct combination, confirming that it is equivalent to the defining form of the self energy integrand, after integration.
We start with the integrand of the self energy subdiagram, which defines the original integrand, S µ 0 (l, p 1 , k). Again we have included the adjacent ( / p 1 + / k)/(p 1 + k) 2 fermion propagator in the expression for the self energy, just as we did for the vertex. In this standard form, the factor S 0 has a double pole at (p 1 + k) 2 = 0, times the truncated self energy, which, after integration over l, leaves a single pole in the combination.
To get started, we easily symmetrize the integral over a change of variable, in which the single pole is explicit, and proportional to the scalar self energy integral. Using the symmetries of the integral and the Dirac equation, there are many ways of finding an integrand for the l integral that is equivalent to the ones found from (5.13) and (5.14). We shall find an expression that preserves the Ward identity for the exchange vector k in regions (1, H) and (2, H). (We do not demand that the Ward identity be exact in other, finite regions.) In our next step, we provide a form equivalent to (5.14), which is a simple combination of two terms, related by a change of variables, To verify that the integral of S µ 1 equals that of S µ 0 , it is only necessary to change variables to l = l + k in the first term of (5.15). It is straightforward to check that this expression, together with the modified vertex expression, Eq. (5.12), obeys the modified Ward identity in region (2, H) as in Fig. 15.
As it stands, however, Eq. (5.15) is not quite what we want, because it introduces a new problem of loop polarizations, previously encountered in the unmodified type V diagrams, in region (1, H), where k is parallel to p 1 . These arise from the components of l µ in both the p 2 and transverse directions.
To solve this problem, we will analyze the loop momentum in the numerators of Eq. (5.15) in terms of components parallel and transverse to the incoming momenta p 1 and p 2 , We will perform a symmetrization only for the p 2 light-cone component in the two numerators of Eq. (5.15). We leave the p 1 light-cone component and the transverse components intact. This partial tensor reduction introduces modified versions of the loop momentum vector,l µ + andl µ − , defined bȳ in the two numerators. It is a simple exercise to show that replacing / l by / l + (or / l − ) in the first (or second) term in the square bracket of Eq. (5.15) does not change the result after integration over l. The resulting expression is, We are again using T to denote the transverse part of the global symmetrization Eq. (5.3). This is our final form for the type S self energy diagram of Fig. 12(a)-(b).
Eliminating the p 2 light-cone component of l removes its contributions to "loop polarizations". The transverse component of l, which also gives rise to loop polarizations, is eliminated in the k p 1 limit by averaging with a global reflection of the amplitude on the transverse plane in Eq. (5.3). This reflection is implicitly applied to every term of the amplitude integrand. Under this transverse reflection and averaging, the denominators in Eq. (5.18) are unchanged, but the / l T parts of the numerators cancel.
The transverse symmetrization is specifically designed to eliminate the "loop polarization" problem of the type S. Although this appears to be of a limited scope, the transverse symmetrization needs to be applied globally, i.e. to all diagrams as in Eq. (5.3), to preserve the interconnected web of Ward identity relations that enable factorization of divergences in a variety of other collinear limits.

Summary and validation of modified amplitude integrand in singlecollinear limits
We have now arrived at an alternative expression for the two-loop amplitude which solves the challenges of the original form. We first specified the assignments of loop momenta in all Feynman diagrams of the two-loop integrand M (2) 0 (k, l), derived in the Feynman gauge, in Eq. (5.1). Then we manipulated the integrand in the following manner.
1. For type V diagrams, i.e. vertex corrections adjacent to external fermion lines as illustrated in Fig. 12(c)-(d), we modified the integrand expression V µ 0 in Eq. (5.10) into V µ mod in Eq. (5.12) which is identical after integration.
2. For type S diagrams, i.e. self energy corrections adjacent to external fermion lines as illustrated in Fig. 12 The outcome is a modified integrand, which is a sum of three contributions, as in Eq. (5.7), now written as with the modified self energy and vertex integrands specified by Eqs. (5.18) and (5.12), respectively. The symmetrizations of Eq. (5.4) are understood. This will be the starting point for subtraction of infrared and ultraviolet divergences at the integrand level. The same modifications will be applied to two-loop form factor subtractions. Before we proceed with the subtraction of infrared divergences, we will verify that indeed the new integrand possesses the promised factorization properties. Specifically, below, we will show analytically that in the regions (H, 1) and (2, H), which were problematic originally, Ward identities and factorizations are now realized locally in momentum space. Factorization in the (1, H) region is self-evident by the absence of loop-polarizations in Eq. (5.12) and (5.18) and we do not discuss it any further.

Factorization in region (H, 1)
To test for the desired property in region (H, 1) (l collinear to p 1 ), we can replace the loop momentum vector l µ by −(1 − x)p µ 1 with 1 > x > 0 in Eq. (5.12), except in the singular denominators, l 2 and (p 1 + l) 2 . At fixed, off-shell k, we expect the photon carrying momentum l to satisfy a tree-level Ward identity, giving one term that factorizes l and k dependence, and another term in which they are linked, as shown in Fig. 16. The latter term will cancel against other diagrams that are singular when l is collinear to p 1 . Explicit calculation confirms that in the l → −(1 − x)p 1 limit, singular behavior follows the expected pattern.
As noted above, the vertex V µ mod , Eq. (5.12), is identical to the original form in this region, so we are just exhibiting how this works. Applying the Dirac equation and p 2 1 = 0 to the averaged vertex, Eq. (5.12), the singular part of the vertex diagram is the collection of terms proportional to the symmetric factor 1/l 2 (p 1 + l) 2 . These terms (which are even in l T ) are given, in detail, by where the second form results from the cancellation of all terms proportional to p µ 1 , and where we have observed that the Dirac equation implies p / 1 γ µ k / u(p 1 ) = −p / 1 k /γ µ u(p 1 ). A little further rewriting shows that this result equals the sum of terms associated with the Ward identity for a scalar-polarized vector particle, connecting a lowest-order jet with one of the diagrams of the one-loop hard subdiagram, This Ward identity is the expected one, and is represented in Fig. 16.

Region (2, H) and the k Ward identity
The final step is to verify that the Ward identity for the exchanged photon, Fig. 15, is respected by our modified integrand, Eq. (5.12), once it is added to the appropriatelymodified self energy, Eq. (5.18). Since we are interested in the behavior of the integral in the k collinear to p 2 region, (2, H), we can approximately set k 2 = 0, which simplifies the algebra considerably. Straightforward algebra and the Dirac equation then yield the relatively simple expression To this expression, we will add the self energy integrand. Neglecting terms that vanish at k 2 = 0, we find Given that ( / l − / l ± ) / p 2 = 0, (5.24) we verify, and thus the factorization is complete on a point-by-point basis, as illustrated in Fig.  15. This takes care of the (2, H) region. Of course, to get the point-by-point Ward identity, we had to use consistent routing of momenta in the vertex and self energy diagrams, and that is what our modifications were designed to achieve.

Two-loop "form factor" infrared counterterms
The amplitude integrands prepared in the previous subsections have the property that they factorize in all infrared limits. As a result of this factorization, infrared singular factors are common to all e − (p 1 ) + e + (p 2 ) → γ * (q 1 ) + . . . + γ * (q n ) amplitudes for arbitrary numbers of photons in the final state. In particular, the same singular factors are found in the 2 → 1 "form factor" amplitude, which will enable us to use this simpler process to organize the infrared structure of the entire class of processes. As discussed in Sec. 2, this universality is known to hold once loop momentum integrations are carried out. In this discussion, we have made the ingredients for factorization manifest at the level of integrands in the class of diagrams that we discuss.
To exploit this universality, we begin as in the one-loop case, and generate integrands for the form factor, now with modified vertex and self energy integrands, joining them to the tree level amplitudes for e − (p 1 ) + e + (p 2 ) → γ * (q 1 ) + . . . + γ * (q n ). We depict again the one-loop and two-loop form factor functionals, following the notations of Eqs. (2.14) and (3.25), as in Figs. 1 and 2, which we reproduce here as equations for ease of discussion, This result is the form we would like to generalize, including the treatment of its induced ultraviolet singularities.
We are now ready to discuss the subtraction of infrared singularities for the twoloop amplitude. Starting from the modified two-loop integrand M (2) constructed in Subsection 5.1, we will first subtract a "global" counterterm that simultaneously cancels all singularities in "double-IR" limits in which both loop momenta k and l are soft or collinear to incoming fermion lines, e.g., the double-soft limit with k µ , l µ much less than the hard scales, the soft-collinear limit with k µ soft, l p i (here i = 1 or 2), the limit of two collinear pairs with k p 1 , l p 2 , and the two-loop collinear limit with k, l p i . This global double-IR counterterm is precisely the set of diagrams given in Eq. (5.27), where we use the two-loop form factor as a local approximation of all double-infrared singularities. The subtraction gives in which we recognize the beginning of the pattern of the two-loop finite remainder, Eq. (2.22), so far with a subtraction only for the double-IR limits. The modification of type S and type V diagrams is understood in these expressions to match the behavior of the modified amplitude, as described above.
Next we must subtract a counterterm that cancels all remaining singularities in "single-IR" limits, where one of the two loop momenta becomes either soft (singlesoft) or collinear (single-collinear) to the incoming electron or positron, and one loop becomes hard while flowing through the short-distance subdiagram. (Type S and type V diagrams, in which a potentially hard loop is disconnected from the hard scattering are already taken care of.) For this, we are seeking an approximation of the two-loop integrand remainder after the subtraction of the double-IR singularities, Due to the local factorization properties of our amplitude integrand, the single-IR singularities of the two-loop amplitude factorize in terms of a one-loop amplitude with hard momentum and a universal factor, which is the same for any number of photons in the final-state and can be approximated with a one-loop form factor integrand. The factorization can be easily derived by using Ward identities, analogous to the one-loop case. For example, Fig. 17 illustrates the factorization of singularities in the single-collinear limit k p 1 in the sum of a class of diagrams. Factorization implies that the full contribution of regions in the amplitude where one loop is infrared and one hard is given by the product of the one-loop integrand in Eq. (2.21) and the one-loop form factor integrand, Similarly, we remove the single-IR limit from the two-loop form factor counterterm, Figure 17: Illustration of factorization of single-collinear singularities of the twoloop amplitude.
The above compact notation corresponds to the following momentum representation where M IR−finite is given by Eqs. (3.32) using (3.24), removing the external spinors. This provides the second IR subtraction in Eq. (2.22).
We recall from the discussion at the end of Sec. 2 that the P 1 projectors play an important role in Eq. (5.35) in preventing the introduction of spurious IR singularities. Let us examine their action here, in the momentum representation of the above counterterm, which reads IR-finite , which is defined without external spinors, is actually divergent. The divergence arises from terms of the amplitude that vanish in the Ward identity by the equation of motion for the external fermions. Such non-factoring terms cannot be canceled by a form factor subtraction. However, given that M (1) IR-finite is positioned in between a pair of projectors which play the role of the external spinors, in the sense that / p 1 P 1 = P 1/ p 2 = 0, these singularities are guaranteed to vanish. In the absence of the projectors, such singularities vanish only in regions where the other loop momentum l is soft or collinear to one of the incoming particles.
In conclusion, we have arrived at a compact preliminary result for the finite remainder, subtracted for all infrared limits of our modified two-loop amplitudes, Effectively, this is Eq. (2.22) before subtractions for UV behavior are carried out. Eq. (5.37) is valid locally, and the right-hand side is free of nonintegrable infrared behavior at all points in the integration domain of the loop momenta. It is worth pointing out that the number of IR counterterms in this expression is smaller than the number of diagrams that contribute to the amplitude.

Ultraviolet singularities of the photonic two-loop amplitude M (2)
The integrand M IR−finite defined by Eq. (5.37) can be thought of as an unrenormalized two-loop amplitude with its infrared singularities subtracted. However, the integral of this amplitude remainder is still singular in the ultraviolet limits k → ∞ and/or l → ∞. To complete our work, we need to remove these ultraviolet singularities from the integrand as well.
Our aim will be to construct ultraviolet counterterms that reflect the structure on the right-hand side of Eq. (5.37). As above, we denote them by the action of an operator R UV , so that, schematically, IR-finite P 1 . Similarly, the second term in Eq. (5.38) cancels the ultraviolet singularities in the two-loop form factor integrals of F (2) which serves as a counterterm of double infrared singularities, The last term in Eq. (5.38) must remove ultraviolet singularities from both the oneloop form factor integral of F (1) and its one-loop infrared-subtracted amplitude kernel IR-finite P 1 , We will produce the desired ultraviolet approximations in the equations above analogously to the quantum field theoretical procedure of ultraviolet renormalization, by introducing diagrammatic counterterms. These will have the effect of removing ultraviolet singularities from one-loop and two-loop propagator and vertex subgraphs whenever they appear in diagrams of the one-loop M (1) and two-loop M (2) amplitudes as well as the one-loop F (1) and two-loop F (2) form factor integrals. Our diagrammatic counterterms will be local in both momentum and coordinate space. Their integrands are functions of the loop momenta, which will combine with corresponding two-loop integrands to give convergent results. However, upon analytic integration of their loop momenta, their integrated values will also match the known poles of QED counterterms in the dimensional regulation parameter (common to any renormalization scheme). This local renormalization procedure completes the integrand decomposition of the two-loop bare amplitude into where all ultraviolet and infrared singularities reside in the first term, M singular , of the right-hand side and the second term M (2) finite is completely free of non-integrable singularities in d = 4 dimensions. The finite remainder is given by Eq. (2.22), reproduced here, where all terms on the right-hand side are convergent in the k → ∞ and/or l → ∞ limits while their sum is free of any infrared singularities. M singular is a sum of ultraviolet finite form factor one-loop and two-loop amplitudes whose analytic integration involves well-known integrals [192,193] and standard reduction methods. The integration of M (2) singular over the loop momenta in dimensional regularization will be presented in a separate publication.
In Section 3, we explained in detail the derivation of local ultraviolet counterterms for the diagrams in M (1) and F (1) . We use these counterterms at two loops as well, in order to cancel ultraviolet singularities from one-loop propagator and vertex subgraphs in diagrams of M (2) and F (2) . We will elaborate below essential infrared properties of their corresponding counterterm diagrams in configurations where one of the loop momenta is in the ultraviolet region and the other is in an infrared singular region. We will then extend the derivation of one-loop counterterms of one-loop self energy and vertex subgraphs for the modified integrands as they appear in type S and type V diagrams. Finally, we will discuss how to derive counterterms which respect infrared factorization, when both loop momenta tend to infinity.

One-loop ultraviolet counterterms in regular diagrams
In Section 3, we derived one-loop local ultraviolet counterterms for the one-loop QED vertex and electron self energy graphs, Eqs. (3.37) and (3.38), respectively. Our approach for their construction aimed to address the following problem. When a loop momentum, for example k, is taken to an infinite value, the second loop momentum l is unrestricted and can assume values which give rise to infrared singularities. For such loop momentum configurations, simultaneous subtractions are required for both ultraviolet and soft or collinear singularities. However, a naive introduction of ultraviolet counterterms could spoil the factorization of collinear singularities.
To address this complication, in Section 3 we introduced local ultraviolet counterterms for the one-loop electron propagator and the one-loop electron-electronphoton vertex that respect the QED Ward identity. Let us examine the function of the counterterms of Eq. (3.37) and (3.38), for diagrams in M (2) regular such as the ones illustrated in Fig. 18. The corresponding counterterm diagrams, in which the ultraviolet limit is approximated by Eqs. (3.37) and (3.38), exhibit the same cancellations of collinear singularities that occur in the original diagrams. The pattern of the Ward identity for the collinear regions of counterterm diagrams is illustrated in Fig. 19.
Therefore, our construction of ultraviolet counterterms ensures that collinear singularities associated with the non-UV loop remain in a factorized form. This factorized singularity can then be precisely canceled by the form factor amplitude F finite P 1 in Eq. (5.43). We remark that for this cancellation to occur, it is important that the local vertex of the form factor counterterm, P 1 M (1) finite P 1 , is UV regularized by the use of the same counterterms for one-loop self energy and vertex graphs as those we have used in M Figure 18: Illustration of ultraviolet counterterms, as defined in Eq. (3.37) and (3.38), for one-loop vertex and self energy subgraphs in regular diagrams (as opposed to type S and type V diagrams).

One-loop ultraviolet counterterms specific to type S and type V diagrams
To complete the derivation of one-loop ultraviolet counterterms for M (2) we have to discuss the ultraviolet limit of the integrand of one-loop vertex and self energy subgraphs in type V and type S diagrams, which we have cast in a modified form.
In the large l limit, the modified vertex expression of Eq. (5.12) is logarithmically divergent, so we can easily write down its UV subdivergence counterterm by taking the leading large-l limit of the expression, (5.44) Though the original type S electron self energy expression Eq. (5.13) is linearly divergent in the limit of large l, our equivalent form of Eq. (5.18) is only logarithmically divergent. This allows us to easily write down the UV counterterm Figure 19: Illustration of a Ward identity that must be satisfied by the UV counterterms for one-loop self energy and vertex subgraphs. Note that we are referring to counterterms at the integrand level, and the counterterm designations in the figure correspond to functions whose dependence on the loop momentum l is suppressed in the notation.
The symmetrization of the transverse components of the loop momentum ensures that those momenta remain UV finite when k is collinear to either p 1 or p 2 . Importantly, the above counterterms Eqs. (5.44) and (5.45) satisfy the necessary Ward identity for factorization (or rather, cancellation) of the singularity in the limit k p 2 , analogous to Eq. (5.9), (5.46) Because of this result, the UV counterterms do not re-introduce unfactorized collinear singularities when k p 2 . It is also easy to check that in the limit k p 1 , the above counterterms do not re-introduce the "loop polarization" problem. In these results the transverse symmetrizations in Eqs. (5.44) and (5.45) are necessary. Type S and type V diagrams are also present in the two-loop form factor amplitude F (2) , which we have used as an infrared counterterm on the right-hand side of Eq. (5.37). Consistently, we use the expressions of Eq. (5.44) and Eq. (5.45) as ultraviolet counterterms in the form factor F (2) as well. Figure 20: A sample graph that contributes to the second term on the right-hand side of Eq. (5.37), using the same notation as Eq. (3.25).

5.4.3
One-loop ultraviolet counterterms specific to the form factor F (1) and F (2) amplitudes In subsections 5.4.1 and 5.4.2, we constructed ultraviolet counterterms for one-loop vertex and self energy graphs in regular diagrams as well as type S and type V diagrams. Additional ultraviolet counterterms are needed for one-loop vertex subgraphs in the form factor F (1) and F (2) amplitudes on the right-hand side of Eq. (5.37) The ultraviolet singularity of F (1) (similar to the usual QED one-loop vertex) has already been encountered in Eq. (3.42). The same singularity appears in F (2) , for which a sample diagram is shown of Fig. 20. The large-l limit in Fig. 20 is singular and the singularity can be also subtracted by a counterterm that is identical to the square bracket of Eq. (3.42).
(5.47) The subtraction of the one-loop singularity as l becomes large, by using the counterterms of Eq.  A further counterterm for the one-loop singularity as k becomes large is required and it is similar to Eq. (3.42). The ultraviolet approximation of k → ∞, following the l → ∞ subtraction, reads

Double ultraviolet counterterms
After discussing the subtraction of ultraviolet singularities in single limits where only one of the loop momenta approaches infinity, we turn our attention to double ultraviolet singularities in the region where both loop momenta k, l become infinitely large. As in the standard renormalization procedure at two loops, the subtraction of such singularities is performed after one-loop singularities have already been subtracted. This order of subtractions for singularities in ultraviolet regions is necessary to avoid double counting, and leads to counterterms for the singularities in the double-ultraviolet region which produce integrands that are also convergent in single-ultraviolet regions. Double-ultraviolet singularities are found in two-loop one-particle-irreducible (1PI) propagator and vertex graphs or subgraphs in the diagrams of M (2) and the two-loop form factor amplitude F (2) of Eq. (5.27). Counterterms for these singularities can be constructed straightforwardly, with the following steps.
1. We first perform a series expansion of the integrands of the singular 1PI graphs/subgraphs in the limit that both loop momenta become large uniformly, and truncate the series to the order that contributes to logarithmic UV singularities. As a result, we obtain tadpole-type integrands whose only denominators are k 2 , l 2 (k ± l) 2 , each raised to some non-negative integer power.
2. We then insert masses into all the denominators to eliminate infrared singularities without changing the UV behavior, as was done at one loop in Eqs. (3.37) and (3.38). Some of the denominators are already massive by this step, because of the previous subtraction of one-loop UV sub-divergences, and these denominators are kept massive.
The above procedure leads to counterterms that preserve the Ward identities of vertex and propagator graphs at two loops. This is a natural consequence of the Taylor expansion around the ultraviolet limit being truncated consistently at the same order and the introduction of infrared mass regulators being performed in the same manner. The benefit of Ward identity-preserving ultraviolet counterterms is two-fold. First, the integration of the counterterms leads to simple results, making manifest cancellations of double-ultraviolet counterterms of diagrams with two-loop vertex subgraphs against diagrams with two-loop self energy subgraphs. Second, it permits an extension of our study to three-loop amplitudes in the future, since our double-ultraviolet counterterms will not destroy the factorization of collinear singularities associated with a third loop.
Taken together with the modified integrands of the two-loop amplitudes, and the results of Section 4, the UV counterterms described above enable us to construct the fully-subtracted two-loop amplitude given in Eq. (5.43) or equivalently (2.22). The number of UV counterterms is similar to the number for the original diagrams, because the number of form factor subtractions is fixed, independent of the number of final-state photons in the process. In the following section, we will verify the integrability of our construction for Eq. (2.22).
6 Numerical checks of cancellation of singularities for offshell multi-photon production In this section, we numerically check the cancellation of all IR / UV singularities at the integrand level in e − (p 1 ) + e + (p 2 ) → γ * (q 1 ) + γ * (q 2 ). This is an example from the class of processes for which our formalism applies. Our analysis includes all the terms that make up M (2) finite in the full expression of Eq. (2.22). We generate the Feynman diagrams using the program QGRAF [194]. Then we use custom Mathematica code to insert the Feynman rules, perform loop momentum symmetrizations, and modify certain diagram expressions as specified in the preceding text. Among the full set of two-loop diagrams, we may summarize the following modifications.
• The finite part of diagrams M (2) 2 , with vacuum polarizations, including both IR and UV counterterms, is generated as specified by Eq. (4.8).
• Fermion loop diagrams, M (2) c with c ≥ 4, are evaluated with modified vertices adjacent to the incoming electron and/or positron lines, given by Eq. (4.18).
• The diagrams for M (2) 4 , which include light-by-light scattering, have UV subtractions, as described in Sec. 4.3, which add to a finite result, but which produce convergence in otherwise divergent individual diagrams.
• For diagrams without fermion loops (photonic contributions), we use the modified two-loop integrands for type V and type S subdiagrams, specified by Eqs. (5.12) and (5.18), respectively.
• The type S and type V modified integrands are also employed for the form factor subtractions, as described in Sec. 5.3.
• Ultraviolet counterterms that ensure convergence for two-loop photonic amplitudes and their IR subtractions are constructed as in Sec. 5.4.
We check that the integrand has the correct scaling behavior to be convergent in all possible infrared and ultraviolet limits at arbitrary phase space points. For illustration, we present results for one such point, A basis of vectors transverse to the p 1 -p 2 plane is defined as The numerical checks explained below have been repeated with various random choices of external / internal momenta, polarization vectors, and spinors, for us to be confident of the cancellation of IR / UV singularities.

Photonic contributions
The M (0) tree amplitude for e − +e + → γ * +γ * has a "t−channel" and a "u−channel" Feynman diagram which are symmetric under the exchange of the two final state photons. Of course, this exchange symmetry is valid for the amplitude at all orders in perturbation theory. At any given number of loops, we can identify a diagram as belonging to the t or the u channel by eliminating all virtual photons from it and mapping it correspondingly onto one of the two tree Feynman diagrams.
Our form factor subtractions act on the t and u channels independently of each other, rendering their corresponding amplitude contributions finite separately. To verify the cancellation of singularities it is sufficient to consider only one of the channels. In the following, we present results only for one of the two orderings for the final-state photons.
To test convergence in a singular region, we assign numerical values for the exponents ω + k , ω − k , ω T k , ω + l , ω − l , ω T l which parameterize it appropriately [184][185][186]. For example, in the single-soft limit where k is soft and l is hard, we set Since all computations are over rational numbers (recall that all momenta, polarizations, and spinors are chosen to have rational numerical components), the above numerical value is computed with no rounding error, and the final results (involving very large numerators and denominators) have been converted into floating point numbers here only for readability. Therefore we are confident that the expansion coefficients are strictly zero at O(1/δ 4 ). Since the integration measure scales as d 4 k d 4 l ∼ δ 4 · δ 0 = δ 4 , the integral scales as δ and is convergent in the single-soft limit.
We repeat the same procedure for all IR / UV limits and we present the scaling of the t-channel contribution as δ → 0 in Table 1.   (2) in various limits. The scaling of the integration measure has been taken into account, and the amplitude is convergent if it scales as any positive power of δ. The loop momenta k and l are parameterized as Eq. (6.5), with the scaling exponents given in this table. Since the integrals always scale as δ raised to a positive power, there is no divergence in any of the IR / UV limits.
For brevity, only a set of independent limits is presented. For example, due to the global symmetrization k ↔ l, the k p 1 limit is identical to the l p 2 limit. Also, since p 1 and p 2 are essentially treated in the same way, the k p 2 limit is not shown separately.
We note that the subtracted integrand converges faster than expected in certain IR limits. For example, in the double-soft limit, as shown in the 3rd row of Table  1, the four-dimensional finite integral is suppressed by δ 2 instead of δ. This will be beneficial for high-precision numerical integration, and we leave it to future work to fully understand the unexpected fast convergence properties.
We have also checked that when a loop momentum tends to infinity (the UV limit) while the second loop momentum becomes soft or collinear, the integrand also scales in a way that guarantees finiteness of the integral. This is simply a consequence of the Ward-identity preserving nature of our UV subtraction terms, which guarantees the factorization of infrared divergences with respect to one of the loop momenta point by point in the space of the value of the other loop momentum. Table 2 shows the relevant IR / UV scaling behavior of M  in various limits. The loop momenta k and l are parameterized as Eq. (6.5), with the scaling exponents given in this table. The scaling of the integration measure has been taken into account. Since the integrals always scale as δ raised to a positive power, there is no divergence in any of the IR / UV limits. Table 3 shows the relevant IR / UV scaling behavior of M (2) 4 which receives contributions from six Feynman diagrams with fermion box subdiagrams. We consider the sum of all these diagrams in our numerical check without deleting those that can be obtained from others by crossing.  in various limits. The loop momenta k and l are parameterized as Eq. (6.5), with the scaling exponents given in this table. The scaling of the integration measure has been taken into account. Since the integrals always scale as δ raised to a positive power, there is no divergence in any of the IR / UV limits.

Conclusions
A striking property of gauge theory scattering amplitudes is the factorization of their infrared singularities. In this article, we have taken a first step towards understanding this structure better at a practical level, and developing a subtraction algorithm for the evaluation of multi-loop amplitudes numerically in exactly four dimensions. We have worked in an Abelian gauge theory, which we may think of as the perturbative expansion of massless quantum electrodynamics, and in Feynman gauge.
In loop amplitudes, infrared factorization is generally not manifest point by point in momentum space, or on a diagram by diagram basis. Factorization becomes manifest for the amplitude when gauge symmetry is exploited, as Ward identities eliminate nonfactorizable collinear singularities in combinations of Feynman diagrams. At the integral level, one can make simplifications that permit these group cancellations. First, "hard" contributions can be explicitly integrated out and UV renormalized. This operation renders them finite and projects them to point-like tensor structures, which in turn generate only scalar-polarizations for virtual photons in collinear regions. Second, approximations of Feynman integrals in the relevant divergent limits can be manipulated independently of all other limits. It is therefore permitted and helpful to shift loop momenta differently for each one of these limits, which is seemingly essential in order to demonstrate the cancellations of non-factorizing terms. These manipulations of the integrand are helpful in demonstrating general properties of amplitudes, such as the factorization of universal infrared functions from processdependent hard scattering functions. They do not, however, provide an algorithmic scheme for the numerical calculation of those functions.
An aim of this work has been to investigate whether infrared factorization can be implemented in an alternative way, by exclusively making use of procedures that are strictly local, i.e. valid not only for the integral but also for the integrand. By committing to implementing factorization at the integrand level, the two tools in the previous paragraph cannot be employed. That is, we cannot rely on the results of sub-integrals, or independent shifts of loop momenta. Indeed, integrating out the hard regions leads to a loss of locality, which prevents the construction of local infrared counterterms. Locality is also incompatible with performing independent shifts of the momenta in diverse singular regions of the amplitude. Rather, a global choice of loop momentum routings should achieve the cancellation of non-factorizing divergences in all possible singular limits simultaneously when diagrams are combined without further shifts in momentum.
In this article, we have demonstrated that this goal can be achieved through two loops for a large class of processes in a realistic gauge theory. In our case study, we considered a class of QED amplitudes in Feynman gauge, for the production of multiple off-shell photons in electron-positron annihilation.
For these amplitudes, we generated a two-loop integrand in a form in which factorization is manifest for all infrared singular limits. Our integrand is an average over four copies of a usual Feynman diagrammatic integrand with symmetric, judiciously chosen, loop momenta routings.
Terms in the integrand originating from one-loop vertex and self energy subgraphs require renormalization and contribute to the amplitude at "hard" scales, comparable to the hard momentum transfer and beyond. Using the symmetries of these integrations, however, we know that once they are fully carried out in dimensional regularization they produce only unphysical scalar polarizations on collinear photons that connect the vertex correction to the hard scattering. For such photons, QED Ward identities ensure that non-factorizing contributions cancel at the pointlike level in a sum over diagrams. This reasoning, however, does not yet tell us how to actually carry out the relevant integrals numerically. Indeed, for generic values of the vertex or self energy loop momentum, these subdiagrams generate "loop" polarizations on their external photons, in arbitrary directions, as opposed to scalar polarizations in the collinear direction.
In this article, we developed alternative local representations of these subgraphs, different than, although of course derived from, those obtained by a direct application of Feynman rules, with properties required for factorization in all collinear and soft regions. Our novel representation of the amplitude integrand achieves the goal of making factorization in all singular limits manifest locally. It allows only scalar polarizations for those collinear virtual photons that give rise to singular contributions. This ensures a simultaneous cancellation of non-factorizing IR divergent contributions in all possible singular limits and preserves, of course, the integrated value of the original amplitude.
With local factorization manifest, we were able to remove infrared singularities with local, universal counterterms for all e + + e − → γ * 1 + . . . + γ * n amplitudes, irrespective of the photon multiplicity, n, in the final state. Specifically, we used a simple form factor amplitude through two loops as a "generating functional" for our infrared counterterms. Each form factor integrand is defined by a finite local vertex, either the tree amplitude or the finite remainder of the one-loop amplitude, between Dirac projection matrices. The expression for the infrared finite remainder of the two loop amplitudes in Eq. (5.37) is an elegant outcome of factorization, and it constitutes a main result of this article. We anticipate that this form will generalize to higher orders. We remark that Eq. (5.37) holds generally for final-states where photons are replaced by other (W, Z, H) massive bosons. It is also valid for the Abelian contributions to QCD amplitudes for the production of generic colorless final-states with off-shell photons or other massive bosons in quark anti-quark annihilation.
To remove ultraviolet singularities, we introduced local counterterms for one-and two-loop Green's functions with two, three and, in the case of diagrams with closed fermion loops, four legs. A delicate issue in constructing our ultraviolet counterterms has been the preservation of factorization in singular regions where one of the loop momenta lies in the ultraviolet regime while the second loop momentum is collinear to an incoming particle. To achieve this, we constructed integrands for ultraviolet counterterms for self energy and vertex graphs with consistent expansions around the ultraviolet limit, which preserve Ward identities and, consequently, local collinear factorization.
With the factorization-based method of this article we were able to separate explicitly, for the first time in an amplitude with arbitrary numbers of external particles, and in a realistic gauge theory, the finite remainder of two-loop amplitude integrands from all of their infrared and ultraviolet singularities. As a test, we computed numerically the degree of divergence for all pinch singularities anticipated in the two-loop e + + e − → γ * + γ * amplitude and demonstrated that the subtracted remainder is indeed integrable, and can in principle be evaluated with numerical integration in exactly d = 4 dimensions.
We have carried out the analytic integration of the form factor infrared counterterms and of the ultraviolet counterterms with a straightforward application of tensor reduction techniques to well known master integrals. Our analytic results for the divergent part of the amplitude will be presented in a forthcoming publication.
Once local subtractions are in place, numerical integration of the finite remainders is possible, but presents additional challenges. An important step in this direction has already been taken in Ref [97] and in Ref. [188] where amplitudes, from the class of processes with equivalent subtractions for pinch singularities as described here, were evaluated efficiently at one loop. Ref. [97,188] presents an algorithm in which, after a first analytic integration of the energy components of the loop momenta, a novel deformation of the integration path for the remaining integrals away from non-pinched singularities is constructed. The efficiency of this numerical integration method was tested in Ref. [97,188] with examples of finite integrals with many external legs and beyond two loops. The application of the same technique for the integration of the finite remainders of the two-loop amplitudes presented in this article is an exciting possibility.
In this article, we formulated a novel subtraction method based on local infrared factorization at two loops for a class of QED amplitudes for the production of off shell photons. We are looking forward to extending our method to generic QCD amplitudes for multi-particle scattering processes at two loops and higher orders in future work.