  In this paper, we study the azimuthal asymmetry in the $J/\psi$
leptoproduction in unpolarized $ep$ collisions. There are two independent
azimuthal asymmetry modulations, namely $\mathrm{cos}(\psi)$ and
$\mathrm{cos}(2\psi)$, where $\psi$ is the azimuthal angle of the lepton
scattering plane with respect to the hadron-interacting plane. We calculate the
two modulations as functions of four kinematic variables, and find that they
provide a very good laboratory to distinguish several models describing the
heavy quarkonium production, including the color-singlet (CS) model, the
nonrelativistic QCD (NRQCD) associated with the $^1S_0^{[8]}$ dominance
picture, and the NRQCD in which the values of all the three color-octet (CO)
long-distance matrix elements are of the same order. In order to make definite
conclusions, we restrict our calculation in a specific kinematic region, where
the CS and CO mechanisms can be distinguished by scrutinizing the values of the
$\mathrm{cos}(\psi)$ modulation, while the $^1S_0^{[8]}$ dominance picture can
be tested by measuring the values of the $\mathrm{cos}(2\psi)$ modulation.
Calculating their values and carrying out a meticulous statistical analysis, we
find that at an integrated luminosity $\mathcal{L}=1000\mathrm{pb}^{-1}$, the
statistical uncertainties of the two quantities are small enough to tell the
three models apart. When this experiment is implemented at the future $ep$
colliders such as the EIC, crucial information for the $J/\psi$ production
mechanism might be discovered.
