Abstract: We present a complete set of analytic helicity amplitudes for top quark pair production via gluon fusion at two-loops in QCD. For the first time, we include corrections due to massive fermion loops which give rise to integrals over elliptic curves. We present the results of the missing master integrals needed to compute the amplitude and obtain an analytic form for the finite remainders in terms of iterated integrals using rationalised kinematics and finite field sampling. We also study the numerical evaluation of the iterated integrals.