system of equations but would allow to simulate the full 6D phase space and hence allow simulation of alternative cosmologies, such as the one including WDM or free-streaming neutrino contributions. The principal difficulty is then to find a good mapping between the distribution function $f(x,v)$ of interest and its “quantum” equivalent $\psi(x)$ and vice-versa. This is achieved by using the so-called *Wigner distribution function*

$$f(x,v) \simeq  \int e^{\frac{i}{\hbar}vy}\psi^*\left(x+\frac{y}{2}\right)\psi\left(x-\frac{y}{2}\right)
d^3y.$$

which obeys an equation similar to the Vlasov equation but is constructed from wave functions. The main feature of this mapping is that the density field can simply be expressed as

$$\label{eq:1}
  \rho(x) = |\psi(x)|^2.$$

However, the limitation of this approach is that one single wavefunction is in general not sufficient to encode all the complexity of the distribution function and we would then use the more general version:

$$\label{eq:2}
\rho(x) = \sum_n |\psi_n(x)|^2.$$

The summation index $n$ can, as a first thought, be understood as a sum over the velocities $v$ that appear in the distribution function $f(x,v)$. We somehow trade a 6D function for a (finite) set of 3D (complex valued) functions. We will, however, demonstrate that the number of wavefunctions required can be very low (of order unity 