All (un)polarized transverse momentum dependent functions (TMDs), both distribution and fragmentation functions, are defined with the same universal soft function, which cancels spurious rapidity divergences within an individual TMD and renders them well-defined hadronic quantities. Moreover, it is independent of the kinematics, whether it is Drell-Yan, deep inelastic scattering, or ${e}^{+}{e}^{\ensuremath{-}}\ensuremath{\rightarrow}2$ hadrons. In this paper, we provide this soft function at next-to-next-to-leading order (NNLO), necessary for the calculation of all TMDs at the same order, and to perform the resummation of large logarithms at next-to-next-to-next-to-leading-logarithmic accuracy. From the results we obtain the $D$ function at NNLO, which governs the evolution of all TMDs. This work represents the first independent and direct calculation of this quantity. Given the all-order relation through a Casimir scaling between the soft function relevant for gluon TMDs and the one for quark TMDs, we also obtain the first at NNLO. The used regularization method to deal with the rapidity divergences is discussed as well.