  Cosmic rays (CRs) at sub-TeV energies play a fundamental role in the chemical
and dynamical evolution of molecular clouds, as they control the ionisation,
dissociation, and excitation of H$_{2}$. Their characterisation is important
both for the interpretation of observations and for the development of
theoretical models. The methods used so far for estimating the CR ionisation
rate ($\zeta$) in molecular clouds have several limitations due to
uncertainties in the adopted chemical networks. We refine and extend the method
proposed by Bialy (2020) to estimate $\zeta$ by observing rovibrational
transitions of H$_{2}$ at near-infrared wavelengths, which are mainly excited
by secondary CR electrons. Combining models of interstellar CR propagation and
attenuation with the calculation of the expected secondary electron spectrum
and updated H$_{2}$ excitation cross sections by electron collisions, we derive
the intensity of the four H$_{2}$ rovibrational transitions observable in
dense, cold gas: (1-0)O(2), (1-0)Q(2), (1-0)S(0), and (1-0)O(4). The proposed
method allows the estimation of $\zeta$ for a given observed line intensity and
H$_{2}$ column density. We are also able to deduce the shape of the low-energy
CR proton spectrum impinging upon the molecular cloud. We present a look-up
plot and a web-based application that can be used to constrain the low-energy
spectral slope of the interstellar CR proton spectrum. We comment on the
capability of the James Webb Space Telescope to detect these near-infrared
H$_{2}$ lines, making it possible to derive for the first time spatial
variation of $\zeta$ in dense gas. Besides the implications for the
interpretation of the chemical-dynamic evolution of a molecular cloud, it will
be possible to test competing models of CR propagation and attenuation in the
interstellar medium, as well as compare CR spectra in different Galactic
regions.
