I attempted to do a N-body-simulation and compare it to the analytical Zeldovich map until shell-crossing appears. In order to do so, I placed them uniformly with the Zeldovich initial condition, e.g they have some offset and are not entirely uniformly distributed. $$ \Delta x = \frac{L}{n_\text{grid}} $$ where L is the size of the box in an $n_\text{grid} \times n_\text{grid}$ mesh with the initial conditions $q$ $$ q = (q_x, q_y) \in [0,L)^2 $$ and their initial positions $$ x(q,a) = (q - a A (2 \pi) \sin(2xq)) \operatorname{mod} L $$ where $a$ is the cosmological time and $A=\frac{1}{4\pi^2}$. The velocity is then of course given by the derivative.
Afterwards, I "split" their mass with a classic cloud-in-cell approach, and solved the Poisson-equation with an FFT . Note that my particle density is literally the density and not some fancy power spectrum. $$ \delta_{ij} = \frac{p_{ij}}{\bar{p}}-1 $$ Afterwards, I applied my solved potential with a leapfrog integrator (conserves the Hamiltonian) to the inverse cloud-in-cell masses.
Now, for 1D, the plot of q vs x at some certain point I receive a sinusoidal wave for particle mesh as well as for Zeldovich, which is correct I think 
, but I would like to test this for 2D as well, e.g $x_1(q_1,a) = (q_1 - a A (2 \pi) \sin(2 \pi q_1),...$. Which curve should I expect from a plot of q vs x or y in 2D? Does it even make sense to compare these with each other? What should I expect from the plot of x vs y before shell crossing?
