The far field potential of an infinite checkerboard with squares of alternate surface charge density I am trying to find the potential far away from an infinite checkerboard having squares of size $x^2$ and carrying alternate surface charge $\pm \sigma$ (for example, all black squares carry $+\sigma$ and all white squares carry $-\sigma$).
I have tried to parametrise the charge distribution by expanding it into a Fourier series (of a square wave in both $x$ and $y$ directions on the plane) and replacing it into the formula for the potential $V = \frac{1}{4\pi\epsilon_0}\int \frac{\rho(\vec{r'})}{|\vec{r} - \vec{r'}|}d^3r'$. But I cannot integrate the resulting expression and I am very confused if this is the correct way to do it or not.  
 A: I think a possible way to find the solution is the following.


*

*The checkerboard function (proportional to $\sigma$) can be written as$$
C(x,y)=\frac{16}{\pi^2} \sum_{n,m \,\,\mbox{odd}} \frac{1}{n m}\sin \left(\ \frac{\pi n x}{L} \right)
\sin \left(  \frac{\pi m y}{L} \right)
$$

*You should use the symmetries of your system to deduce that the potential is an even function of $z$, if you take a coordinate system with $x$ and $y$ on the checkerboard. As the discontinuity of the normal electric field crossing the checkerboard is proportional to $\sigma$, then $$
\left. \frac{\partial \phi}{\partial z}(x,y,z) \right|_{z=0^+} = \frac{\sigma}{2 \epsilon_0}
$$ must be proportional to che checkerboard function. 

*By direct check the function $$e^{-\beta z} \sin p x \sin q y$$ is a solution of the equation $\nabla^2 \phi_{pq}=0$ if $\beta^2=p^2+q^2$. Now you can construct a series of terms which solve the Poisson equation for $z>0$ by writing$$\phi = \frac{16}{\pi^2} \sum_{n,m \,\,\mbox{odd}} \frac{K_{mn}}{n m}\sin \left(\ \frac{\pi n x}{L} \right)
\sin \left(  \frac{\pi n y}{L} \right) \exp \left(-\frac{\pi\sqrt{n^2+m^2}}{L}z  \right)$$
where $K_{mn}$ are arbitrary coefficients.

*The derivative of this equation at $z=0$ must be proportional to the checkerboard function. This fixes the coefficients and we obtain the final result$$\phi = -\frac{8}{\pi^2} \frac{\sigma L}{\epsilon_0 \pi} \sum_{n,m \,\,\mbox{odd}} \frac{1}{n m\sqrt{m^2+n^2}}\sin \left(\ \frac{\pi n x}{L} \right)
\sin \left(  \frac{\pi n y}{L} \right) \exp \left(-\frac{\pi\sqrt{n^2+m^2}}{L}z  \right)$$ 
For large values of $z$ the dominant term is for $m=n=1$ (details of the checkerboard are "blurred"). Note that the field decreases faster than any power law. This depend from the fact that there is not a good multipole approximation that could work: two squares give a dipole, but four give an octupole canceling the dipole and so on.
$$\phi \sim -\frac{8}{\pi^2} \frac{\sigma L}{\epsilon_0 \pi}  \frac{1}{\sqrt{2}}\sin \left(\ \frac{\pi x}{L} \right)
\sin \left(  \frac{\pi y}{L} \right) \exp \left(-\frac{\pi\sqrt{2}}{L}z  \right)$$ 
n.b.: the checkerboard size is $L$.
