To compute the evolution of the number density operator, let us apply the covariant form of the Liouville’s operator to the corresponding phase space distribution function:
where and C[f] is the collisional operator, which takes into account the processes which change the number of particles (like annihilations or decays).
In the case of a FRW Universe for which , we have
Integrating over phase space, we can relate this to the time evolution of the number density,
Regarding the collisional operator, let us concentrate on annihilation processes, where SM particles (A, B) can annihilate to form DM particles (1, 2) or vice-versa. The phase space corresponding to each particle is defined as
, from where
.
The terms account for viable phase space of the produced particles, taking into account whether they are fermions (-) or bosons (+). Assuming no CP violation in the DM sector (T invariance), . Also, energy conservation in the annihilation process allows us to write , thus,
Since the SM particles are in equilibrium, , where we have defined the thermally-averaged cross-section as
.
Hence, we are left with the familiar form of Boltzmann equation,
Notice that this is an equilibrium restoring equation. If the RHS dominates, then traces its equilibrium value . However, when , the RHS can be neglected, and the resulting differential equation implies that . This equivalent to saying that the DM particles do not annihilate anymore and their number density decreases only because the scale factor of the universe increases.