Numerically solving 2D poisson equation by FFT, proper units The 2D Poisson equation  is:
(1)$$\frac{d^2\varphi(x,y)}{dx^2}+\frac{d^2\varphi(x,y)}{dy^2}=-\frac{\varrho(x,y)}{\epsilon_0\epsilon}$$
And in $k$-space it is in form of:
(2)$$(k_x^2+k_y^2) \phi(k_x,k_y)=-\frac{\rho(k_x,k_y)}{\epsilon_0\epsilon}.$$
To solve numerically I use complex FFT (FFTW library in C). For area of physical size L and grid of size N (grid const. h=L/N), discrete coordinates and periodic boundaries I have:
(3)$$\rho[k_x,k_y]=\frac{1}{N^2}\sum_{x=0}^{N-1}\sum_{y=0}^{N-1}\varrho[x,y]e^{-j2\pi(xk_x+yk_y)/N}$$
I can multiply both sides by $-\frac{1}{\epsilon\epsilon_0}$, then dividing $\rho[k_x,k_y]$ by $(k_x^2+k_y^2)$ at every point (taking into account that FFT output is symmetric over $k_{x}$ = N/2 and $k_y$ = N/2) I get $\phi[k_x,k_y]$. By inverse FFT:
(4)$$\varphi[x,y]=\sum_{k_x=0}^{N-1}\sum_{k_y=0}^{N-1}\phi[k_x,k_y]e^{j2\pi(xk_x+yk_y)/N}$$
I get 2D potential coming from density of charge $\varrho$.
But I am missing some scaling factor in above definitions and I can not figured this out. In short what to do to match the SI units in above definitions so they fit, both in terms of Fourier transform and real result of $\varphi$.
My main concern is that if I input the test charge of density $\varrho(x,y)=\frac{e}{h^2}\delta[x,y]$ in units of $C/m^2$  in the middle of area with $L=10^{-6}m$  I expect output similar to Coulomb potential $\frac{1}{r}\frac{e}{4\pi\epsilon\epsilon_0}$, but the solution differs in many orders. Why?
So I think it may be a problem with dividing by k. What should be k for this definition of discrete transform? For example $k=2\pi h/L$ or $k=2\pi x/N$ 
 A: Ok I think I solved the problem. So to divide the density FFT by $k^2$ I need actual values of $k$ in $k$-space for my system. FFTW orders the result of transformation in so called "in-order" output, that means in first quadrant of FFT the first pixel corresponds to DC frequency and next to $k/L$ frequency ($k$ is from $0$ to $N-1$) where $L$ is length of whole system. Smallest wave length of the system is $h=L/N$, then $1/h$ is highest frequency. So $k_x$ in above equations should be changed to $2 \pi*k/L$. The coordinate $k$ should also depend on which quadrant of output fft we are, for symmetric quadrant we should use $N-k$ instead of $k$. Probably that is it.
A: Substitute
$$
\varphi(x,y) = \int dx dy \phi(k_x, k_y) e^{i k_x x + i k_y y},~~\varrho(x,y) = \int dx dy \rho(k_x, k_y) e^{i k_x x + i k_y y}
$$
in the first equation and this should immediately give the equation you desire. Also, just for potential future purposes, note
$$
\int dx e^{i k x } = 2\pi \delta (k)
$$
(See here for discussion.)
