$U(x)$, we can write a closed system of equations for the formation of structures [@Bertschinger1995; @Bernardeau2002]:

$$\begin{aligned}
    {\frac{\partial f}{\partial \tau}} + \frac{v}{a(\tau)}{\frac{\partial f}{\partial x}} - a(\tau)\nabla U {\frac{\partial f}{\partial v}} &= 0,  \\
    \nabla^2 U &= 4\pi Ga^2(\tau) \delta\rho,
  \end{aligned}
\label{eq:VP}$$

where $x$ and $v$ are comoving coordinates and velocities, $a(\tau)$ is the scale factor and $\tau$ is the conformal time (We will use this convention throughout this paper). This Vlasov-Poisson system has no solution in the general case and the only way to handle it is to use numerical techniques. For completeness, we also give the expressions for the density and density contrast:

$$\begin{aligned}
 \rho(x,\tau) &=& \frac{1}{a^3(\tau)}\int f(x,v,\tau)d^3v, \\
 \delta\rho(x,\tau) &=& \frac{1}{a^3(\tau)}\left(\int f(x,v,\tau)d^3v - \frac{M_{\rmn{tot}}}{V_{\rmn{tot}}} \right),\end{aligned}$$

where $V_{\rmn{tot}}$ is the total comoving volume over which we average.

Structure formation simulations {#sec:theory}
===============================

The numerical analysis of the Vlasov-Poisson system of equations (\[eq:VP\]) is very challenging. The first reason is that the system is six-dimensional. Recent simulations can only handle up to $64$ resolution elements in each spatial and velocity space direction [@Yoshikawa2013] due to memory restrictions. Even the use of the biggest supercomputers would not allow to go much beyond this figure.

The second shortcoming of such technique is 