  We describe the unitarity approach for the numerical computation of two-loop
integral coefficients of scattering amplitudes. It is well known that the
leading propagator singularities of an amplitude's integrand are related to
products of tree amplitudes. At two loops, Feynman diagrams with doubled
propagators appear naturally, which lead to subleading pole contributions. In
general, it is not known how these contributions can be directly expressed in
terms of a product of on-shell tree amplitudes. We present a universal
algorithm to extract these subleading pole terms by releasing some of the
on-shell conditions. We demonstrate the new approach by numerically computing
two-loop four-gluon integral coefficients.
