  The emission of $\gamma$-rays after a neutron capture in a cryogenic detector
can generate mono-energetic nuclear recoils in the sub-keV regime, of direct
interest for the calibration of Dark Matter and Coherent Elastic Neutrino
Nucleus Scattering experiments. Here we show that accurate predictions of the
nuclear recoil spectra induced by neutron captures require taking into account
the interplay between the development in time of the de-excitation
$\gamma$-cascade of the target nucleus and that of the associated atomic
collisions in matter. We present detailed simulations coupling the FIFRELIN
code for the description of the $\gamma$-cascades and the IRADINA code for the
modelling of the fast atomic movements in matter. Nuclear recoil spectra are
predicted, and made available to the community, for concrete cases of
Al$_2$O$_3$, Si, Ge and CaWO$_4$ crystals exposed to a low intensity beam of
thermal neutrons. We find that timing effects cause new calibration peaks to
emerge in the recoil spectra and also impact the shape of the continuous recoil
distribution. We discuss how they could give access to a rich physics program,
spanning the accurate study of the response of cryogenic detectors in the
sub-keV range, tests of solid state physics simulations and tests of nuclear
models.
