integral of spherical harmonics over cube The complex solid spherical harmonics can be defined as
$$
U_n^m(\boldsymbol{r}) = r^n P_n^m(\cos{\theta}) e^{im\phi},
$$
where $r,\theta,\phi$ are the usual spherical coordinates of $\boldsymbol{r}=(x,y,z)$. Note that $U_n^m(\boldsymbol{r})$ is a homogeneous polynomial in $x$, $y$, and $z$ of degree $n$, i.e. $U_n^m(\alpha\boldsymbol{r})=\alpha^n U_n^m(\boldsymbol{r})$. I want to work out the integrals
$$
I_n^m := \int_{-1}^1 dx \int_{-1}^1 dy \int_{-1}^1 dz \;U_n^m(\boldsymbol{r})
$$
but cannot get it (nor did I manage to have maple get it for me). Note that since the $U_n^m(\boldsymbol{r})$ are homogenous, we can reduce this to a surface integral:
$$
I_n^m = \frac{1}{n+3} \int \partial \Omega \;U_n^m(\boldsymbol{r})
$$
where the intregral is over the surface of the unit cube. Obviously, $I_n^m=0$ for any odd $n$ or $m$, but I need the full story. (I want to get the multipole moments for a uniform density within a cube).
 A: I've solved the problem numerically using a C++ code and found that not only to odd n's vanish, all m's not divisible by 4 vanish - and to my surprise all quadrupole moments vanish. So the first non-zero term past n=0 is n=4 with m=-4,0,+4, then n=6 with m=-4,0,+4, then n=8 with m=-8,-4,0,+4,+8, etc. 
Just curious - do you need this to do multipole calculations in a grid code? 
EDIT - also, all imaginary parts vanish for all.
A: Here are details of the strategy outlined in the comments for recursively calculating the coefficients of $U_n^m(\mathbf{r})$. This doesn't fully answer the question, but perhaps it will be useful nonetheless.
To fix notation, let $(r, \phi, \theta)$ denote spherical coordinates with the physicists' convention:
$$
x = r\cos\phi \sin\theta,\quad
y = r\sin\phi \sin\theta,\quad
z = r\cos\theta.
$$
For non-negative integers $n$ and $m$, write $U_n^m(\mathbf{r}) = r^n P_n^m(\cos\theta) e^{im\phi}$ as
$$
\sum_{i, j, k = 0}^\infty (a_n^m)_{i,j,k} x^i y^j z^k,
$$
with the understanding that the coefficient $(a_n^m)_{i,j,k}$ is zero if any index is negative, or is larger than $n$. (Particularly, each sum is finite.)
The first recursion relation for the Legendre polynomials becomes, upon replacing "$x$" by $\cos\theta$, setting $\ell = n + 1$,
$$
(n - m + 2) P_{n+2}^m(\cos\theta)
  = (2n + 3) \cos\theta P_{n+1}^m(\cos\theta)
  - (n + m + 1) P_n^m(\cos\theta).
$$
Multiplying through by $r^{n+2} e^{\sqrt{-1}m\phi}$,
$$
(n - m + 2) U_{n+2}^m(\mathbf{r})
  = (2n + 3) z U_{n+1}^m(\mathbf{r}) - (n + m + 1) r^2 U_n^m(\mathbf{r}).
$$
Writing each polynomial $U_n^m$ in terms of its coefficients, shifting indices, and equating coefficients of $x^i y^j z^k$ gives the coefficient relation
\begin{align*}
(a_{n+2}^m)_{i,j,k}
  &= \frac{2n + 3}{n - m + 2} (a_{n+1}^m)_{i,j,k-1} \\
  &\quad- \frac{n + m + 1}{n - m + 2}
  \bigl((a_n^m)_{i-2,j,k} + (a_n^m)_{i,j-2,k} + (a_n^m)_{i,j,k-2}\bigr).
\end{align*}
Since $U_0^0(\mathbf{r}) = 1$, we have $(a_0^0)_{0,0,0} = 1$ and $(a_0^0)_{i,j,k} = 0$ for all other $i$, $j$, $k$.
Since $U_1^0(\mathbf{r}) = z$, we have $(a_1^0)_{0,0,1} = 1$ and $(a_1^0)_{i,j,k} = 0$ for all other $i$, $j$, $k$.
The preceding recurrence determines the coefficients $(a_n^0)_{i,j,k}$.
Similarly, $U_0^1(\mathbf{r}) = 0$, so $(a_0^1)_{i,j,k} = 0$ for all $i$, $j$, and $k$; and $U_1^1(\mathbf{r}) = -(x + \sqrt{-1} y)$, so $(a_1^1)_{1,0,0} = -1$, $(a_1^1)_{0,1,0} = -\sqrt{-1}$, while $(a_1^1)_{i,j,k} = 0$ for all other $i$, $j$, $k$.
The preceding recurrence determines the coefficients $(a_n^1)_{i,j,k}$.
The second recursion relation for the Legendre polynomials becomes, upon replacing "$x$" by $\cos\theta$ and $m$ by $m + 1$, and rearranging,
$$
\sin\theta P_n^{m+2}(\cos\theta)
  = -2(m + 1) \cos\theta P_n^{m+1}(\cos\theta)
  - (n + m + 1)(n - m) \sin\theta P_n^m(\cos\theta).
$$
Multiplying through by $r^{n+1} e^{\sqrt{-1}(m+1)\phi}$ and using $r\sin\theta e^{\pm \sqrt{-1}\phi} = x \pm \sqrt{-1} y$ gives
$$
(x - \sqrt{-1} y) U_n^{m+2}(\mathbf{r})
  = -2(m + 1) z U_n^{m+1}(\mathbf{r})
  - (x + \sqrt{-1} y)(n + m + 1)(n - m) U_n^m(\mathbf{r}).
$$
Writing each polynomial $U_n^m$ in terms of its coefficients, shifting indices, and equating coefficients of $x^i y^j z^k$ gives the coefficient relation
\begin{align*}
&(a_n^{m+2})_{i-1,j,k} - \sqrt{-1} (a_n^{m+2})_{i,j-1,k} \\
  &\qquad= -2(m + 1) (a_n^{m+1})_{i,j,k-1} \\
  &\qquad\quad- (n + m + 1)(n - m) \bigl((a_n^m)_{i-1,j,k} + \sqrt{-1}(a_n^m)_{i,j-1,k}\bigr).
\end{align*}
Since the $(a_n^0)_{i,j,k}$ and $(a_n^1)_{i,j,k}$ are known for all $n \geq 0$, the coefficients $(a_n^m)_{i,j,k}$ may be found recursively, and the integrals $I_n^m$ computed as noted in the comments.
