  We set up a formalism for calculating transverse-momentum-dependent parton
distribution functions (TMDs) using the tools of saturation physics. By
generalizing the quasi-classical Glauber-Gribov-Mueller/McLerran-Venugopalan
approximation to allow for the possibility of spin-orbit coupling, we show how
any TMD can be calculated in the saturation framework. This can also be applied
to the TMDs of a proton by modeling it as a large "nucleus." To illustrate our
technique, we calculate the quark TMDs of an unpolarized nucleus at large-x:
the unpolarized quark distribution and the quark Boer-Mulders distribution. We
observe that spin-orbit coupling leads to mixing between different TMDs of the
nucleus and of the nucleons. We then consider the evolution of TMDs: at
large-x, in the double-logarithmic approximation, we obtain the Sudakov form
factor. At small-x the evolution of unpolarized-target quark TMDs is governed
by BK/JIMWLK evolution, while the small-x evolution of polarized-target quark
TMDs appears to be dominated by the QCD Reggeon.
