Not yet reviewed by Pith; the record is open.

We propose to use transverse momentum $p_T$ distribution of $J/\psi$ production at the future Electron Ion Collider (EIC) to explore the production mechanism of heavy quarkonia in high energy collisions. We apply QCD and QED collinear factorization to the production of a $c\bar{c}$ pair at high $p_T$, and non-relativistic QCD factorization to the hadronization of the pair to a $J/\psi$. We evaluate $J/\psi$ $p_T$-distribution at both leading and next-to-leading order in strong coupling, and show that production rates for various color-spin channels of a $c\bar{c}$ pair in electron-hadron collisions are very different from that in hadron-hadron collisions, which provides a strong discriminative power to determine various transition rates for the pair to become a $J/\psi$. We predict that the $J/\psi$ produced in electron-hadron collisions is likely unpolarized, and the production is an ideal probe for gluon distribution of colliding hadron (or nucleus). We find that the $J/\psi$ production is dominated by the color-octet channel, providing an excellent probe to explore the gluon medium in large nuclei at the EIC.