  We consider the transverse momentum spectrum and the $\cos 2\phi$ azimuthal
distribution of $J/\psi$ mesons produced in semi-inclusive, deep-inelastic
electron-proton scattering, where the electron and the proton are unpolarized.
At low transverse momentum, we propose factorized expressions in terms of
transverse momentum dependent gluon distributions and shape functions. We show
that our formulae, at the order $\alpha_s,$ correctly match with the collinear
factorization results at high transverse momentum. The latter are computed at
the order $\alpha_s^2$ in the framework of nonrelativistic QCD (NRQCD), with
the inclusion of the intermediate $^3S_1^{[1]}$ color-singlet Fock state, as
well as the subleading color-octet ones that are relatively suppressed by a
factor $v^4$ in the NRQCD velocity parameter $v$. We show that the
$^1\!S_0^{[8]}$ and $^3\!P_J^{[8]}$ ($J = 0,1,2$) contributions diverge in the
small transverse momentum region and allow us to determine the perturbative
tails of the shape functions, which carry the same quantum numbers. These turn
out to be identical, except for the overall magnitude given by the appropriate
NRQCD long distance matrix element.
