  Traditional computations of the dark matter (DM) relic abundance, for models
where attractive self-interactions are mediated by light force-carriers and
bound states exist, rely on the solution of a coupled system of classical
on-shell Boltzmann equations. This idealized description misses important
thermal effects caused by the tight coupling among force-carriers and other
charged relativistic species. We develop for the first time a comprehensive
ab-initio derivation for the description of DM long-range interactions in the
presence of a hot and dense plasma background directly from non-equilibrium
quantum field theory. Most importantly, the scattering and bound states get
strongly mixed in the thermal plasma environment, representing a characteristic
difference from a pure vacuum theory computation. The main result of this work
is a novel differential equation for the DM number density, written down in a
form which is manifestly independent under the choice of what one would
interpret as a bound or a scattering state at finite temperature. The collision
term, unifying the description of annihilation and bound state decay, turns out
to have in general a non-quadratic dependence on the DM number density. This
generalizes the form of the conventional Lee-Weinberg equation which is
typically adopted to describe the freeze-out process. We prove that our general
number density equation is consistent with previous literature results under
certain limits.
