  We compute the quark and gluon transverse momentum dependent parton
distribution functions at next-to-next-to-next-to-leading order (N$^3$LO) in
perturbative QCD. Our calculation is based on an expansion of the differential
Higgs boson and Drell-Yan production cross sections about their collinear
limit. This method allows us to employ cutting edge techniques for the
computation of cross sections to extract the universal building blocks in
question. The corresponding perturbative matching kernels for all channels are
expressed in terms of simple harmonic polylogarithms up to weight five. As a
byproduct, we confirm a previous computation of the soft function for
transverse momentum factorization at N$^3$LO. Our results are the last missing
ingredient to extend the $q_T$ subtraction methods to N$^3$LO and to obtain
resummed $q_T$ spectra at N$^3$LL$^\prime$ accuracy both for gluon as well as
for quark initiated processes.
