Electric potential of a spheroidal gaussian I'm looking for results that compute the electrostatic potential due to a spheroidal gaussian distribution. Specifically, I'm looking for solutions of equations of the form
$$
\nabla^2\Phi=N\exp\left({-\rho^2/2\sigma_r^2}\right)\exp\left({-z^2/2\sigma_\ell^2}\right),
$$
where $\rho^2=x^2+y^2$, possibly with some polynomial factor of the form $\rho^{|m|}\cos(m\phi)z^n$ in front.
This has an easy solution when the two variances $\sigma_r$ and $\sigma_\ell$ are equal, in which case there is spherical symmetry applies and Gauss's law readily yield $\Phi$ in terms of error functions.
However, I'm having trouble finding results for the general case. Using spheroidal coordinates doesn't seem to help, as the constant-$\mu$ surfaces are confocal spheroids that taper out to spheres as they get larger, instead of maintaining their ellipticity as the charge density does.
This looks standard enough that it ought to have been done before (right?), but it is also messy enough that I'm unsure it made it to any textbook. Has anyone seen the likes of this before?
 A: $\def\a{\alpha}
\def\p{\pi}
\def\s{\sigma}
\def\f{\varphi}
\def\S{\Sigma}
\def\r{\varrho}
\def\vr{\mathbf{r}}
\def\ur{\mathbf{\hat{r}}}
\def\o{\cdot}
\def\RR{\mathbb{R}}$Not a complete solution, but too long for a comment. 
Consider the charge density 
$$\r(x,y,z) = Q \f_1(x)\f_2(y)\f_3(z),$$
where 
$$\f_i(u) = \frac{1}{\sqrt{2\p}\s_i} e^{-u^2/(2\s_i^2)}.$$
Note that $\f_i$ is a normal distribution with zero mean and standard deviation $\s_i$. 
Thus, 
$\int \r(\vr) dV = Q.$
We find 
\begin{align*}
\int \frac{\r(\vr')}{|\vr-\vr'|}dV'
&= \sum_{n=0}^\infty \frac{1}{r^{n+1}} \int_{r'\le r} \r(\vr') P_n(\ur\o\ur'){r'}^n dV' \\
&\quad + \sum_{n=0}^\infty r^n \int_{r'>r} \r(\vr') P_n(\ur\o\ur')\frac{1}{{r'}^{n+1}} dV',
\end{align*}
where $P_n$ is the $n$th Legendre polynomial. 
For $r\gg \max(\s_i)$ this will be well-approximated by 
$$I = \sum_{n=0}^\infty \frac{1}{r^{n+1}} \int_{\RR^3} \r(\vr') P_n(\ur\o\ur'){r'}^n dV'.$$
Since $P_n$ is odd if and only if $n$ is odd we find 
$$I = \sum_{m=0}^\infty \frac{c_{2m}}{r^{4m+1}},$$
where 
$$c_{2m} = \int_{\RR^3} \r(\vr') P_{2m}(\ur\o\ur')(r {r'})^{2m} dV'.$$
(Note, for example, that $P_1(\ur\o\ur'){r r'}=\vr\o\vr'=x x'+y y'+z z'$, which integrates to zero. Thus, the dipole, octupole, ... contributions are zero.)
The first few $c$s are given as 
\begin{align*}
c_0 &= Q \\
c_2 &= \frac{Q}{2} (3\S_{2,2}-\S_{2,0}r^2) \\
c_4 &= \frac{Q}{8} 
    \left(105 \S_{2,2}^2
    -30(\S_{2,2}\S_{2,0}+2\S_{4,2}) r^2
    +3(\S_{2,0}^2+2\S_{4,0}) r^4\right),
\end{align*}
where 
\begin{align*}
\S_{m,n} &= \sum_{i=1}^3 \s_i^m x_i^n.
\end{align*}
This gives the monopole, quadrupole, and hexadecapole contributions to $\Phi$ for $r\gg\max(\s_i)$. 
If $\s_1=\s_2 = \s$ we find
\begin{align*}
c_2 &= \frac{Q}{2}(\s^2-\s_3^2)(x^2+y^2-2z^2) \\
c_4 &= \frac{3Q}{8} (\s^2-\s_3^2)^2(3(x^2+y^2)^2-24(x^2+y^2)z^2+8z^4)
\end{align*}
In principle the $c$s may be calculated generally.
The coefficients of the Legendre polynomials are known and we may expand powers of $\vr\o\vr'$ and ${r'}^{2m}$ as trinomial series. 
The problem is reduced to finding moments of the normal distribution.
We expect corrections due to the assumption that $r\gg\max\s_i$ to be of order 
$\frac{1}{r}\left(1-\mathrm{erf}
\sqrt{\sum_i \frac{x_i^2}{2\s_i^2}}\right)$.
