  We present a new and fully general algorithm for the automated construction
of the integrands of two-loop scattering amplitudes. This is achieved through a
generalisation of the open-loops method to two loops. The core of the algorithm
consists of a numerical recursion, where the various building blocks of
two-loop diagrams are connected to each other through process-independent
operations that depend only on the Feynman rules of the model at hand. This
recursion is implemented in terms of tensor coefficients that encode the
polynomial dependence of loop numerators on the two independent loop momenta.
The resulting coefficients are ready to be combined with corresponding tensor
integrals to form scattering probability densities at two loops. To optimise
CPU efficiency we have compared several algorithmic options identifying one
that outperforms naive solutions by two orders of magnitude. This new algorithm
is implemented in the OpenLoops framework in a fully automated way for two-loop
QED and QCD corrections to any Standard Model process. The technical
performance is discussed in detail for several $2\to2$ and $2\to 3$ processes
with up to order $10^5$ two-loop diagrams. We find that the CPU cost scales
linearly with the number of two-loop diagrams and is comparable to the cost of
corresponding real-virtual ingredients in a NNLO calculation. This new
algorithm constitutes a key building block for the construction of an automated
generator of scattering amplitudes at two loops.
