Average value divergence in spectral method for Poisson equation

Question

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.)

Answer 1

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.