  In this paper we calculate analytically the perturbative matching
coefficients for unpolarized quark and gluon Transverse-Momentum-Dependent
(TMD) Parton Distribution Functions (PDFs) and Fragmentation Functions (FFs)
through Next-to-Next-to-Next-to-Leading Order (N$^3$LO) in QCD. The N$^3$LO TMD
PDFs are calculated by solving a system of differential equation of Feynman and
phase space integrals. The TMD FFs are obtained by analytic continuation from
space-like quantities to time-like quantities, taking into account the
probability interpretation of TMD PDFs and FFs properly. The coefficient
functions for TMD FFs exhibit double logarithmic enhancement at small momentum
fraction $z$. We resum such logarithmic terms to the third order in the
expansion of $\alpha_s$. Our results constitute important ingredients for
precision determination of TMD PDFs and FFs in current and future experiments.
