  We present a systematic study of deuterated molecular hydrogen (HD) at high
redshift, detected in absorption in the spectra of quasars. We present four new
identifications of HD lines associated with known $\rm H_2$-bearing Damped
Lyman-$\alpha$ systems. In addition, we measure upper limits on the $\rm HD$
column density in twelve recently identified $\rm H_2$-bearing DLAs. We find
that the new $\rm HD$ detections have similar $N({\rm HD})/N(\rm H_2)$ ratios
as previously found, further strengthening a marked difference with
measurements through the Galaxy. This is likely due to differences in physical
conditions and metallicity between the local and the high-redshift interstellar
media. Using the measured $N({\rm HD})/N({\rm H_2})$ ratios together with
priors on the UV flux ($\chi$) and number densities ($n$), obtained from
analysis of $\rm H_2$ and associated CI lines, we are able to constrain the
cosmic-ray ionization rate (CRIR, $\zeta$) for the new $\rm HD$ detections and
for eight known HD-bearing systems where priors on $n$ and $\chi$ are
available. We find significant dispersion in $\zeta$, from a few $\times
10^{-18}$ s$^{-1}$ to a few $\times 10^{-15}$ s$^{-1}$. We also find that
$\zeta$ strongly correlates with $\chi$ -- showing almost quadratic dependence,
slightly correlates with $Z$, and does not correlate with $n$, which probably
reflects a physical connection between cosmic rays and star-forming regions.
