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$