We present a recent estimate [8] of azimuthal asymmetry $(cos~2\phi)$ in $J/\psi$ production in ep collision within the kinematics of the planned electron-ion collider (EIC). We find a small but sizable asymmetry in the small-$x$ region which can be observed at the EIC. This asymmetry directly probes the Weizs\"acker-Williams (WW) type linearly polarized gluon distribution in the unpolarized proton. We calculate the asymmetry at next-to-leading order in the kinematical region which satisfies, energy fraction of $J/\psi$, $z < 1$ and the subprocess $\gamma^∗ + g \rightarrow c + \bar{c}+ g$ gives the leading contribution. We use color singlet (CS) model for $J/\psi$ production. We estimate the asymmetry in two models for TMDs: (i) Gaussian type parameterization and (ii) McLerran-Venugopalan (MV) model at small-$x$ and also we present the upper bound of the asymmetry.