  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.
