title: Locally finite two-loop amplitudes for off-shell multi-photon production in electron-positron annihilation
authors: Charalampos Anastasiou Rayan Haindl George Sterman Zhou Yang Mao Zeng
text: 
  =1
   
  
  
  
  
  
  
  
  
  arrows,shapes
  
  
  
  figures/
  
  
      type Stype VType SType Va]Charalampos Anastasiou,a]Rayan Haindl,b]George Sterman,a]Zhou Yang,a]Mao Zeng[a]Institute for Theoretical Physics, ETH Zurich, 8093 Zürich, Switzerland[b]C.N. Yang Institute for Theoretical Physics and Department of Physics and Astronomy, Stony Brook University, Stony Brook, NY 11794,
      USAbabis@phys.ethz.chhaindlr@phys.ethz.chgeorge.sterman@stonybrook.eduyangzho@student.ethz.chmzeng@phys.ethz.ch 
      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.
      YITP-SB-2020-25Locally finite two-loop amplitudes for off-shell multi-photon production in electron-positron annihilation
          [
          December 30, 2023
      ==========================================================================================================
      
      
      
      
      
      § 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” <cit.>:
      all important  hard-scattering LHC processes can be
      computed, in an automated way, at next-to-leading-order in the strong-coupling expansion <cit.>.
      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) <cit.>, while hadron collider processes with simpler kinematics can be
      calculated at N^3LO <cit.>. 
      A first cross-section for a 2 → 3 process has been computed at NNLO in the approximation of a large number of colors <cit.>. 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 <cit.>.  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^3LO cross-section computations necessary to match the projected experimental uncertainties.
      
      
      We believe that the biggest stumbling block
      towards this goal is the computation of new two-loop and three-loop amplitudes. 
      This may be surprising given that we have seen an enormous progress in the development of
      methods for multi-loop amplitude computations in recent years. On the analytic front,
      powerful algebraic Feynman diagram simplifications with reductions to “master integrals” have been
      devised <cit.>. The master integrals themselves have been objects of an
      intense field of research and their evaluation is nowadays made with
      powerful mathematical methods which exploit newly discovered 
      algebraic properties of polylogarithmic  functions and iterated
      integrals <cit.>. 
      On the numerical side, methods such as 
      solving master integral differential equations numerically <cit.>, sector decomposition <cit.>,
      Mellin-Barnes integration <cit.>, and direct parametric integration <cit.> using (quasi-) finite bases <cit.> have proven to be very powerful. 
      
      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 [In supersymmetric gauge
        theories, scattering amplitudes  of even higher external state
        multiplicities and at four loops have been computed
        recently, e.g. in <cit.>.]
      beyond the massless case (for which we note recent spectacular breakthroughs <cit.>) 
      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 <cit.>.
      
      
      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 <cit.> 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 (≤ 4 N_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 <cit.>. 
      
      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 <cit.>. 
      
      
      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. <cit.>. 
      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 <cit.>, it was shown
      that such singularities can be successfully removed in practice at higher loops by an amplitude-level “nested subtraction” approach <cit.>. 
      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 𝒩=4 super-Yang-Mills theory <cit.> and to less supersymmetric theories, including QCD <cit.>.
      
      In this article, we demonstrate for the first time that two-loop amplitudes for the 
      processes e^-+ e^+ →γ^* 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
       <cit.>
      which is established for amplitudes after integration of the loop momenta.
      Unlike Ref. <cit.>, 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 <cit.>.   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 <cit.>.   
      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. <ref> 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. <ref>), and then for diagrams with loops formed by additional photons (Sec. <ref>). 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. <ref> 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,
      
          L = ∑_f=1^N_f i ψ̅_f ( ∂ -
            i q_f e  A
          ) ψ_f - 1/4F_μν F^μν -1/2 ξ(∂· A )^2,     ξ =1.
      
      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)  ,
      
      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_i^2 ≠ 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 on-shell photons and charged particles in the final state.[Beginning at two loops, even for outgoing on-shell photons we will encounter “transient" final-state singularities, which disappear in the fully integrated amplitude, but which require special treatment to achieve local integrability.]
      
       §.§ Expansions
      
      
      The infrared structure of a probability amplitude for a process in
      this class, Eq. (<ref>), 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
      
          M     = e^n [  
          M^(0) + e^2 M^(1) + e^4 M^(2)_ full + …],
      
      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. <ref> below).
      The two-loop amplitude can  then be  further decomposed as 
      
          M^(2)_ full =  M^(2) + ∑_c=2^n+2( ∑_f q_f^c ) M_c^(2)  ,
      
      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_c^(2) 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 
      
          M^(1) = ∫d^dl/(2π)^d ℳ^(1)(l),    M^(2) = ∫d^d k/(2π)^dd^d l/(2π)^d ℳ^(2)(k,l)  ,…
      
      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  ℳ will be used to
      denote unintegrated quantities.   For the tree amplitude, there is no 
      loop integration, so M^(0) and ℳ^(0) will be used 
      interchangeably. The tree-level and one-loop amplitudes, ℳ^(0) and ℳ^(1), respectively, will be written down from diagrams in the Feynman gauge; the same is true for a starting expression of the two-loop ℳ^(2), but a modified version of amplitude ℳ^(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
      
          _rℳ
      
      to represent a local (i.e. integrand-level) approximation of some amplitude ℳ in a certain IR or UV limit r. For example, we will write _softℳ^(1) to denote the soft limit of the one-loop amplitude. We stress the inherent freedom and ambiguity left in this notation: _rℳ only needs to match the divergent part of ℳ 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
      
          _r M
      
      to denote the result of integrating Eq. (<ref>) over loop momenta.
      
      
      Due to the electron-positron initial state, the amplitude takes the
      form 
      
          ℳ=v(p_2) ℳ (p_1, p_2, {q_j} ) u(p_1)  ,  { q_j} =  q_1, … , q_n  ,
      
      where ℳ  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) (ℳ^(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,
      
          𝐏_1    ≡ p_1  p_2/2p_1 · p_2,  
            𝐏_2 ≡1 -𝐏_1
      
      which satisfy
      
          𝐏_1^2 = 𝐏_1,    𝐏_2^2 = 𝐏_2,
      
      and 
      
          𝐏_1  u(p_1) = u(p_1),   v̅(p_2) 𝐏_1 = v̅(p_2), 
            𝐏_2  u(p_1) = 0,   v̅(p_2) 𝐏_2 = 0   .
      
      The first of these projectors, also commonly used in SCET <cit.>,
      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 <cit.>. 
      
      
      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. (<ref>).   (We suppress the subscript “full" in Eq. (<ref>), but the
      following discussion applies independently to each of the terms in the expansion in the order of fermion loops, Eq. (<ref>).)    
      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
      
          M ( p_1,p_2,{q_j})     =    ⟨ 0 | ψ̅(0)  [ 𝐏_1 T(p_1,p_2, {q_j}) 𝐏_1  ]  ψ(0)  | p_1,p_2 ⟩
             =    F  [ 𝐏_1 T( {q_j}) 𝐏_1  ]  
            ,
      
      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  𝐏_1 T( {q_j}) 𝐏_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. (<ref>), whose purpose we discuss at the close of this subsection.
      Like F, the matrix 𝐏_1 T 𝐏_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. (<ref>) is that all IR singularities arise from form factor integration, and are
      independent of the final state momenta.  Correspondingly, the matrix 𝐏_1 T( {q_j}) 𝐏_1 is IR finite in four dimensions.   This relation holds to all orders in perturbation theory <cit.> for renormalized amplitudes.
      
      Equation (<ref>) 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. (<ref>), 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
      
          M^(a)( {q_j} )    =   ∫_k_1 … a  M^(a)( {q_j}, k_1, … k_a ) 
          
             =    ∑_b=0^a-1∫_k_b+1… a  F^(a-b) [ k_b+1, … k_a  ; 𝐏_1 ∫_k_1… b  T^(b)( {q_j}, k_1, … k_b) 𝐏_1  ] 
          
                    +  F^(0) [ 𝐏_1 ∫_k_1… a T^(a)( {q_j}, k_1, … k_a  ) 𝐏_1  ]  .
      
      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 diagrammatically in Fig. <ref>.   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. (<ref>) 
      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 𝐏_1 T 𝐏_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. (<ref>), 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 spinors
      v̅(p_2) and u(p_1).   We then find, using the properties of the projectors in Eq. (<ref>), that at the level of integrands as well as integrals, we may identify
      
          F^(0) [ 𝐏_1  ∫_k_1… a  T^(a)( {q_j}, k_1, … k_a  ) 𝐏_1  ]        
                  =v̅(p_2)  [ 𝐏_1  ∫_k_1… a  T^(a)( {q_j}, k_1, … k_a  ) 𝐏_1  ]   u(p_1)
          
                  = v̅(p_2) ∫_k_1… a  T^(a)( {q_j}, k_1, … k_a  )   u(p_1)
          
                 ≡ ∫_k_1… a  M^(a)_ finite ({q_j}, k_1, … k_a  )
          
                 ≡ M^(a)_ finite ({q_j} )  .
      
      In the final definitions we recognize from Eq. (<ref>)  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 (<ref>), 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 (<ref>), we can now identify, following Eq. (<ref>) for amplitudes without external spinors, that
       
          T^(a) = M^(a)_ finite  .
      
       This result enables us to solve iteratively for the finite remainder at order a+1 in terms of those at lower order,
       
          M^(a+1)_ finite = M^(a+1) - ∑_b=0^aF^(a+1-b) [ 𝐏_1 M_ finite^(b)𝐏_1  ]  .
      
      At the lowest orders, we can immediately write down expressions that we will encounter again below,
      
          M^(0)_ finite    =    M^(0)  ,
          
          M^(1)_ finite    =    M^(1) - F^(1) [ 𝐏_1 M^(0)𝐏_1  ]   ,
          
          M^(2)_ finite    =    M^(2) - F^(1) [ 𝐏_1 M^(1)_ finite𝐏_1  ] 
          - F^(2) [ 𝐏_1 M^(0)𝐏_1  ]  .
      
      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. (<ref>) and (<ref>) 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 (<ref>) and (<ref>) as
      
          M^(0)_ finite    =     M^(0)  ,
          
           M^(1)_ finite    =     M^(1)_ UV-finite -  F_ UV-finite^(1) [ 𝐏_1  M^(0)𝐏_1  ]   ,
          
           M^(2)_ finite    =     M^(2)_ UV-finite -  F^(1)_ UV-finite [ 𝐏_1  M^(1)_ finite𝐏_1  ] 
          -  F^(2)_ UV-finite [ 𝐏_1  M^(0)𝐏_1  ]  .
      
      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),
      
          M^(a)_ UV-finite    =     M^(a) - ∑_ UV regions i _i  M^(a) + ⋯  ,
          
           F^(a)_ UV-finite    =     F^(a) - ∑_ UV regions i _i  F^(a) + ⋯  ,
      
      where _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. (<ref>) and (<ref>), 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, 𝐏_1.
      
      In our treatment, the ultraviolet counterterms in Eq. (<ref>)
      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 
      𝐏_1 T 𝐏_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. (<ref>).   
      
      Finally, before going on to the treatment of the one-loop amplitudes, Eq. (<ref>), we will address the motivation for the use of the projectors 𝐏_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. <ref>. 
      Because the
      integral in F^(1) is fully decoupled from M^(0), the fermion propagators that are contracted with the local vertex in Fig. <ref> 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,
      
          q_e^ν ∑_ q_e insertions M_ν^(0,Π)   =     p_2  ( p_2-  q_e  )^-1∏_c=2, e^n ϵ^*(q_c)  ( p_1 - ∑_d=1, e^c-1 q_d  )^-1
                    + ∏_c=2, e^n ϵ^*(q_c)  ( p_1 - ∑_d=1, e^c-1 q_d  -  q_e )^-1  ( p_1 -  q_e  )^-1 p_1  .
       
      If this expression  is inserted in Eq. (<ref>)  without the projectors 𝐏_1, both terms will vanish only if both virtual fermion lines in Fig. <ref> 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 gauge-dependence vanishes identically because p_1^2= 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. (<ref>) and (<ref>), in which all infrared singularities are
      subtracted as form factor integrands (the F^(L) in Eq. (<ref>)), 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. (<ref>).
      M^(1) is singular in d=4 space-time 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 fully-integrable and process-dependent
      result, M^(1)_ finite in Eq. (<ref>), 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. (<ref>).
      
      Our aim is to identify a set of counterterms  that are local in momentum space, and which add up to
      a single function  M_ singular^(1)(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,
      
          M^(1)   = ∫d^d l/(2 π)^d[
               M^(1)(l)
          -  M_ singular^(1)(l)
          ]
          +∫d^d l/(2 π)^d M_ singular^(1)(l)  .
      
      The first integral on the right-hand side, which we will identify with M_ finite^(1) in Eq. (<ref>),
      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. <cit.>. 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. (<ref>) 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. <cit.>.  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. <ref>, 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. (<ref>).
      We will study their UV structure in the following subsection.
      
      We may decompose the one-loop amplitude into four classes,
      
          M^(1)≡[.28]c19-4 = M^(1,A) + M^(1,B_1) + M^(1,B_2) +M^(1,C)  ,
      
      each of which is characterized by the following diagrams:
      
        *   the virtual photon is adjacent to both incoming fermions,
      
          M^(1,A)≡[.25]M1Adiags+ external photon permutations ,
      
        *  the virtual photon is adjacent to the incoming electron but not
        to the incoming positron, 
      
          M^(1,B_1)≡[.25]M1B1_diag1 +[.25]M1B1_diag2+ … ,
      
      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, 
      
          M^(1,B_2)≡[.25]M1B2_diag1 +[.25]M1B2_diag2+ … ,
      
        *  and, finally, the virtual photon is not adjacent to either of the incoming fermions
      
          M^(1,C)≡[.25]M1C_diag1 +[.25]M1C_diag2+ … .
      
      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
      
        *  a soft singularity l = 0 in the graphs of M^(1,A), 
      
        *  a collinear singularity l = -x_1 p_1, 0≤ x_1 ≤ 1, in the graphs of
        M^(1,A) and M^(1,B_1),
      
        *  a collinear singularity l =  x_2 p_2, 0≤ x_2 ≤ 1, in the graphs of
        M^(1,A) and  M^(1,B_2).
      
      We refer to these points or lines as “pinch surfaces" <cit.>.
      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. <cit.>.    A recent
      and very general discussion has been given  in Ref. <cit.>.
       
      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. (<ref>).    As we will show, as indicated in Eq. (<ref>), 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. (<ref>), the annihilation form factor e^++e^- →γ^*.
      
      We now seek a function  M_IR^(1) 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 <cit.>. 
      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. <cit.>.  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 
      
          M^(1,A) = -i e^2/l^2 (p_1+l)^2 (p_2-l)^2 N^(1)  ,
      
      where the numerator N^(1) reads
      
          N^(1) =  v̅(p_2) γ^μ( l  -  p_2 ) 
          ℳ^(0) (p_1 + l, p_2 - l, q_1, …, q_n)
          ( p_1+  l ) 
          γ_μ u(p_1).
      
      In the soft limit l^μ→ 0, we neglect l inside ℳ^(0), and approximate M^(1) with
      
          M^(1)   _ soft M^(1)
              = i e^2
            v̅(p_2) γ^μ( l -  p_2  ) 
          ℳ^(0) (p_1, p_2, q_1, …, q_n)
          ( p_1+  l ) 
          γ_μ u(p_1)
            /l^2 (p_1+l)^2 (p_2-l)^2  .
      
      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 ℳ^(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 <cit.>.
      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
      
          M^(1,A+B_1)   ≡ M^(1,A) +  M^(1, B_1)
              =   -i e/l^2 (p_1+l)^2 v̅(p_2)
          η^μν ℳ^(0,A+B_1)_μ (p_1+l, p_2, l;q_1, …, q_n) (  p_1 +  l ) γ_ν u(p_1) .
      
      In this relation, ℳ^(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. (<ref>).
      The polarization vector of the exchanged photon (momentum l) as well  as the electron and positron spinors are removed from ℳ^(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
      
          q^μ≈q ·η_1 /l ·η_1  l^μ ,
      
      or alternatively,
      
          q^μ≈2 q ·η_1 /(l + η_1)^2 -η_1^2 l^μ.
      
      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. (<ref>), with a denominator linear in
      l, suffices in order to demonstrate the factorization of collinear singularities (see, for example, Ref. <cit.>), in
      this article  we use the equivalent version of
      Eq. (<ref>) 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 
      
          η^μν =l^μη_1^ν/l ·η_1 +l^νη_1^μ/l ·η_1   -η_1^2/(l ·η_1)^2
          l^μ l^ν+ η_⊥^μν  ,
      
      where η_⊥^μν  l_ν = η_⊥^μν η_1_ν =0. 
      Staying always within the collinear approximation, where l^2 → 0, we can equivalently approximate the
      metric with
      
          η^μν≈2 l^μη_1^ν/(l
              +η_1)^2-η_1^2 +2 l^νη_1^μ/(l +η_1)^2-η_1^2   -4η_1^2/(  (l + η_1 )^2 - η_1^2  )^2
          l^μ l^ν+ η_⊥^μν  .
      
      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. (<ref>) can be re-written as, 
      
          (  p_1 +  l ) γ_ν  u(p_1) ≈ 2 (p_1+l)_ν  u(p_1) ∝ p_1 ν  .
      
      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. (<ref>) according to Eq. (<ref>),
      we notice that only one term gives a non-vanishing
      contribution in the collinear limit.
      Thus, in the limit l ≈ -x p_1  we can make the replacement η^μν→ 2 η_1^ν  l^μ / [(l+η_1)^2 -η_1^2] without changing the behavior of M^(1,A+B_1) in the region l || p_1. This allows us to write
      
          M^(1,A+B_1)   ≡ M^(1,A) +  M^(1, B_1)≈
             -i e/l^2 (p_1+l)^2 v̅(p_2) [2η_1^ν  l^μ/(l + η_1)^2  - η_1^2 ℳ^(0,A+B_1)_μ (p_1+l, p_2, l;q_1, …, q_n) ] 
          
             × (  p_1 +  l ) γ_ν u(p_1).
      
      We represent the right-hand side of Eq. (<ref>) graphically by 
      
          e^2 M^(1,A+B_1)   ≈[.25]M1Adiag-coll1_RH + [.25]M1B1diag-coll3 
             + [.25]M1Adiag-coll2_RH + … .
      
      The triangle arrow at the lower end of photon lines denotes the following polarization approximation in the photon propagator: 
      
          [.15]triangleNotation  :     -i/q^2η^μν→-i/q^22 η_1^ν  q^μ/(q+η_1)^2 -η_1^2 ,
      
      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,
      
          [.15]triangleNotation2  :     -i/q^2η^μν→-i/q^2-2 η_2^ν  q^μ/(q-η_2)^2 -η_2^2.
      
      The expression in Eq. (<ref>) is now ready to be simplified with a partial fractioning
      identity, which is a manifestation of the QED Ward identity,
      
          2η_1^ν/(q + η_1)^2 - η_1^2[ i/p-q (-ie  q) i/p]
            = e 2η_1^ν/(q + η_1)^2 - η_1^2[i/ p -  q - i/ p]  .
      
      This identity describes the effect of attaching a photon with a
      scalar polarization on a fermion line.
      We visualize the identity in Fig. <ref>.
      
      We also note the following identities for insertions of scalar photons in external fermion lines:
      
          2η_1^ν/(q + η_1)^2 - η_1^2[
            i/p + q (-ie  q) u(p) ]
               =e 2η_1^ν/(q + η_1)^2 - η_1^2
              u(p)   , 
          2η_1^ν/(q+ η_1)^2 - η_1^2[ v̅(p) (-ie  q)  i/p - q]    = -e 2η_1^ν/(q + η_1)^2 - η_1^2v̅(p),
      
      which are depicted pictorially in Fig. <ref>.
      
      We now see that repeated applications of Eqs. (<ref>) - (<ref>) to
      Eq. (<ref>) result in pairwise cancellations of all terms except one,
      
          e^2 M^(1,A+B_1)≈-i e^2/l^2 (p_1+l)^2v̅(p_2) [2η_1^ν/(l + η_1)^2 - η_1^2 ℳ^(0)]
            (  p_1 +  l ) γ_ν u(p_1)  ,
      
      where ℳ^(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. <ref>.
      
      With the action of the Ward identity, non-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. (<ref>) holds generally.
      
      Interestingly, if we choose η_1 to be equal to -p_2 in Eq. (<ref>), then the three 
      denominators become identical to those of the soft-photon
      approximation given by  the right-hand side of
      Eq. (<ref>).  After
      some simple Dirac algebra, it can be shown that the numerator of Eq. (<ref>) also coincides with that of
      Eq. (<ref>) 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. (<ref>), 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. <ref>, we will place the tree amplitude between the Dirac
      projectors defined in Eq. (<ref>).  We arrive at the following global IR counterterm,
      
          _ IR M^(1)   = i e^2 v̅(p_2) γ^μ(    p_2  -  l )
              𝐏_1   ℳ^(0)(p_1, p_2; q_1, …, q_n)  𝐏_1 
              (  p_1  +  l ) γ_μ   u(p_1)
            /l^2 (p_1+l)^2 (l-p_2)^2  ,
      
      which is equivalent to the soft photon approximation in Eq. (<ref>) because 𝐏_1 acts as the identity in the soft, and collinear, regions.
      This is the usual one-loop form factor, which we will denote by ℱ^(1), but with the QED vertex replaced by the truncated tree-level spinor matrix sandwiched between a pair of projectors, just as in Eq. (<ref>). Therefore, the one-loop IR counterterm is equivalent to
      
          _ IR M^(1)   ≡ℱ^(1)[𝐏_1  ℳ^(0)(p_1, p_2; q_1, …, q_n)  𝐏_1 ]  
          ≡[0.3]c19-5 ,
      
      where we have recalled the corresponding diagrammatic notation of Fig. <ref>. 
      
      
      The IR subtraction with the form factor counterterm of Eq. (<ref>) 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 𝐏_1 projector, as defined by
      Eq. (<ref>), in Eq. (<ref>) 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) or v̅(p_2), the projector 𝐏_1 has the following two alternative representations,
      
          𝐏_1    = 1/2p_1 · p_2 p_1 ( v_L (p_2) v̅_L(p_2) + v_R (p_2) v̅_R(p_2) ), 
          𝐏_1    = 1/2p_1 · p_2( u_L (p_1) u̅_L(p_1) + u_R (p_1) u̅_R(p_1) )  p_2  .
      
      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 replace v̅(p_2) by v̅_R(p_2) in Eq. (<ref>). (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,
      
          1 - γ^5/2 u_L(p_1) = u_L(p_1),   1 + γ^5/2 u_L(p_1) = 0, 
             v̅_R (p_2) 1 + γ^5/2 = v̅_R (p_2),   v̅_R (p_2) 1 - γ^5/2 = 0   .
      
      The IR counterterm _ IR M^(1) in Eq. (<ref>) is rewritten as
      
          _ IR M^(1)    =
             -ie^2/l^2 (p_1+l)^2 (l-p_2)^2 (2 p_1 · p_2)^2 v̅_R(p_2) γ^μ ( p_2 -  l)
             p_1 
               × ( v_L (p_2) v̅_L(p_2) + v_R (p_2) v̅_R(p_2) )
            ℳ^(0)(p_1, p_2, q_1, …, q_n) 
               × ( u_L (p_1) u̅_L(p_1) + u_R (p_1) u̅_R(p_1) )  p_2
            (  p_1 +  l ) γ_μ u_L(p_1)  .
      
      Since an odd number of gamma matrices preserves chirality, the above expression 
      reduces to 
      
          _ IR M^(1)    =
             -ie^2/l^2 (p_1+l)^2 (l-p_2)^2 (2 p_1 · p_2)^2 v̅_R(p_2) γ^μ ( p_2 -  l)
             p_1  v_R (p_2)
               × [ v̅_R(p_2) ℳ^(0)(p_1, p_2, q_1, …, q_n) u_L (p_1)]  u̅_L(p_1) 
             p_2 (  p_1 +  l ) γ_μ u_L(p_1)   .
      
      The square bracket in this equation is precisely the tree amplitude
      contracted with the external spinors v̅_R(p_2) and u_L(p_1),
      so the counterterm of Eq. (<ref>) 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,
      
          ℳ^(1)_IR-finite
          = ℳ^(1) - _IRℳ^(1)  ,
      
      with _IRℳ^(1) given by the explicit subtraction of Eq. (<ref>), or equivalently, in form factor notation by Eq. (<ref>).
      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. (<ref>), 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 _UVℳ_IR-finite^(1) that matches the singular behavior of the integrand in the UV-limit,
      
          _UVℳ_IR-finite^(1)   = _UV( ℳ^(1) - _IRℳ^(1)) 
             = _UVℳ^(1) - _UV_IRℳ^(1)  ,
      
      with _IRℳ^(1) given in Eq. (<ref>).   It is at this point that we will encounter UV divergences associated with
      the factorized amplitude, Eq. (<ref>).
      
      At one loop in the process under study, UV divergences occur only in triangle diagrams and fermion self-energies, 
      which we will denote by Γ_eeγ^(1),ν and Π_e^(1), respectively.  
      Following our convention for the routing of the loop momenta, these
      Green's functions have the integrands
      
          e Γ_eeγ^(1),ν(p, q, l) ≡[0.25]vertex_1L 
          = (- e^3 ) γ^μ(p+l+q)γ^ν(p+l)γ_μ/(p+q+l)^2(p+l)^2l^2  ,
      
      and
      
          Π_e^(1)(p, l) ≡[0.21]selfenergy 
          = -e^2 γ^μ(l+p)γ_μ/l^2(l+p)^2  ,
      
      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 
      
          q_ν Γ_eeγ^(1),ν(p, q, l) 
          = Π_e^(1)(p, l) - Π_e^(1)(p+q, l)   .
      
      For the vertex diagram, we choose, as suggested above, a UV counterterm that is defined by its integrand, given in this case by
      
          e Γ_eeγ^(1,UV),ν(l) ≡[0.25]vertexUV_1L
          = (- e^3 ) γ^μlγ^νlγ_μ/(l^2-M^2)^3  ,
      
      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
      
          Π_e^(1,UV)(p, l)    ≡   [0.12]selfenergyUV_1L
          = (-e^2 ) [ 
          γ^μlγ_μ/(l^2-M^2)^2 - γ^μlplγ_μ/(l^2 - M^2)^3] 
          
                 = -e^2 γ^μlγ_μ/(l^2-M^2)^2 
          -  p_νΓ_eeγ^(1,UV),ν(l)   .
      
      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, 
      
          q_νΓ_eeγ^(1,UV),ν(l) 
          =  Π_e^(1,UV)(p, l) -  Π_e^(1,UV)(p+q, l),
      
      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 sub-divergences. 
      
      Note that the Ward identity for the UV counterterms Eq. (<ref>) is by no means guaranteed by the Ward identity for the original diagrams Eq. (<ref>), 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 <cit.>, is
      
          -e^2 γ^μ ( l +  p) γ_μ/[ (l + p/ 2 )^2 ]^2  .
      
      This alternative counterterm perfectly matches the UV-divergent behavior of Eq. (<ref>) 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. (<ref>).
      We have chosen non-divergent terms so that the UV counterterm for the self energy, Eq. (<ref>), is aligned with our UV counterterm for the vertex, Eq. (<ref>), for both finite and divergent parts, preserving the Ward identity locally.
      [A clean way to arrive at Eqs. (<ref>) and (<ref>) is to perform a series expansion of Eqs. (<ref>) and (<ref>) in the limit of large l, and truncate both series at the order that corresponds to a logarithmic divergence after loop integration, before adding a mass regulator to every propagator. The uniform truncation of the series preserves the Ward identity that relates the vertex and the self energy, and the final step of adding an IR mass regulator again preserves the Ward identities.]
      
      We have now specified the local UV-counterterms of the one-loop amplitude,
      
          e^2 _UVℳ^(1)   =
          [0.25]UVd1 
          +[0.25]UVd2
          +[0.25]UVd3
             +… 
          +[0.25]UVd4 
          +  external photon permutations  .
      
      Finally, we need to find a UV approximation for the infrared
      counterterm, _IRℳ^(1) in
      Eq. (<ref>).  We choose
      
          e^2 _UV_IRℳ^(1)   ≡v̅(p_2) [
          [0.2]UVFF1L] u(p_1)
          
             = v̅(p_2)
          [ 
          (-i e^2 ) γ^μl𝐏_1
            M^(0)(p_1, p_2; q_1, …, q_n ) 
          𝐏_1
            lγ_μ/(l^2-M^2)^3]
          u(p_1)
            .
      
      In the above, we have drawn a square vertex with an extra cross to denote the UV limit of the triangle diagram in Eq. (<ref>) with a form factor vertex 𝐏_1 ℳ^(0)𝐏_1, analogous to how the diagram in Eq. (<ref>) denotes the UV limit of the vertex diagram in Eq. (<ref>), with a label “(1)” indicating that it is a 1-loop counterterm. Again the projector 𝐏_1 is always omitted in the diagrammatic notation.   In terms
      of the relation, Eq. (<ref>) that defines M^(1)_ finite, this counterterm will produce a contribution to F^(1)_ UV-finite.
        
      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, 
      
          ℳ^(1)(l) = ℳ^(1)_ finite(l) +  ℳ^(1)_ singular(l)  .
      
      The first term on the right-hand side of Eq. (<ref>), ℳ^(1)_ finite(l), is integrable in
      four space-time dimensions.  All singularities are now contained in
      
          ℳ^(1)_ singular(l) = _IRℳ^(1)(l) 
          +_UVℳ^(1)(l) 
          -_UV_IRℳ^(1)(l).
      
      Expressions for  the counterterms on the right-hand side of  
      Eq. (<ref>) have been given by
      Eqs. (<ref>) for the IR counterterm, Eqs. (<ref>), (<ref>) and (<ref>) for the UV subtraction 
      and (<ref>) for the combined UV/IR term.     
      The integration of the counterterms in ℳ^(1)_ singular(l) over 
      loop momentum l can be performed simply in d=4-2ϵ
      dimensions. These results can then be combined with the 
      numerical evaluation of ℳ^(1)_ finite(l).
      In terms of the general relation, Eq. (<ref>), we can identify
      
          M^(1)_ UV-finite    =  M^(1) - _UV M^(1)  ,
          
           F^(1)_ UV-finite    = _ IR M^(1) - _UV_IRℳ^(1)
             =  F^(1) - _UVℱ^(1)  ,
      
      where in the second relation for F^(1)_ finite we have used Eq. (<ref>).  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. (<ref>).
      
      Let us consider first the UV counterterms for the one-loop vertex and
      propagator, which appear in Eq. (<ref>). 
      The vertex counterterm integrates to,
      
          ∫d^d l /(2 π)^d  e Γ_eeγ^(1,UV),ν(l)  
          = Z̃_1^(1)( -ie γ^ν),
      
      or, diagrammatically, 
      
          ∫d^d l /(2 π)^d[0.2]vertexUV_1L   =
          Z̃_1^(1)[0.2]vertex1gam  ,
      
      with 
      
          Z̃_1^(1) = e^2/(4 π)^d/2Γ(1+ϵ)  (1-ϵ)^2/ϵ(M^2)^-ϵ  .
      
      For the integral of the electron-propagator counterterm we find 
      
          ∫d^d l /(2 π)^di/ pΠ_e^(1,UV)(p, l) i/ p 
          =  - Z̃_1^(1)i/ p,
      
      or, diagrammatically, 
      
          ∫d^d l /(2 π)^d[0.2]dressedfermionprop   =
          - Z̃_1^(1)[0.2]simplefermionprop.
      
      As expected, every
      vertex  or  propagator counterterm in
      Eq. (<ref>) 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 (<ref>).
      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,
      
          e^2 _UV  M^(1)≡ e^2 ∫d^d l/(2 π)^d_UVℳ^(1)
          = Z̃_1^(1) M^(0).
      
      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. (<ref>),
      
          e^2 _UV_IR  M^(1)   ≡ e^2 ∫d^dl/(2 π)^d _UV_IRℳ^(1)(l)  
             = Z̃_1^(1)/4 (1-ϵ)^2 v̅(p_2) γ^μ γ^ν 𝐏_1M^(0)𝐏_1 γ_ν γ_μ u(p_1) .
      
      Meanwhile, the integration of the IR counterterm yields
      
          _IR  M^(1)   ≡∫d^dl/(2 π)^dℱ^(1)(l)  
             =
          -1/(4π)^d/2Γ() Γ(1-)^2/Γ(2-2) (-s)^-(  [ 2/  + /(1-)]   v̅(p_2) 𝐏_1  M^(0)𝐏_1 u(p_1) . 
             .  - 1/4(1-) v̅(p_2) γ^μγ^ν 𝐏_1  M^(0)𝐏_1  γ_νγ_μ u(p_1)  )  ,
      
      with s=(p_1+p_2)^2.
      Combining this result with the ultraviolet counterterm of Eq. (<ref>) and the UV-IR combination counterterm
      of Eq. (<ref>) to get the full singular contribution to the amplitude in Eq. (<ref>), we find that 
      the divergent part of the bare one-loop amplitude is
      
          M^(1)_ singular   = 
           -1/(4π)^d/21/Γ(1-) (-s)^-
               ×{(2/^2 + 3/ +ln(M^2/-s) + 11 ) v̅(p_2) 𝐏_1  M^(0)𝐏_1 u(p_1) 
          
             -1/4(ln(M^2/-s) + 3) v̅(p_2) γ^μγ^ν 𝐏_1  M^(0)𝐏_1  γ_νγ_μ u(p_1)
          }
          +𝒪() 
           .
      
      We note that 
      
          v̅(p_2) 𝐏_1
            M^(0)𝐏_1 u(p_1) = v̅(p_2) M^(0) u(p_1) = M^(0),
      
      and the divergent part as ϵ→ 0 is proportional to the
      tree-amplitude in agreement with previous expectations <cit.>.  
      
      
      
      
      § 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. (<ref>), 
      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. (<ref>), we start
      with M^(2)_2, consisting of two-loop diagrams with
      a photon vacuum polarization subgraph, 
      as in Fig. <ref>. 
      
      The integrand M^(2)_2 can be generated from the diagrams of
      the one-loop amplitude in Eqs. (<ref>) - (<ref>),
      by replacing the free photon propagator with its one-loop correction,
      
          -i g_μν/ l^2 → [.3]photonprop1L≡-i / l^2( -i Π_μν(l) )  -i /  l^2  ,
      
       with
       
          -i Π_μν(l) =  -  e^2  tr[ ( k +  l
              )  γ_μ  k  γ_ν] /k^2 (k+l)^2  .
      
      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  
      
          -i / l^2( -i Π_μν(l) )  -i /
              l^2∼δ^-4,
      
      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. (<ref>) 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 <cit.>. 
      After this reduction, the integrand for the
      one-loop photon propagator reads,
      
          photonprop1L
            = 
            1/l^22 (d-2)/d-1  e^2 ( η_μν    -  l_μ l_ν/l^2)
          1 /k^2 (k+l)^2  .
      
      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. (<ref>), 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.[Scaleless integrals vanish in dimensional regularization due to cancellation of UV and IR poles. Dropping scaleless integrals mixes UV and IR poles, but this is not a problem because the universal structure of UV and IR divergences of gauge theory amplitudes are known.] After doing so, the tensorial form of the
      effective one-loop virtual propagator matches the one of the tree photon-propagator,
      
          photonprop1L→η_μν/l^22 (d-2)/d-1 1/k^2 (k+l)^2  .
      
      Using this result, it is simple to  generate a
      suitable integrand for M^(2)_2(k,l) from the corresponding
      one-loop amplitude integrand M^(1)(l) as 
      
          M^(2)_2(k,l) 
             =    
          i 2 (d-2)/d-11 /k^2 (k+l)^2  M^(1)(l)    .
      
      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. <ref>.
      
      Infrared singularities in Eq. (<ref>) are
      approximated by the counterterm
      
          _ IR M^(2)_2(k,l)   =   i 2 (d-2)/d-1 1 /k^2 (k+l)^2 _ IR M^(1)(l)
          
             ≡     F^(2)_2(k,l)  ,
      
      where the “form factor” diagram _ IR
        M^(1)(l) =  F^(1)(l)  is given by Eq. (<ref>). 
        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. (<ref>).
        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 ℳ^(2)_2 - _ IRℳ^(2)_2
      as k →∞ or l →∞, as suggested in Eq. (<ref>).
       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, 
      
          M^(2)_2,  finite   =    
          i 2 (d-2)/d-1 [ 
          1 /k^2 (k+l)^2
          - 
          1 /( k^2-M^2)^2] 
           M^(1)_ finite(l)  .
      
      Here, the finite one-loop result is specified by Eqs. (<ref>) and (<ref>),
      
          M^(1)_ finite(l) =  M^(1)(l) 
          - _ UV M^(1)(l) 
          - _ IR M^(1)(l)
          + _ UV_ IR M^(1)(l)  .
      
      Notice that M^(2)_2,  finite in Eq. (<ref>) 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. (<ref>), 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. (<ref>), using (<ref>),
      shows that it is precisely of the form anticipated in Eq. (<ref>) in terms of our modified amplitude and form factor integrands.
      
      In summary, we have decomposed the M^(2)_2
      coefficient to the two-loop amplitude integrand into a remainder  that can 
      be integrated directly in four dimensions,  and a simple singular part 
      that can be integrated  in d=4-2ϵ dimensions,     
      
          M^(2)_2 = M^(2)_2,   finite+ M^(2)_2,  singular  ,
      
      with M^(2)_2,   finite given by Eqs. (<ref>) and (<ref>).
      
      
      
      
       §.§ 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. <ref>. 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.[It is not important which one of the four lines in the box is chosen as k, and the choice does not have to be aligned between different diagrams.]  
      Due to Furry's theorem, amplitude coefficients for odd values of κ vanish.  We do not consider diagrams in which a single virtual
      photon decays into an odd number of final-state photons (κ=c-1), because the IR divergences of this set of diagrams are identical to those of
      a one-loop diagram with the emission of n-κ+1 photons, at least one of which is off-shell.
      
        
      Some diagrams that contribute to M^(2)_κ+2, κ>0, illustrated by 
      Fig. <ref>, 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  
      
          M^(2)_κ+2   =   1/l^2
            (p_1+l)^2v̅ (p_2) A_ν(q_κ+1, …, q_n)
          ( p_1 + l )
          γ_μ u(p_1)
          
                × M_ photons^μν(q_1, …, 
            q_κ,  l,  -∑_j=1^κq_j-l ) .
      
      In the collinear limit l ≈ - x p_1, the last few factors times the Dirac spinor u(p_1) become
      
          ( p_1 + l ) γ_μ u(p_1) ≈ 2(1-x) p_1 μ  u(p_1).
      
      Since this  is proportional to p_1μ, we can apply the approximation of Eqs. (<ref>) and (<ref>) with q replaced by l to rewrite Eq. (<ref>) as
      
          M^(2)_κ+2|_l ∥ p_1   →   1/l^2
            (p_1+l)^2v̅ (p_2) A_ν(q_κ+1, …, q_n)
          ( p_1 + l )
           2 η_1/(l+η_1)^2 - η_1^2 u(p_1)
          
                × l_μ  M_ photons^μν(q_1, …, 
            q_κ,  l,  -∑_j=1^κq_j-l )  .
      
      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. <ref> and Eq. (<ref>), 
      
          2 η_1/(l+η_1)^2 - η_1^2 l_μ M_ photons^μν(k, q_1, q_2, l, l̃) 
              =  [0.25]LbL2 + [0.25]LbL1 + [0.25]LbL3
              + 
             
          = e  ( [0.25]LbL2p_RH - [0.25]LbL2m_RH) 
             
          + e ( [0.25]LbL1p_RH - [0.25]LbL1m_RH)
          
             
          + e ( [0.25]LbL3p_RH - [0.25]LbL3m_RH)
          + 
             
           = e ( [0.25]LbL2p_RH - [0.25]LbL3m_RH)
          +   .
      
      We have used the shorthand notation l̃=-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, 
      
          2 η_1/(l+η_1)^2 - η_1^2 l_μ M_ photons^μν(k, q_1, q_2, l, l̃)     =   2 η_1/(l+η_1)^2 - η_1^2[  
           M_ photons^ν( k, q_1, q_2, -q_1-q_2)
           . 
                . 
          -  M_ photons^ν( k+l, q_1, q_2, -q_1-q_2)
          ]  ,
      
      where M_ photons^ν is defined in Eq. (<ref>).
      This difference integrates to zero in dimensional
      regularization, which respects momentum-shift invariance. Therefore
      
          ∫d^d k/(2π)^d2 η_1/(l+η_1)^2 - η_1^2  l_μ M_ photons^μν(k, q_1, q_2,… , l,
            l̃) = 0.
      
      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. (<ref>) and modify the integrand of Eq. (<ref>) as, 
      
          M^(2)_κ+2   =   1/l^2
            (p_1+l)^2v̅ (p_2) A_ν(q_κ+1, …, q_n)
          ( p_1 + l )
          γ̂_μ(l, η) u(p_1)  
          
                × M_ photons^μν(q_1, …,
            q_κ, l, -∑_j=1^κq_j-l )   .
      
      Here, we have modified the vertex adjacent to the incoming electron line, replacing γ_μ by
      
          γ̂_μ(l, η_1 ) 
          ≡γ_μ - l_μ2 η_1/(l+η_1)^2 - η_1^2  ,     η_1^2 ≠ 0  ,
      
      where η_1 is chosen so that it produces no new pinches from the extra denominator.
      In the collinear limit, l= -x p_1, we have 
      
          ( p_1 + l )
          γ̂_μ(l, η_1 ) u(p_1)  
            0.
      
      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, 
      
          [0.6]C1C2floop_mod
             = v(p_2) γ̂_ν(l_b, -η_2) 
          ( l_b - p_2 ) 
          …( p_1 + l_a ) 
          γ̂_μ(l_a, η_1)   M_ photons^μν .
      
      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. (<ref>).
      We note that in this case, there is no corresponding modification to a
      form factor subtraction term in Eq. (<ref>).   The modified 
      contributions are now infrared-integrable on their own, without further subtraction.
      
      In summary, we have modified the integrand of M_κ+2^(2), to make each of its diagrams
      regular in all its collinear limits, without changing its integral.   This was done with a simple
      modification of the fermion-photon vertex when it occurs adjacent to an
      incoming leg. This class of diagrams has no further infrared divergences. 
      
      
      
       §.§ UV counterterms for c=4
      
      
      In QED,  the fermion loops in M^(2)_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. <cit.>.  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
      
          V^μ_1 μ_2 μ_3 μ_4 (q_1, q_2, q_3, -q_1-q_2-q_3)
             = e^4 (-1) tr[ k γ^μ_1 ( k +  q_1) γ^μ_2 ( k +  q_1 +  q_2) γ^μ_3 ( k +  q_1 +  q_2 +  q_3) γ^μ_4]/k^2 (k+q_1)^2 (k+q_1+q_2)^2 (k+q_1+q_2+q_3)^2  ,
      
      as illustrated in Fig. <ref>.
      
      In the UV limit, we write down the counterterm just as for the vacuum polarization in Eq.  (<ref>), by dropping all external momenta 
      and IR-regulating denominators with a mass M,
      
          _ UV V^μ_1 μ_2 μ_3 μ_4 = -e^4  tr( 
          kγ^μ_1kγ^μ_2kγ^μ_3kγ^μ_4)
          /( k^2-M^2 )^4
            .
      
      The counterterm in this expression corresponds to the specific subdiagram of 
      Eq. (<ref>).   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^(2) is just the sum of these six distinguishable
      diagrams, each minus its counterterm.   Each one
      of these combinations specifies a locally integrable two-loop integral, which is 
      also UV convergent.   
      
      The corresponding contribution to  M_4,  singular^(2) is just the sum of the six counterterms,
      Eq. (<ref>) plus permutations.   Their sum can be combined into a single, local vertex 
      embedded in the one-loop diagram in Fig. <ref>.
      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,
      
          ∑_ permutations_ UV V^μ_1 μ_2 μ_3 μ_4 = 
          8 e^4 ϵ (1-ϵ)/(2-ϵ) (3-ϵ)(k^2)^2/( k^2-M^2 )^4
          
                ×[ 
          η^μ_1 μ_2η^μ_3 μ_4
          + 
          η^μ_1 μ_3η^μ_2 μ_4
          + 
          η^μ_1 μ_4η^μ_2 μ_3]  .
      
       The integral of loop momentum k has a pole, so that the sum of counterterms is finite in d=4, 
       
      
          ∑_ permutations∫d^d k/(2 π)^d_ UV V^μ_1 μ_2 μ_3 μ_4
          =i e^4/12 π^2[ 
          η^μ_1 μ_2η^μ_3 μ_4
          + 
          η^μ_1 μ_3η^μ_2 μ_4
          + 
          η^μ_1 μ_4η^μ_2 μ_3]  .
      
      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. <ref> 
      is locally IR finite and UV convergent.   Figure <ref> simply produces finite
      shifts in M^(2)_4,  finite and M^(2)_4,  singular, with the same magnitudes but with opposite
      signs.   Its only effect is to enable us to express  M^(2)_4,  finite as
      a sum of locally integrable and UV convergent diagrams.
      
      
      
      
      To summarize the results of this section, we have identified integrable 
      expressions for the sets of two-loop diagrams with fermion loops, { M^(2)_c}.  For those with
      a vacuum polarization, M^(2)_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. (<ref>), which is equivalent
      to the subtraction specified in Eq. (<ref>), 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. (<ref>), 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
      #1#2=#2#1p/k/l/q/ϵ/Q/
      
      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. (<ref>), 
      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 <ref> 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,
      
          M^(2)_0(k, l)    =   [.25]m2a1_RH 
          +[.25]m2a2_RH 
          +[.25]m2a3_RH
                
          +[.25]m2a4_RH 
          +[.25]m2a5_RH 
          +[.25]m2a6_RH
                 
          +[.25]m2a7_RH 
          +[.25]m2a8_RH + …  .
      
      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
      
          k↔ l,    (k_⊥, l_⊥) ↔ (-k_⊥, -l_⊥),
      
      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, 
      
          M^(2)_ sym(k, l)    ≡   1/4[ 
            M_0^(2)(k, l)
          + M_0^(2)(l, k)
          +  M_0^(2)(k^*, l^*)
          + M_0^(2)(l^*, k^*) ]  ,
      
      with l^*=l(l_⊥→ -l_⊥), k^*=k(k_⊥→ -k_⊥), 
      exhibits more complete factorization properties than the individual terms that make it up
      on the right-hand side.   Here, the notation M_0^(2) refers to the integrands
      written down directly from the QED Lagrangian, which are symmetrized as shown to define M_ sym. 
      
      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,
      
          M_0^(2)(k, l)    [.17]M2p2a
          +[.17]M2p2b
          + [.17]M2p2c
          + [.17]M2p2d
                  + …
             = [.17]M2p2e + …
      
      In the above,  grey triangles denote the collinear approximation for photon propagators collinear to p_2, defined in Eq. (<ref>).
      The identity follows, as at one loop, by successive application of the lowest order Ward identity, Eq. (<ref>).
      In Eq. (<ref>), 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,
      
          M_0^(2)(k, l)    [.17]M2p1a
          +[.17]M2p1b
          + [.17]M2p1c
          + [.17]M2p1d
                  + …
             = [.17]M2p1e + …  .
      
      In Eq. (<ref>), 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. (<ref>). 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. (<ref>) is actually also needed in Eq. (<ref>), 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_0^(2)(k, l) +  M_0^(2)(l, k) without a change of variables in momenta.  
      Once the momenta are averaged in this manner, each momentum assignment will find
      its matching assignment in every diagram. 
      
      The additional reflections (k_⊥,l_⊥) ↔ (-k_⊥, -l_⊥) in
      the transverse plane in Eq. (<ref>) are required for the factorization of
      collinear singularities in propagators which receive hard self energy
      and vertex corrections. 
      
      The diagrams in question are shown in Fig. <ref>, 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 , 
      and the diagrams with a QED vertex correction outside the gray
      photon emission blob as .
      For concreteness, we will assign l to be the momentum of the virtual photon of the
      self energy or vertex sub-loop in  and  diagrams as in
      figure <ref>. This choice is simply a convention, since the integrand will be
      symmetrized according to Eq. (<ref>) 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,
      
          ℳ_ sym^(2) =  [ ℳ^(2)_ + ℳ^(2)_ + ℳ^(2)_regular ]_ sym  ,
      
      where M_ sym^(2) denotes the full initial 
       integrand with fixed momentum assignments and symmetrizations in Eq. (<ref>) for each of the terms on the right.
      We will not alter the integrand of ℳ^(2)_regular
      beyond these momentum symmetrizations, while ℳ^(2)_ and ℳ^(2)_ will
      require a somewhat more elaborate treatment. 
      
      
      
        §.§.§ Collinear regions in   and  diagrams
      
      
      The diagrams in ℳ^(2)_ and  ℳ^(2)_ of Fig. <ref>, 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.
      
      
        *  Stronger than logarithmic divergences due to doubled fermion propagators,  
      
        *  apparent contributions of unphysical “loop polarizations”  that are not compatible with factorization in collinear limits,
      
        *  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 , 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 [A related discussion of
        double propagators can be found in Ref. <cit.>.]. 
      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. <ref>, 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.
      
      
      The second challenge posed by “loop polarizations” emerges, for
      example, in the  calculation of the collinear limit k → - x p_1
      in  diagrams. Let us focus  on the vertex subgraph given by
      Eq. (<ref>)  multiplied with a common fermion
      propagator and a spinor factor,  as they appear in  diagrams,
      
          [.35]type2v_RH
            = -i   e^3  (p_1 + k)γ^μ (p_1+l+k)γ^ν(p_1+l)γ_μ u(p_1)/(p_1+l+k)^2(p_1+l)^2l^2
               
                           -i  e^3   2  (2-d) 
          1+ 2 k · p_2/s/
          (p_1+l+k)^2
          [
          l^ν/l^2-
          l^ν+p_1^ν/(l+p_1)^2] u(p_1)  .
      
      The tensor structure of the subgraph in the collinear limit  is
      crucial for factorization.
      Analogous to the one-loop case in Eq. (<ref>), 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  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 (<ref>) 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  or ) factorize
      straightforwardly by Ward identities after using the collinear
      polarization approximation of Eq. (<ref>) with
      q=k. But the sum of  and  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. <ref>.
      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. <ref> 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
       and  diagram expressions so that the sum in
      Fig. <ref> vanishes locally in the limit k∥ p_2.
      
      Later we will construct the modified vertex integrand, Eq. (<ref>) and modified self energy integrand, Eq. (<ref>), which exactly satisfy the desired Ward identity,
      
          η_2^ν k_μ/(k-η_2)^2 - η_2^2 V_ mod^μ +
          η_2^ν k_μ/(k-η_2)^2 - η_2^2 S_ mod^μ
          = 0,   when  k ∥ p_2,
      
      as illustrated in
      Fig. <ref>.
      
      
      
        §.§.§ 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. (<ref>).
      
      Each integral in the  and  diagrams of
      Fig. <ref>  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.  The second, B, refers to the
      loop momentum l, and can take the values B=1,S,H only for the singular regions of these diagrams.
      
      For all leading singular regions of the form (A,B) where neither A nor B
      equals H,  there are no hard sub-integrals by definition.   In these regions,
      (1,1), (2,1), (S,1),
      on-shell lines  are attached to tree rather than one-loop subgraphs and the
      necessary Ward identities, based as always on Eq. (<ref>), work point-by-point in loop momentum space.  
      The singularities factorize locally, and  do not require
      any alterations of  the functional form of the integrand.   
      We can thus focus on “single limits", (1,H), (2,H), (S,H),
      (H,1), where only one of the loop momenta becomes soft or collinear.
      We will now show how to treat the single-collinear limits.  The region (S,H) may be treated similarly to (1,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. (<ref>).
      
      
        *  In the same region (1, H), power-like singularities coming from doubled
        propagators in  diagrams, shown in
        Fig. <ref>(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  vertex diagram on the p_1 line or is adjacent to the  self energy diagram, as in Fig. <ref>.
      
      
        *   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. <ref>, 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  diagrams.
        In practice this means that the modification of the  diagram, which
        appears on the left-hand side of Fig. <ref>, should not change its singular
        behavior in the l ∥ p_1 limit.
      
      
        §.§.§ Modification of  integrands
      
      
      For  diagrams in Fig. <ref>(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 
      [.35]type2v_RH   ≡    
      V_0^μ(l,k,p_1) 
             
       =-i e^3 [  (  p_1 +  k ) γ_α  (  p_1 +  k +  l )  γ^μ (  p_1  + l )γ^α/(p_1+k)^2 (p_1+k+l)^2 l^2 (p_1+l)^2]  u(p_1)  .
      
      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 γ_αγ^μγ^νγ^ργ^α = 2ϵγ^μγ^νγ^ρ - 2 γ^ργ^νγ^μ, we find 
       
       V_0^μ(l,k,p_1)    =     2   i   e^3 /(p_1+k)^2 { (1-ϵ)  (  2(l^μ + p_1^μ) /l^2 (p_1+l)^2[ 1 - k^2 + ( p_1 +  l)   k /(p_1+k+l)^2 ]
       
                
      - 2(p_1^μ + l^μ)/(p_1+l)^2 (p_1+k+l)^2  -   2(p_1+k)^μ - γ^μ k /l^2 (p_1+k+l)^2  )
      
                
      -   ( p_1 +  k ) [ -( p_1 +  l) γ^μ k   +  ϵ k γ^μ ( p_1 +  l)  ] /l^2 (p_1+l)^2 (p_1+k+l)^2}  u(p_1)  .
      
        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. (<ref>) 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. (<ref>), 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. (<ref>) 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. (<ref>), 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. (<ref>) applied to both regular diagrams and / 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 (<ref>) is 
       
       V_ mod^μ(l,k,p_1)    ≡      2i e^3 /(p_1+k)^2 {(1-ϵ) /l^2(p_1+l)^2 [ p_1^μ  -  (p_1^μ + l^μ) (k^2 + ( p_1  +  l)  k )/(p_1+k+l)^2
               + l^μ ( k^2 -  l  k )/(k-l)^2  ] +  (1-ϵ)k^μ/(p_1+l)^2(p_1+k+l)^2
                 -       (1- ϵ) ( p_1 +  k) γ^μ/l^2   (p_1+k+l)^2  
         -    ( p_1 +  k)  [-  ( p_1 +  l)  γ^μ k   + ϵ k  γ^μ ( p_1 +  l) ] /l^2  (p_1+l)^2  (p_1+k+l)^2 }_⟨ T⟩  u(p_1)  .
         
        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  acting on the Dirac spinor, or by an explicit factor of k^2.   We emphasize that although V_ mod^μ 
        in (<ref>) 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  integrands
      
      
      The  self energy integrand is treated in a manner analogous to
      the  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. <ref>, 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,
      
          [.4]type2s_RH≡ S_0^μ(l,p_1,k)  
             = -i e^3   (  p_1+  k )/ ( p_1 +k  )^2  γ^α (  p_1 +  k +  l ) γ_α/l^2 (p_1+k+l)^2  (  p_1 + k )/ ( p_1 +k  )^2  γ^μ u(p_1)
          
             ,
      
      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,
      
          ∫ d^d l  S_0(l,p_1,k)    = ∫ d^d l  1/2[ S_0(l,p_1,k) + S_0(-p_1-k-l,p_1,k) ] 
              = i e^3   (1-ϵ) (_1 +  )/(p_1+k)^2 γ^μ  u(p_1) ∫ d^dl 1/l^2(p_1+l+k)^2  ,
      
      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 (<ref>) and (<ref>).    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 (<ref>), which is a simple combination of two terms, related by a change of variables,
        
      S_1^μ(l,p_1,k)    ≡     -e^3 i (1-ϵ) / ( p_1 +k  )^2   [
        k -  l /(p_1+l)^2 (k-l)^2 - _1 +  -  l  / l^2 (p_1+k+l)^2   ] γ^μ u(p_1)  .
       
      
      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 (<ref>).    It is straightforward to check that this expression, together with the modified vertex expression, Eq. (<ref>), obeys the modified Ward identity in region (2,H) as in Fig. <ref>.
      
      As it stands, however, Eq. (<ref>) is not quite what we want, because it introduces a new problem of loop polarizations, previously encountered in the unmodified  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. (<ref>) in terms of
      components parallel and transverse to the incoming momenta p_1 and
      p_2,
      
          l^μ = l · p_2/p_1 · p_2  p_1^μ
          +l · p_1/p_1 · p_2  p_2^μ
          +l_T^μ  .
      
      We will perform a  symmetrization only for
      the p_2 light-cone component in the two numerators  of 
      Eq. (<ref>).   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̅_+^μ and l̅_-^μ, defined by
      l̅_±^μ = l· p_2/p_1· p_2  p_1^μ ± 1/2 k· p_1/p_1· p_2  p_2^μ  +   l_T^μ
      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. (<ref>) does not change the
      result after integration over l. 
      The resulting expression is, 
      
      S_ mod^μ(l,p_1,k)    ≡     -ie^3  (1-ϵ) / ( p_1 +k  )^2   [
        k -  l_+ /(p_1+l)^2 (k-l)^2 - _1 +  -  l_- / l^2 (p_1+k+l)^2   ]_⟨ T ⟩γ^μ u(p_1)   .
       
      
      We are again using ⟨ T ⟩ to denote the transverse part
      of the global symmetrization Eq. (<ref>). 
      This is our final form for the  self energy diagram of Fig.<ref>(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. (<ref>).  This reflection is implicitly applied to
      every term of the amplitude integrand. 
      Under this transverse reflection and averaging, the denominators in 
      Eq. (<ref>) 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 . 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. (<ref>), 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
        single-collinear 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. (<ref>).  
      Then we manipulated the integrand in the following manner.
      
       
        *  For  diagrams, i.e. vertex corrections adjacent to
        external fermion lines as illustrated in
        Fig. <ref>(c)-(d), we modified the integrand expression
        V_0^μ in Eq. (<ref>) into V^μ_ mod in
        Eq. (<ref>) which is identical after integration. 
      
      
        *  For  diagrams, i.e. self energy corrections adjacent to
        external fermion lines as illustrated in
        Fig. <ref>(a)-(b), we modified the integrand expression
        S_0^μ in Eq. (<ref>) into S^μ_ mod in
        Eq. (<ref>). 
      
        *  For analogous  and  diagrams with self energy or vertex corrections adjacent to positron lines, the corresponding modifications are obtained by reflection.
      
        *  A four-fold symmetrization of the integrand (including all diagrams, some of which have been modified) is performed as in Eq. (<ref>).
      
      The outcome is a modified integrand, which is a sum of three contributions, as in Eq. (<ref>), now written as
      
          M^(2)(k,l) = [  ℳ^(2)_ + ℳ^(2)_]_ modified + ℳ^(2)_regular  ,
      
      with the modified self energy and vertex integrands specified by Eqs. (<ref>) and (<ref>), respectively.
      The symmetrizations of Eq. (<ref>) 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. (<ref>) and (<ref>) 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. (<ref>), 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. <ref>.  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. (<ref>), is identical to the original form in this region, so we are just exhibiting how this works.
        Applying the Dirac equation and p_1^2=0 to the averaged vertex, Eq. (<ref>), 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
        lim_l → -(1-x)p_1[ (p_1+l)^2 l^2  V_ mod^μ  ]     =    2 i e^3 /(p_1+k)^2 { p_1^μ -  x  (k^2 + x_1) p_1^μ/k^2+2xp_1· k
               -  x  (k^2 +  (1-x) _1) p_1^μ/k^2+2 (1-x) p_1· k   +  (_1 + ) (x _1) γ^μ/ k^2 + 2xp_1· k } u(p_1)
        
              =  -2i e^3 /(p_1+k)^2 { x
          ( _1  )/( k+x p_1)^2}γ^μ  u(p_1)  ,
        
        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 _1γ^μ  u(p_1) = - _1 γ^μ  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,
      
          lim_l → -(1-x)p_1[ (p_1+l)^2 l^2  V_ mod^μ  ]    ≡lim_l → -(1-x)p_1[ (p_1+l)^2 l^2  V_0^μ  ] 
             = -i e^3 2x/1- x [  1/(x_1 + ) - 1/ ( _1 +  ) ] γ^μ u(p_1) .
      
        This Ward identity is the expected one, and is represented in Fig. <ref>.
        
        
      
      
        §.§.§ Region (2,H) and the k Ward identity
      
      
      The final step is to verify that the Ward identity for the exchanged
      photon, Fig. <ref>, is respected by our
      modified integrand, Eq. (<ref>), once it is added to
      the appropriately-modified self energy, Eq. (<ref>).
      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
        
         k_μ V_ mod^μ(l,k,p_1)_k^2=0    =     -i e^3  (1-ϵ)/(p_1+k)^2  [  (_1 - l) / (p_1+k+l)^2 l^2   +  l /(k-l)^2 (p_1+l)^2  ]_⟨ T ⟩u(p_1) .
         
      
      To this expression, we will add the self energy integrand.  Neglecting
      terms that vanish at k^2=0, we find 
      
          k_μ V_ mod^μ(l,k,p_1) |_k^2=0 
          +k_μ S_ mod^μ(l,k,p_1) |_k^2=0   = 
            -i e^3  (1-ϵ)/(p_1+k)^2 
                ×[ 
           ( l_- -  l)  k/ (p_1+k+l)^2 l^2   +
             ( l- l_+ ) /(k-l)^2 (p_1+l)^2  ]_⟨ T ⟩ u(p_1)  .
      
      Given that 
      
          ( l- l_± ) p_2  = 0,
      
      we verify, 
      
          k_μ V_ mod^μ(l,k,p_1) |_k^2=0 
          +k_μ S_ mod^μ(l,k,p_1) |_k^2=0   = 0,    k ∥ p_2   ,
      
      and thus the factorization is complete on a point-by-point basis,  
      as illustrated in Fig. <ref>.
      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. <ref>, 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. (<ref>) and (<ref>), as in Figs. <ref> and <ref>, which we
      reproduce here as equations for ease of discussion,
      
          ℱ^(1)[𝐏_1 T 𝐏_1 ] = [.20]c19-5-alt
      
          ℱ^(2)[
          𝐏_1 T 𝐏_1
            ]
               =   [.2]19a 
          + [.2]19b 
          + [.2]19c
                 + [.2]19d 
          + [.2]19e 
          + [.2]19f .
       
      The loop momenta in the form factor amplitude are assigned and symmetrized with an identical procedure as that described in Subsection <ref> for the generic e^-(p_1)+ e^+(p_2) →γ^*(q_1)+… +
      γ^*(q_n) amplitudes.  
      For concreteness, we have made specific assignments on the  k,l
      labels to the virtual photons exchanged among fermion lines in the form factor
      diagrams. The averaging over loop momentum routings and
      transverse loop momentum reflections of Eq. (<ref>) is 
      understood implicitly.  
       self energy and  vertex subgraphs in the last four diagrams
      of Eq. (<ref>) are modified with the extra symmetrizations
      explained in subsection <ref>, just as the diagrams in Fig. <ref>.
      
      The squared vertex in each figure represents a matrix in spinor space, in between
      a pair of 𝐏_1≡ p_1  p_2/2p_1
        · p_2  projectors, Eq. (<ref>). 
      T captures contributions from amplitude propagators which remain 
      off-shell in infrared limits where either k and/or l vanish. 
      As such, and due to factorization, it is independent of the loop momenta
      assigned to the virtual photons outside the squared vertex in the above diagrams.  
      For our purposes, T will stand for  a tree amplitude integrand,
      independent of both k and l,
      inserted in the one-loop and two-loop form factor of 
      Eqs. (<ref>) and (<ref>), or a one-loop
      amplitude integrand, depending on the loop momentum l, inserted in
      the one-loop form factor of Eq. (<ref>). 
      
      We have explained in Section <ref> how an infrared finite remainder of the one-loop
      amplitude integrand is derived with the use of a one-loop form factor
      subtraction at the level of the integrand. We recall here the final result  at one loop of
      Eqs. (<ref>), equivalent to Eqs. (<ref>) and (<ref>),
      
          M^(1)_ IR-finite =
             M^(1) - ℱ^(1)[ 𝐏_1ℳ^(0)𝐏_1].
      
      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 two-loop amplitude.  
      Starting from the modified two-loop integrand ℳ^(2)
      constructed in Subsection <ref>, 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. (<ref>),
      
          _double-IRℳ^(2) = ℱ^(2)[ 𝐏_1ℳ^(0)𝐏_1]   ,
      
      where we use the two-loop form factor as a local  approximation of all
      double-infrared singularities.   
      The subtraction gives
      
          ℳ^(2)_double-IR-finite = ℳ^(2) -
            ℱ^(2)[  𝐏_1ℳ^(0)𝐏_1]  ,
      
      in which we recognize the beginning of the pattern of the two-loop finite remainder, Eq. (<ref>), so far with a subtraction
      only for the double-IR limits.   The modification of  and  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 (single-soft) 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  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,  
      
          _single-IR ℳ^(2)_double-IR-finite  
          = _single-IR(ℳ^(2) - ℱ^(2)[  𝐏_1ℳ^(0)𝐏_1] ).
      
      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. <ref> 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. (<ref>) and the one-loop form factor integrand,
      
          _single-IRℳ^(2) 
           = ℱ^(1)[ 𝐏_1ℳ^(1)𝐏_1] .
      
      Similarly, we remove the single-IR limit from the two-loop form factor
      counterterm, 
      
          _single-IR ℱ^(2)[  𝐏_1ℳ^(0)𝐏_1] 
          =ℱ^(1)[ 𝐏_1ℱ^(1)
              [ 𝐏_1ℳ^(0)𝐏_1] 
              𝐏_1] .
      
      The above compact notation corresponds to  
      the following momentum representation
      
          ℱ^(1)[ 𝐏_1ℱ^(1)
              [ 𝐏_1ℳ^(0)𝐏_1] 
              𝐏_1] = -e^4/4v̅ (p_2) γ^μ1/ k - p_2
                × 𝐏_1γ^ν1/ l - p_2𝐏_1ℳ^(0)𝐏_11/ l + p_1γ_ν 𝐏_1
                ×1/ k + p_1γ_μ u(p_1)
          + .
       
      Combining, Eq. (<ref>) and Eq. (<ref>) into
      Eq. (<ref>) we have that
      
          _single-IR ℳ^(2)_double-IR-finite   = _single-IR(ℳ^(2) - ℱ^(2)[  𝐏_1ℳ^(0)𝐏_1] ) 
             = ℱ^(1)[ 𝐏_1ℳ^(1)𝐏_1]                                                                                -
               ℱ^(1)[ 𝐏_1ℱ^(1) [ 𝐏_1ℳ^(0)𝐏_1] 𝐏_1] 
             = ℱ^(1)[ 𝐏_1(
           ℳ^(1) - ℱ^(1)
              [ 𝐏_1ℳ^(0)𝐏_1]  )
              𝐏_1] 
             = ℱ^(1)[ 𝐏_1ℳ^(1)_IR-finite𝐏_1]   ,
      
      where M^(1)_ IR-finite is given by Eqs. (<ref>) using (<ref>), removing
      the external spinors.    This
       provides the second IR subtraction in Eq. (<ref>).
      
      We recall from the discussion at the end of Sec. <ref> that the 𝐏_1 projectors play an important role in
      Eq. (<ref>) in preventing the introduction of spurious
      IR singularities. Let us examine their action here, in the momentum representation of the 
      above counterterm, which reads
      
          ℱ^(1)[ 𝐏_1ℳ^(1)_IR-finite𝐏_1]  
          = -ie^2/4v̅ (p_2) γ^μ1/ l - p_2 𝐏_1ℳ^(1)_IR-finite(k) 
          𝐏_11/ l + p_1γ_μ u(p_1)
          
           [2mm]
                + .
       
      In a collinear limit such as k ∥ p_1 of
      Eq. (<ref>),  the factor ℳ^(1)_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 
        ℳ^(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 𝐏_1 = 𝐏_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,
      
          ℳ^(2)_IR-finite = ℳ^(2) 
            - ℱ^(2)[  𝐏_1ℳ^(0)𝐏_1]
              - ℱ^(1)[  𝐏_1ℳ^(1)_IR-finite𝐏_1] .
      
      Effectively, this is Eq. (<ref>) before
      subtractions for UV behavior are carried out.
       Eq. (<ref>) 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^(2)_ IR-finite defined by Eq. (<ref>)
      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. (<ref>).
      As above, we denote them by the action of an operator _ UV, so that, schematically,  
      
          _ UV M^(2)_ IR-finite   =   _ UV M^(2)
          - _ UVℱ^(2)[  𝐏_1ℳ^(0)𝐏_1]
          
          
          
          - _ UVℱ^(1)[  𝐏_1ℳ^(1)_IR-finite𝐏_1]  .
      
      The first term in Eq. (<ref>) removes the ultraviolet
      singularities of the two-loop amplitude integrand M^(2)
          M^(2)-_ UV M^(2)
          = 
            M^(2)_ UV-finite .
       
      Similarly, the second term in Eq. (<ref>) cancels the
      ultraviolet singularities in the two-loop form factor integrals of ℱ^(2)  which serves as  a counterterm of double infrared
      singularities, 
      
          ℱ^(2)[  𝐏_1ℳ^(0)𝐏_1]
          -_ UVℱ^(2)[  𝐏_1ℳ^(0)𝐏_1]
          = ℱ^(2)_ UV-finite[  𝐏_1ℳ^(0)𝐏_1].
      
      The last term in Eq. (<ref>) must remove ultraviolet
      singularities from both the one-loop form factor integral of
      ℱ^(1) and its one-loop infrared-subtracted amplitude
      kernel 𝐏_1ℳ^(1)_IR-finite𝐏_1,  
      
          ℱ^(1)[  𝐏_1ℳ^(1)_IR-finite𝐏_1] 
          - 
          _ UVℱ^(1)[  𝐏_1ℳ^(1)_IR-finite𝐏_1]
          =ℱ^(1)_ UV-finite[  𝐏_1ℳ^(1)_finite𝐏_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 ℱ^(1) and two-loop ℱ^(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  
      
          M^(2) =  M^(2)_ singular 
          +   M^(2)_ finite ,
      
      where  all ultraviolet and infrared singularities reside in the first
      term, M^(2)_ 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. (<ref>), reproduced here,
      
          ℳ^(2)_finite = ℳ^(2)_UV-finite 
            - ℱ^(2)_UV-finite[  𝐏_1ℳ^(0)𝐏_1]
              - ℱ^(1)_UV-finite[  𝐏_1ℳ^(1)_finite𝐏_1],
      
      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^(2)_ singular  is a sum of  ultraviolet finite
       form factor one-loop and two-loop amplitudes whose analytic
       integration involves well-known 
      integrals <cit.> 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 <ref>, we explained in detail the derivation of local
      ultraviolet counterterms for the diagrams in M^(1) and 
      ℱ^(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 
      ℱ^(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  and  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 <ref>, we derived one-loop local ultraviolet
      counterterms for the one-loop QED vertex and electron self energy graphs, Eqs. (<ref>) 
      and (<ref>), 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 <ref> we
      introduced local ultraviolet counterterms for the one-loop electron
      propagator and the one-loop electron-electron-photon vertex that
      respect the QED Ward identity. Let us examine the function of the
      counterterms of  Eq. (<ref>) and (<ref>), for
      diagrams in M^(2)_ regular such as the ones illustrated 
      in Fig. <ref>.
      
      The corresponding counterterm diagrams,  in which the ultraviolet
      limit is approximated by Eqs. (<ref>) and (<ref>),
      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. <ref>.
      
      
      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 ℱ^(1)_UV-finite[  𝐏_1ℳ^(1)_finite𝐏_1] in
      Eq. (<ref>). We remark that for this
      cancellation to occur, it is important that  the local vertex of the
      form factor counterterm, 𝐏_1ℳ^(1)_finite𝐏_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^(2)_ regular.
      
      
      
        §.§.§ One-loop ultraviolet counterterms specific to    and   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  and diagrams, which we have cast in a modified form.  
      
      In the large l limit, the modified vertex expression of Eq. (<ref>) is logarithmically divergent, so we can easily write down its UV subdivergence counterterm by taking the leading large-l limit of the expression,
      
          _l →∞ V^μ_ mod. =  2i e^3   (1-ϵ)/(p_1+k)^2 (l^2-M^2)^2[
              p_1^μ + k^μ - ( p_1 +  k) γ^μ  -  2 l^μ l  k /l^2 - M^2]_⟨ T ⟩ u(p_1)   .
      
      Though the original  electron self energy expression
      Eq. (<ref>) is linearly divergent in the limit of
      large l, our equivalent form  of Eq. (<ref>) is only
      logarithmically divergent.  This allows us to easily write down the UV counterterm
      
          _l →∞ S^μ_ mod.   = e^3 i (1-ϵ)/(p_1+k)^2(l^2-M^2)^2[
           p_1 +  l̅_+ - l̅_- 
          + (l̅_+ + l̅_-) 2k· l/(l^2-M^2)]_⟨ T ⟩ γ^μ u(p_1) .
      
      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. (<ref>) and
      (<ref>) satisfy the necessary Ward identity for
      factorization (or rather, cancellation) of the singularity 
      in the limit k ∥ p_2, analogous to Eq. (<ref>),
      
          k_μ_l →∞ S_ mod.^μ  |_k → yp_2 = -  k_μ_l →∞ V_ mod.^μ  |_k → yp_2.
      
      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. (<ref>) and (<ref>) 
      are necessary.
      
        and 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. (<ref>).  
      Consistently, we use the expressions of 
      Eq. (<ref>)  and Eq. (<ref>) 
      as ultraviolet counterterms in the form factor F^(2) as well. 
      
      
      
        §.§.§ One-loop ultraviolet counterterms specific to the form factor
        F^(1) and  F^(2)  amplitudes
      
      
      In subsections <ref> and <ref>, we
      constructed ultraviolet  counterterms for one-loop vertex and
      self energy graphs in regular diagrams as well as  and   
      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. (<ref>)
      
      The ultraviolet singularity of  F^(1) 
       (similar to the usual QED one-loop
      vertex) has already been  encountered in Eq. (<ref>). The same singularity appears in 
      F^(2), for which 
      a sample diagram is shown of Fig. <ref>. 
      The large-l limit in Fig. <ref> is singular
      and the singularity can be also subtracted by a counterterm that is
      identical to 
      the square bracket of Eq. (<ref>).
      
      
      Next we construct a one-loop ultraviolet counterterm for  the last term in
      Eq. (<ref>) with a momentum representation [As pointed out earlier, the k ↔ l symmetrization in Eq. (<ref>) makes it rather unimportant how we assign k and l to the two photons.]
          ℱ^(1)[  𝐏_1ℳ^(1)_IR-finite𝐏_1] =
            v̅(p_2) [ 
              (-i e^2 ) γ^μ ( k -  p_2)
                𝐏_1
                ℳ^(1)_IR-finite(l)
                𝐏_1
                ( k +  p_1) γ_μ/k^2 (k+p_1)^2 (k-p_2)^2]  u(p_1).
      
      The subtraction of the one-loop singularity as l becomes large,  
      by using the counterterms of  Eq. (<ref>) and
      Eq. (<ref>), 
      turns ℳ^(1)_IR-finite(l) 
      into ℳ^(1)_finite(l) of Eq. (<ref>).
      A further  counterterm for the one-loop singularity as k
      becomes large is required and it is similar to
      Eq. (<ref>). 
      The ultraviolet approximation of k →∞, following  the l →∞
      subtraction, reads
      
          _k →∞ℱ^(1)[  𝐏_1ℳ^(1)_finite(l) 𝐏_1] =
            v̅(p_2) [
              (-i e^2 ) γ^μk𝐏_1
                ℳ^(1)_finite(l)
                𝐏_1
                kγ_μ/(k^2-M^2)^3]  u(p_1).
      
      
        §.§.§ 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. (<ref>).  Counterterms for these singularities 
      can be constructed straightforwardly, with the following steps. 
      
        *  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.
      
        *  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. (<ref>) and (<ref>). 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 <ref>, the
      UV counterterms described above enable us to construct the fully-subtracted
      two-loop amplitude given in Eq. (<ref>) or equivalently (<ref>).   
      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. (<ref>).
      
      
      
      
      § NUMERICAL CHECKS OF CANCELLATION OF SINGULARITIES FOR OFF-SHELL 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_ finite^(2) in the
      full expression of Eq. (<ref>).   
      
      We generate the Feynman diagrams using the program QGRAF
       <cit.>. 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. (<ref>).
      
      
        *   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. (<ref>).
      
      
        *  The diagrams for M^(2)_4, which include light-by-light scattering, have UV subtractions, as 
      described in Sec. <ref>, 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  and  subdiagrams, specified by
      Eqs. (<ref>) and (<ref>), respectively.   
      
      
        *  The  and  modified integrands are also employed
      for the form factor subtractions, as described in Sec. <ref>.
      
      
        *  Ultraviolet counterterms that ensure convergence for two-loop photonic amplitudes and
      their IR subtractions are constructed as in Sec. <ref>.
      
      
      
      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, 
      
          p_1^μ   = (1, 0, 0, 1), 
          
            p_2^μ   = (1, 0, 0, -1), 
          
            q_1^μ   = (1, 0, 1/3, 1/7), 
          
            q_2^μ   = (1, 0, -1/3, -1/7) ,
      
      with polarization vectors (satisfying the transversality condition ϵ_i^* · q_i = 0)
      
          ϵ_1^*  μ   = (4, -4, 9, 7), 
          ϵ_2^*  μ   = (2, -1, -3, -7),
      
      and external spinors in the Weyl basis (satisfying the Dirac equations p_1 u(p_1) = v̅(p_2)  p_2 = 0)
      
          u(p_1) = [  0; -7;  3;  0 ],   v̅(p_2) = [  0, 13, -9,   0 ]  .
      
      
      A basis of vectors transverse to the p_1-p_2 plane is defined as
      
          n_x^μ = (0, 1, 0, 0),      n_y^μ = (0,0, 1, 0)   .
      
      
      We also use a parameterization of  the loop momenta
      k^μ, l^μ, which depends on a small parameter δ and
      scaling exponents ω_k,   l^i, i={+,-,T}, starting
      from an arbitrary choice of loop momenta at δ =1,
      
          k^μ   = 33/17δ^ω_k^+  p_1^μ - 48/89δ^ω_k^-  p_2^μ +  δ^ω_k^T( 21/23  n_x^μ + 21/41  n_y^μ), 
          
            l^μ   = 47/23δ^ω_l^+  p_1^μ - 7/61δ^ω_l^-  p_2^μ + δ^ω_l^T( -37/73  n_x^μ -39/67  n_y^μ)   .
      
      By selecting
      appropriately the values of the exponents we can approach particular IR / UV
      regions. 
      The proof of  UV finiteness by testing convergence as loop momenta are scaled to infinity is well established  <cit.>.   The IR scalings of Eq. (<ref>) shown in  in Tables <ref> - <ref> serve the analogous purpose of testing the behavior of integrals in the neighborhoods of regions where loop momenta are pinched between coalescing poles.  The demonstration of IR finiteness or divergence for wide-angle scattering and production amplitudes by use of these scalings is discussed in Refs. <cit.>.   The explicit tests here demonstrate the adequacy of our IR counterterms in regulating singular behavior in these amplitudes.[For other processes, typically involving momentum pinches in so-called Glauber regions <cit.>, additional scaling tests become necessary to verify finiteness.]
      The mass regulator for the UV counterterms is arbitrarily chosen as M=5/3.
      
      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^-,ω_k^T,ω_l^+,ω_l^-,ω_l^T
      which parameterize it appropriately <cit.>.
      For example, in the single-soft limit where k is soft and
      l is hard,  we set
      ω_k^+ = ω_k^- = ω_k^T =1,  ω_l^+ = ω_l^- =
      ω_l^T = 0 in Eq. (<ref>). 
      Now we can perform an analytic Laurent expansion in δ. 
      For example, in the same region, we obtain 
      
          . v̅(p_2)  M_ finite^(2)  u(p_1) 
          |_t-channel = 1/δ^3( -630.517 - 24.0761 i ) + 𝒪( 1/δ^2)   .
      
      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 𝒪(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
      <ref>.
      
      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 <ref>,
      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.
      
      
      
       §.§ Fermion loop contributions
      
      Table <ref> shows the relevant IR / UV scaling
      behavior of M^(2)_2 which receives contributions 
       from Feynman diagrams with a vacuum polarization correction to an  
      internal photon propagator. Again, only the t-channel ordering of the two
      final-state photons is presented.
      
      
      Table <ref> 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.
      
      
      § 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 process-dependent
      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 point-like 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. (<ref>) 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. (<ref>) 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 <cit.> and in Ref. <cit.> where amplitudes, from the class of
      processes with equivalent subtractions for pinch singularities as described here, 
      were evaluated efficiently at one loop.  Ref. <cit.> 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. <cit.> 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.  
      
      
      
      
      § ACKNOWLEDGMENTS
      
      
      We are grateful to Alexander Ochirov, Alexander Penin, Andrea Pelloni, Armin Schweitzer, Ben Ruijl,
      Bernhard Mistlberger, Claude Duhr, Dario Kermanschah, Valentin Hirschi, Vittorio del Duca,
      Yao Ma, Zeno Capatti and Zoltan Kunszt for inspiring discussions.  
      This research was supported in part by the Pauli Centre for
      Theoretical Physics, the National Science Foundation
      under Grants PHY-1620628 and PHY-1915093,  by the Swiss National Science Foundation
      under contract SNF200021_179016 and by the European Commission through the ERC grant “pertQCD”.
      
      JHEP
  