We present the two-loop virtual amplitudes for heavy-quark pair production in light quark-antiquark annihilation and gluon fusion channels, including full spin and color dependence. We use expansions around kinematical limits and numerical integration to obtain results for the involved master integrals. From these, we determine the renormalised infrared finite remainders of the coefficients of amplitude decompositions in terms of color and spin structures. The remainders are given in form of numerical interpolation grids supported by expansions around the production threshold and the high energy limit. Finally, we provide the spin density matrix, which encodes the heavy-quark spin correlations and is sufficient for phenomenological applications. Our results are necessary for the derivation of top-quark pair production cross sections in hadron collisions in the narrow width approximation with next-to-next-to-leading order accuracy in QCD.