In this proceeding we present our derivation of a compact representation of the renormalized $N$-jettiness soft function free of infrared and collinear divergences through next-to-next-to-leading order (NNLO) in perturbative QCD. The number $N$ of hard partons enters as a parameter in the formula of the finite remainder. We demonstrate analytically the cancellation of all infrared and collinear singularities between the bare soft function and its renormalization matrix in color space and we compare our results with the ones available in the literature.