Gluon TMDs and NRQCD matrix elements in $J/\psi$ production at an EIC

In this paper we analyze azimuthal asymmetries in the processes of unpolarized and polarized $J/\psi \,(\Upsilon)$ production at an Electron-Ion Collider. Apart from giving access to various unknown gluon transverse momentum distributions, we suggest to use them as a new method to extract specific color-octet NRQCD long-distance matrix elements, i.e.\ $\langle0|{\cal O}_{8}^{J/\psi}(^{1}S_{0})|0\rangle$ and $\langle0|{\cal O}_{8}^{J/\psi}(^{3}P_{0})|0\rangle$, whose values are still quite uncertain and for which lattice calculations are unavailable. The new method is based on combining measurements of analogous asymmetries in open heavy-quark pair production which can be performed at the same energy. To enhance the gluon contribution one can consider smaller values of $x$ and, in order to assess the impact of small-$x$ evolution, we perform a numerical study using the MV model as a starting input and evolve it with the JIMWLK equations.