simple technique to generate the initial wave function. They use a set of particles sampling the phase space distribution function exactly as in the case of [$N$-body ]{}simulations. They can then turn each particle into a Gaussian in phase space by smoothing it and use this set of wave packets as their initial wave function.

In our approach, we depart from this need of an initial [$N$-body ]{}sampling by considering other techniques to generate the set of wave functions. By doing so, we allow for a completely generic distribution function and should, in principle, not experience the consequences of an *a priori* artificial Monte-Carlo sampling of $f(x,v)$.

The second feature of our framework is the replacement of the Poisson equation by a Klein-Gordon equation for the potential $U(x)$: $$-\frac{1}{c^2}\frac{\partial^2}{\partial t^2}U + \nabla^2_xU = 4\pi G\delta\rho,$$

where $c$ is the numerical speed of gravity. This scalar gravity equation is, once again, purely a mathematical trick to reduce the complexity of the original system (\[eq:VP\]) and not an attempt to modify Newton’s gravity. Such a replacement makes the framework entirely local and does not require complicated integration methods for the Poisson equation. The complexity of the scheme is then formally reduced to 