is way too large to track individually each of them on a computer. Therefore most of the cosmological [$N$-body ]{}simulations use macroscopically large simulation “particles” (with their masses ranging from masses much larger than DM particles up to the size of a small galaxy, $10^8-10^9 {M_\odot}$).

The problem of dark matter evolution in the Universe can be formulated as an evolution of a collisionless self-gravitating fluid.

The main tool used to describe this dark matter fluid is the *phase space density distribution* $f(x,v,t)$, defined such that $f(x,v,t)dqdv$ represents the mass of material at position $x$ moving at velocity $v$ at time $t$. This function is usually normalized such that its integral over all positions and velocities gives the total mass

$$\int d^3x\int d^3v f(x,v,t) = M_{\rmn{tot}}.$$

Notice that one could also normalize this integral to one or to the total number of particles in the system. When integrating over velocity space only, one gets the usual mass density $\rho(x)$, whereas integrating over all space returns the velocity distribution $d_v(v)$:

$$\rho(x) = \int  f(x,v,t)d^3v, \qquad d_v(v) = \int  f(x,v,t) d^3x.$$

This distribution function obeys the Liouville theorem [@Binney2008] and if the only force acting on the particle is the gravitational potential 