# Average value divergence in spectral method for Poisson equation

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:

1. Top ground plane: $$\quad\tilde{\phi}_1\left(\mathbf{k},L_z\right)=0$$
2. Bottom ground plane: $$\quad\tilde{\phi}_2\left(\mathbf{k},-L_z\right)=0$$
3. Continuity of potential at $$z=0$$: $$\quad\tilde{\phi}_1\left(\mathbf{k},0\right)=\tilde{\phi}_2\left(\mathbf{k},0\right)$$
4. 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.)

I have an answer, but I'm still curious if there are simpler answers.

The basic idea is that we break the charge density into two parts: one with an average density of zero, and one with a constant value. We use the above procedure on the part with zero average charge density, and we know the answer for a constant charge density. Then we just add the two potentials together.

The surface charge density $$\sigma\left(x,y\right)$$ can be broken into two parts:

$$\sigma\left(x,y\right) = \left(\sigma\left(x,y\right) - \bar{\sigma}\right) + \bar{\sigma},$$

where $$\bar{\sigma}$$ is the average value of $$\sigma\left(x,y\right)$$ over the unit cell.

First, deal with the term in parenthesis. Define $$\sigma_1\left(x,y\right) = \sigma\left(x,y\right) - \bar{\sigma}$$. Now, $$\sigma_1$$ has zero average value, so $$\tilde{\sigma}_1\left(\mathbf{0}\right)=0$$. That makes it fair to set $$\tilde{\phi}_1\left(\mathbf{0},z\right)=0$$, and we can solve for $$\phi_1$$ using the above method.

Now, $$\sigma_2\left(x,y\right) = \bar{\sigma}$$ is just a uniform sheet charge, and we know that it will produce $$\mathbf{E}=\pm \frac{\bar{\sigma}}{2\epsilon_0}$$ (+ for $$z>0$$, - for $$z<0$$). This results in $$\phi_2=\frac{\bar{\sigma}}{2\epsilon_0}\left(\pm L_x-z\right)$$.

Then we put everything together to find $$\phi = \phi_1 + \phi_2$$.

I'd be interested in any other approaches, because I'm interested in tackling more complicated problems (multiple charged sheets and dielectric layers) with a similar method, and extending the above method will certainly work, but maybe it's not the simplest answer.