  We compute the top quark contribution to the two-loop amplitude for on-shell
$Z$ boson pair production in gluon fusion, $gg \to ZZ$. Exact dependence on the
top quark mass is retained. For each phase space point the integral reduction
is performed numerically and the master integrals are evaluated using the
auxiliary mass flow method, allowing fast computation of the amplitude with
very high precision.
