I'd like to know how to deal with a divergence when trying to solve the Poisson equation for electrostatics with a simple spectral method. I'm not sure how to best state my problem, so I'll explain how I got to the problem. I'm sorry that this makes the post long.
Suppose that I have two conducting ground planes at $z=\pm L_z$, an infinite domain in $x$ and $y$, and some surface charge at $z=0$ (i.e. $\rho\left(x,y,z\right)=\sigma\left(x,y\right)\delta\left(z\right)$). I can take Poisson's equation
$$\nabla^2\phi\left(x,y,z\right)=-\frac{\rho\left(x,y,z\right)}{\epsilon_0}=\frac{\sigma\left(x,y\right)\delta\left(z\right)}{\epsilon_0}$$
and Fourier transform both sides of the equation in $x$ and $y$ (but not $z$) to get
$$-k^2\tilde{\phi}\left(k_x, k_y, z\right) + \frac{\partial^2}{\partial z^2}\tilde{\phi}\left(k_x, k_y, z\right)=\frac{\tilde{\sigma}\left(k_x,k_y\right)\delta\left(z\right)}{\epsilon_0},$$
where $k^2=k_x^2+k_y^2$, and for brevity, I'll use $\mathbf{k}$ for $\left(k_x,k_y\right)$. When $z\neq0$, the above reduces to
$$-k^2\tilde{\phi}\left(\mathbf{k}, z\right) + \frac{\partial^2}{\partial z^2}\tilde{\phi}\left(\mathbf{k}, z\right)=0,$$
which has the solution of the form
$$\tilde{\phi}\left(\mathbf{k},z\right)=A e^{-kz} + B e^{kz}.$$
So, I can look for two solutions $\tilde{\phi}_1\left(\mathbf{k},z\right)=A_1 e^{-kz} + B_1 e^{kz}$ and $\tilde{\phi}_2\left(\mathbf{k},z\right)=A_2 e^{-kz} + B_2 e^{kz}$ for above and below $z=0$, respectively. I can solve for $A_{1,2},B_{1,2}$ using four boundary conditions:
- Top ground plane: $\quad\tilde{\phi}_1\left(\mathbf{k},L_z\right)=0$
- Bottom ground plane: $\quad\tilde{\phi}_2\left(\mathbf{k},-L_z\right)=0$
- Continuity of potential at $z=0$: $\quad\tilde{\phi}_1\left(\mathbf{k},0\right)=\tilde{\phi}_2\left(\mathbf{k},0\right)$
- Effect of surface charge at $z=0$: $\quad\left.\frac{d\tilde{\phi}_1}{dz}\right|_{z=0} - \left.\frac{d\tilde{\phi}_2}{dz}\right|_{z=0} = -\frac{\tilde{\sigma}\left(\mathbf{k}\right)}{\epsilon_0}$
For example, I find
$$A_1 = \frac{\tilde{\sigma}\left(\mathbf{k}\right)}{2k\epsilon_0\left(e^{-2kL_z} + 1\right)}.$$
Once I know $\tilde{\phi}\left(\mathbf{k},z\right)$, I can do some inverse Fourier transforms to get $\phi\left(x,y,z\right)$.
Suppose now that I want to do this for a system that is periodic in $x$ and $y$ (say, with a square domain), and I want to do transforms numerically. So, I discretize $\sigma\left(x,y\right)$ on a square mesh and replace the Fourier transforms with discrete Fourier transforms, and almost everything works. The only problem is the $\frac{1}{k}$ that show up in $A_{1,2},B_{1,2}$. It causes $\tilde{\phi}\left(\mathbf{k},z\right)$ to blow up at $\mathbf{k}=0$. I can just set $\tilde{\phi}\left(\mathbf{0},z\right)=0$, but that sets the average value of $\phi$ to zero, and I lose information I'd like to have.
Is there a way to get around that?
(I'll note that, regardless of the value of $\sigma\left(x,y\right)$, each unit cell has zero net charge due to the conducting ground planes, which will have a net charge equal and opposite to the net charge from $\sigma$. I'll also note that I can solve this problem numerically with FEM, and it produces good results, but I'd like to solve the problem with a lighter-weight method.)