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