  Linear cosmological observables can be used to probe elastic scattering of
dark matter (DM) with baryons. Availability of high-precision data requires a
critical reassessment of any assumptions that may impact the accuracy of
constraints. The standard formalism for constraining DM-baryon scattering
pre-recombination is based on assuming a Maxwell-Boltzmann (MB) velocity
distribution for DM. This assumption is not always justified and does not allow
for probing DM self-interactions in addition to its interactions with baryons.
Lifting the MB assumption requires solving the full collisional Boltzmann
equation (CBE), which is highly non-trivial. Earlier work proposed a more
tractable Fokker-Planck (FP) approximation to the CBE, but its accuracy
remained unknown. In this work, we numerically solve the exact CBE for the
first time, in a homogeneous expanding background. We consider DM-baryon
scattering cross-sections that are positive power-laws of relative velocity. We
derive analytical expressions for the collision operator in the case of
isotropic differential scattering cross-sections. We then solve the background
CBE numerically and use our solution for the DM velocity distribution to
compute the DM-baryon heat-exchange rate, which we compare against those
obtained with the MB assumption and FP approximation. Over a broad range of
DM-to-baryon mass ratios, we find that the FP approximation leads to a maximum
error of 17%, significantly better than the up to 160% error introduced by the
MB assumption. While our results strictly apply only to the background
evolution, the accuracy of the FP approximation is likely to carry over to
perturbations. This motivates its implementation into cosmological Boltzmann
codes, where it can supersede the much less accurate MB assumption, and allow
for a more general exploration of DM interactions with baryons and with itself.
