  We compute the total cross-section and invariant mass distribution for
heavy-quark pair production in $e^+e^-$ annihilation at the
next-to-next-to-next-to-leading order in QCD. The obtained results are
expressed as piecewise functions defined by several deeply expanded power
series, facilitating a rapid numerical evaluation. Utilizing top-pair
production at a collision energy of 500 GeV as a benchmark, we observe a
correction of approximately $0.1\%$ for the total cross-section and around
$10\%$ for the majority of the invariant mass distribution range. These results
play a crucial role in significantly reducing theoretical uncertainty: the
scale dependence has been diminished to $0.06\%$ for the total cross-section
and to $5\%$ for the invariant mass distribution. This reduction of uncertainty
meets the stringent requirements of future lepton colliders.
