How to find that triple integral? How to find the triple integral of $$ \frac{(z-z_0)z}{\sqrt{(x-x_0)^2+(y-y_0)^2+(z-z_0)^2}}$$ over the sphere $ \{(x,y,z):x^2+y^2+z^2 \le 1 \}  $ under the assumption $x_0^2+y_0^2+z_0^2 \le 1?$ Its physical interpretation suggests
 the integral can be expressed through elementary functions of the parameters.

Let $\vec{p} = (x,y,z)$ and $\vec{p}_0 = (x_0,y_0,z_0)$. Let $(r,\theta,\phi)$ and $(r_0,\theta_0,\phi_0)$ be their spherical polar coordinates. More precisely,
$$\begin{array}{lll}
\vec{p}   &= (x,y,z)       &= r(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)\\
\vec{p}_0 &= (x_0,y_0,z_0) &= r_0(\sin\theta_0\cos\phi_0,\sin\theta_0\sin\phi_0,\cos\theta_0)
\end{array}$$
Let $\gamma$ be the angle between $\vec{p}$ and $\vec{p}_0$. i.e
$$\cos\gamma = \cos\theta\cos\theta_0 + \sin\theta\sin\theta_0\cos(\phi-\phi_0)$$
The integral $\mathcal{I}$ at hand can be rewritten as
$$\mathcal{I} = \int_{r\le 1} \frac{r^2\cos\theta^2 -rr_0\cos\theta\cos\theta_0}{\sqrt{r^2 + r_0^2 - 2rr_0\cos\gamma}} r^2 dr d\Omega$$
where $d\Omega = \sin\theta d\theta d\phi$ is the differential element for solid angle.
In terms of Legendre polynomials $P_\ell(x)$, we have
$$r^2\cos\theta^2 -rr_0\cos\theta\cos\theta_0 = 
r^2\frac{2P_2(\cos\theta) + P_0(\cos\theta)}{3} - rr_0 P_1(\cos\theta)P_1(\cos\theta_0)$$
and
$$\frac{1}{\sqrt{r^2 + r_0^2 - 2rr_0\cos\gamma}} = \sum_{\ell=0}^\infty \frac{r_<^\ell}{r_>^{\ell+1}} P_\ell(\cos\gamma)
\quad\text{ where }\quad\begin{cases}r_< &= \min(r,r_0)\\ r_> &= \max(r,r_0)\end{cases}
$$
Using following identity for Legendre polynomials
$$\int P_\ell(\cos\gamma)P_{\ell'}(\cos\theta) d\Omega = \frac{4\pi}{2\ell+1}\delta_{\ell\ell'} P_{\ell'}(\cos\theta_0)$$
We find
$$\mathcal{I} = 
\frac{8\pi}{15}A_{2,4}(r_0) P_2(\cos\theta) 
+ \frac{4\pi}{3}A_{0,4}(r_0) P_0(\cos\theta)
- \frac{4\pi}{3} A_{1,3}(r_0) r_0 P_1(\cos\theta_0)^2$$
where
$$\begin{align}
A_{\ell,s}(r_0) 
&= \int_0^1 \frac{r_<^\ell}{r_>^{\ell+1}} r^s dr
= \int_0^{r_0} \frac{r^{\ell+s}}{r_0^{\ell+1}} dr + \int_{r_0}^1 r_0^\ell r^{s-\ell-1} dr\\
&= \frac{r_0^s}{\ell+s+1} + r_0^\ell \left(\frac{1 - r_0^{s-\ell}}{s-\ell}\right)
= \frac{r_0^\ell}{s-\ell} + \frac{r_0^s(2\ell+1)}{\ell(\ell+1) - s(s+1)} 
\end{align}
$$
As a result,
$$\begin{align}
\mathcal{I} &=
\frac{8\pi}{15}\left(\frac{r_0^2}{2} - \frac{5r_0^4}{14}\right) \frac{3 \cos\theta_0^2 - 1}{2}
+ \frac{4\pi}{3}\left(\frac14 - \frac{r_0^4}{20}\right)
- \frac{4\pi}{3} \left(\frac{r_0^2}{2} - \frac{3r_0^4}{10}\right) \cos\theta_0^2\\
&= \frac{4\pi}{210}\left(7-5r_0^2\right)\left(3z_0^2 - r_0^2\right)
+ \frac{4\pi}{60}\left(5 - r_0^4\right)
- \frac{4\pi}{30}\left(5-3r_0^2\right)z_0^2
\end{align}
$$
