  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.
