  In most direct detection experiments, the free nuclear recoil description of
dark matter scattering breaks down for masses $\lesssim$ 100 MeV, or when the
recoil energy is comparable to a few times the typical phonon energy. For dark
matter lighter than 1 MeV, scattering via excitation of a single phonon
dominates and has been computed previously, but for the intermediate mass range
or higher detector thresholds, multiphonon processes dominate. We perform the
first calculation of the scattering rate via multiphonon production for the
entire keV-GeV dark matter mass range, assuming a harmonic crystal target. We
provide an analytic description that connects the single phonon, multiphonon,
and the nuclear recoil regimes. Our results are implemented in the public
package $\texttt{DarkELF}$.
