  Low-energy cosmic-rays (CRs) control the thermo-chemical state and the
coupling between gas and magnetic fields in dense molecular clouds, the sites
of star-formation. However, current estimates of the low-energy CR spectrum ($E
\lesssim 1$ GeV) and the associated CR ionization rate are highly uncertain. We
apply, for the first time, a new method for constraining the CR ionization rate
and the CR spectral shape using H$_2$ rovibrational lines from cold molecular
clouds. Using the MMIRS instrument on the MMT, we obtained deep near-infrared
(NIR) spectra in six positions within four dense cores, G150, G157, G163, G198,
with column densities $N_{\rm H_2} \approx 10^{22}$ cm$^{-2}$. We derive
3$\sigma$ upper limits on the H$_2$ $(1-0)$S(0) line (2.22 $\mu$m) brightness
in the range $I = 5.9 \times 10^{-8}$ to $1.2 \times 10^{-7}$ erg cm$^{-2}$
s$^{-1}$ sr$^{-1}$ for the different targets. Using both an analytic model and
a numerical model of CR propagation, we convert these into upper limits on the
CR ionization rate in the clouds' interior, $\zeta = 1.5$ to $3.6 \times
10^{-16}$ s$^{-1}$, and lower limits on the low-energy spectral slope of
interstellar CR protons, $\alpha = -0.97$ to $-0.79$. We show that while MMT
was unable to detect the H$_2$ lines due to high atmospheric noise,
JWST/NIRSpec will be able to efficiently detect the CR-excited H$_2$ lines,
making it the ideal method for constraining the otherwise elusive low-energy
CRs, shedding light on the sources and propagation modes of CRs.
