  Semi-inclusive $J/\psi$ production in electron-proton collisions is a
promising process to study gluon transverse momentum distributions (TMDs) at
the future Electron-Ion Collider. In this article, we improve on previous
studies of the $\cos 2 \phi$ azimuthal asymmetry that arises from the linear
polarization of gluons inside unpolarized protons by including TMD evolution.
We find that in the TMD regime the asymmetry grows monotonically with
increasing transverse momentum of the outgoing $J/\psi$, in contrast to tree
level calculations with Gaussian TMDs. Our predictions for the asymmetry at EIC
can become very large at larger $x$, $Q$, and transverse momenta, even larger
than the positivity bound. This problem stems from the very small $b$ region
and implies a range of validity of TMD factorization that is more restricted
than usually expected. We also include an estimate of the nonperturbative
uncertainty from the large $b$ region and we conclude that it is smaller than
the largest source of uncertainty, which stems from the choice of Color Octet
Long-Distance Matrix Elements.
