Analytic NNLO transverse-momentum-dependent soft function for heavy quark pair hadroproduction at threshold

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 $S$-wave quarkonium hadroproduction--including charmonium, bottomonium, and toponium--using the $q_T$-slicing formalism. It also enables next-to-next-to-next-to-leading-logarithmic (N$^3$LL) resummation at small transverse momentum for the same process.