<p>The transverse-momentum-dependent (TMD) soft function for non-relativistic heavy quark pair production at hadron colliders is analytically computed at next-to-next-to-leading order (NNLO) in the strong coupling expansion. We present the details of our computational approach and analyze the general two-loop structure of the soft function. The final result, which takes a particularly simple form, provides the last missing ingredient for a complete NNLO calculation of color-octet <i>S</i>-wave quarkonium hadroproduction — including charmonium, bottomonium, and toponium — using the <i>q</i>_<i>T</i>-slicing formalism. It also enables next-to-next-to-next-to-leading-logarithmic (N³LL) resummation at small transverse momentum for the same process.</p>