Solution (source code)

= Solution

Let $\mathbf A_i=\sum_{j\ne i}K_{ij}(\mathbf v_j-\mathbf v_i)$. Conservation of <probability> in the full <phase space> gives the <dissipative Liouville equation>
$$
\boxed{\partial_t\varrho+\sum_{i=1}^N\nabla_{\mathbf x_i}\cdot(\mathbf v_i\varrho)+\sum_{i=1}^N\nabla_{\mathbf v_i}\cdot(\mathbf A_i\varrho)=0.}
$$
Normalize $\int\varrho\,\prod_i d^3x_i\,d^3v_i=1$. Decay at infinity makes all boundary fluxes vanish, so the equation preserves this <normalization of a probability distribution>.

Unlike the <Hamiltonian Liouville equation>, this transport has nonzero <phase-space divergence>. In three <dimensions>,
$$
\nabla_{\mathbf v_i}\cdot\mathbf A_i=-3\sum_{j\ne i}K_{ij},\qquad
\partial_t\varrho+\sum_i\left(\mathbf v_i\cdot\nabla_{\mathbf x_i}\varrho+\mathbf A_i\cdot\nabla_{\mathbf v_i}\varrho\right)=3\varrho\sum_i\sum_{j\ne i}K_{ij}.
$$
The <probability density> increases along a contracting <characteristic flow map>; setting its material derivative to zero would lose conservation of <probability>. The $j=i$ force is identically zero, including its derivative, which is why the last sums exclude it.