use a particle method in which the distribution function is sampled by a finite number $N$ of particles such that

$$f(x,v) \cong \frac{1}{N} \sum_{i=1}^{N} m_i \delta\left(x - x_i\right)\delta\left(v - v_i\right).
\label{eq:nbody}$$

Each particle or body is then evolved according to Newton’s law under the influence of the gravitational potential created by all the others as described by Poisson’s equation. In other words, [$N$-body ]{}simulations solve the Vlasov equation via its characteristics by sampling the initial phase space distribution with a discrete number of particles. The number of bodies is typically chosen as large as computationally feasible. The [$N$-body ]{}formalism is thus a Monte-Carlo approximation of the Vlasov-Poisson system. The advent of large supercomputers combined with the development of more efficient numerical algorithms has enabled the field of cosmological simulations to make considerable progress over the last decades. Simulations such as the *Millennium run* [@Springel2005] or *Bolshoi simulation* [@Klypin2011] are able to follow as many as a few billion particles.

The complicated part of the [$N$-body ]{}simulation is the evaluation of the forces between pairs of particles. Over the years, many ingenious techniques (see [@Dehnen2011] for a review) have been invented to reduce the algorithms complexity for the force integration to 