  We solve coupled momentum-dependent Boltzmann equations for the phase space
distribution of cosmic relic particles, without resorting to approximations of
assuming kinetic equilibrium or neglecting backscattering or elastic
interactions. Our method is amendable to precision numerical computations. To
test it, we consider two benchmark models where the momentum dependence of dark
matter distribution function is potentially important: a real singlet scalar
extension near the Higgs resonance and a sterile neutrino dark matter model
with a singlet scalar mediator. The singlet scalar example shows that the
kinetic equilibrium may hold surprisingly well even near sharp resonances.
However, the integrated method may underestimate the relic density by up to 40%
in extreme cases. In the sterile neutrino dark matter model, we studied how the
inclusion of previously ignored elastic interactions and processes with initial
state sterile neutrinos could affect the nonthermal nature of their resulting
distributions. Here the effects turned out to be negligible, proving the
robustness of the earlier predictions.
