Abstract We study the polar and azimuthal decay angular distributions of J/ψ mesons produced in semi-inclusive, deep-inelastic electron-proton scattering. For the description of the quarkonium formation mechanism, we adopt the framework of nonrelativistic QCD, with the inclusion of the intermediate color-octet channels that are suppressed at most by a factor v4 in the velocity parameter v relative to the leading color-singlet channel. We put forward factorized expressions for the helicity structure functions in terms of transverse momentum dependent gluon distributions and shape functions, which are valid when the J/ψ transverse momentum is small with respect to the hard scale of the process. By requiring that such expressions correctly match with the collinear factorization results at high transverse momentum, we determine the perturbative tails of the shape functions and find them to be independent of the J/ψ polarization. In particular, we focus on the cos 2ϕ azimuthal decay asymmetry, which originates from the distribution of linearly polarized gluons inside an unpolarized proton. We therefore suggest a novel experiment for the extraction of this so-far unknown parton density that could be performed, in principle, at the future Electron-Ion Collider.