  We consider transverse momentum dependent gluon distributions inside both
unpolarized and transversely polarized protons and show how they can be probed
by looking at azimuthal modulations in $e \, p \to e \, J/\psi \, \mathrm{jet}
\, X$. We find that the contribution due to quark induced subprocesses is
always suppressed in the considered kinematic regions, accessible in principle
at a future Electron-Ion Collider. Our model-independent estimates of the
maximal values of these asymmetries allowed by positivity bounds suggest the
feasibility of their measurement. In addition, by adopting the
McLerran-Venugopalan model for the unpolarized and linearly polarized gluon
densities, we study the behavior of the $\cos2\phi$ asymmetries in the
small-$x$ limit.
