We consider the transverse momentum spectrum and the cos 2φ azimuthal distribution of J/ψ 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 α s , correctly match with the collinear factorization results at high transverse momentum. The latter are computed at the order α 2 s in the framework of nonrelativistic QCD (NRQCD), with the inclusion of the intermediate 3 S [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 [8] 0 and 3 P [8] J (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.