In this paper, we study the azimuthal asymmetry in the J/ψ leptoproduction in unpolarized ep collisions. There are two independent azimuthal asymmetry modulations, namely cos(ψ) and cos(2ψ), where ψ 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). In addition, this process can be used to test the magnitude of the 1 S [8] 0 long-distance matrix element, which is extraordinarily large from various phenomenological extraction. 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 cos(ψ) modulation, while the 1 S [8] 0 dominance picture can be tested by measuring the values of the cos(2ψ) modulation. Calculating their values and carrying out a meticulous statistical analysis, we find that at an integrated luminosity L = 1000pb −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/ψ production mechanism might be discovered.