  We present a calculation of the $cos (2 \phi_h)$ asymmetry in $J/\psi$
production in electron-proton collision at the future electron-ion collider
(EIC), a useful channel to probe the gluon TMDs. We calculate the asymmetry at
next-to-leading order (NLO) in $\alpha_s$ in the framework of generalized
factorization. The dominating sub-process is $\gamma^* +g \rightarrow
J/\psi+g$. The production of $J/\psi$ is calculated in the non-relativistic
QCD(NRQCD) framework with the inclusion of both color singlet and color octet
contributions. Numerical estimates of the $cos(2\phi_h)$ asymmetry are given in
the kinematical region to be accessed by the future EIC. The asymmetry depends
on the parameterization of the gluon TMDs, as well as on the long distance
matrix elements (LDMEs). We use both Gaussian-type parameterization and
McLerran-Venugopalan model for the TMDs in the kinematical region of small-$x$,
where the gluons play a dominant role. We obtain sizable asymmetry, which may
be useful to probe the ratio of the linearly polarized and the unpolarized
gluon distribution in the proton.
