$\mathcal{O}(M)$, where $M$ is the number of mesh points in real space used in the simulation. In this respect our approach also differs from the original work by [@Widrow1993], who stick to the classical Poisson form of gravity.

We stress that this step is not formally necessary. The well-known techniques used to solve Poisson’s equation on a mesh (FFT, Gauss-Seidel relaxation, etc.) can also be used in our framework. This change of equation for gravity does just make the computations slightly faster in the cases where our approximation is valid. However, in the case of cosmological simulations with vastly different scales interacting, it is unclear how the Klein-Gordon equation for gravity would behave and defaulting to standard mesh techniques might be required.

The algorithm in brief {#ssec:algorithm}
======================

Here we present the main algorithm of our framework, decomposed in a few simple steps. A formal derivation and a discussion of the convergence and accuracy of the method will be presented in the next section.

Choose the parameters of your simulation. The precision and speed of the method is governed by three parameters $\hbar$, $c$ and $N$. The algorithm of choosing them is the following:\
The parameters $c$ and $\hbar$ are 