  We show that the gauge-invariant transverse-momentum-dependent (TMD) quark
distribution function can be expressed as a sum of all higher-twist collinear
parton matrix elements in terms of a transport operator. From such a general
expression, we derive the nuclear broadening of the transverse momentum
distribution. Under the maximal two-gluon correlation approximation, in which
all higher-twist nuclear multiple-parton correlations with the leading nuclear
enhancement are given by products of twist-two nucleon parton distributions, we
find the nuclear transverse momentum distribution as a convolution of a
Gaussian distribution and the nucleon TMD quark distribution. The width of the
Gaussian, or the mean total transverse momentum broadening squared, is given by
the path integral of the quark transport parameter $\hat q_F$ which can also be
expressed in a gauge invariant form and is given by the gluon distribution
density in the nuclear medium. We further show that contributions from
higher-twist nucleon gluon distributions can be resummed under the extended
adjoint two-gluon correlation approximation and the nuclear transverse momentum
distribution can be expressed in terms of a transverse scale dependent quark
transport parameter or gluon distribution density. We extend the study to hot
medium and compare to dipole model approximation and ${\cal N}=4$
Supersymmetric Yang-Mills (SYM) theory in the strong coupling limit. We find
that multiple gluon correlations become important in the strongly coupled
system such as ${\cal N}=4$ SYM plasma.
