1
$\begingroup$

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

$\endgroup$
0
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.