What's the Quickest Way to Integrate this Formidable Vector Integral?

[fluxions]

Suppose $\vec{M}, \vec{P}$ are arbitrary, constant vectors, and $\hat{r}$ is the (unit) position vector in spherical polar coordinates.

I need to integrate the vector function $\frac{1}{r^6}[\hat{r} \times ((\vec{P}\cdot \hat{r})\vec{M} - ((\vec{M} \cdot \hat{r})\vec{P})]$ over the entire exterior of the sphere of radius R centered at the origin of coordinates. In other words, I need to compute:

$$\int_{\phi = 0}^{2\pi} \int_{\theta = 0}^{\pi} \int_{r = R}^{\infty} \frac{1}{r^4}[\hat{r} \times ((\vec{P}\cdot \hat{r})\vec{M} - ((\vec{M} \cdot \hat{r})\vec{P})] sin\theta dr d\theta d \phi$$

I'm looking for a cute and clever way to do this, instead of the straightforward and tedious method. Any ideas or hints?

 

[fzero]

The slickest way to do this is to write

$$[r \times ((\vec{P}\cdot r)\vec{M} - ((\vec{M} \cdot r)\vec{P})]_i = \sum_{jklmn} \epsilon_{ijk} r_j r_l P_m M_n ( \delta_{lm}\delta_{kn} - \delta_{ln}\delta_{km} )$$

and write the $$r_i$$ in terms of spherical harmonics:

$$\frac{r_1}{r} = - \sqrt{\frac{2\pi}{3}} \left( Y_1^1 + Y_1^{-1} \right),$$

$$\frac{r_2}{r} = i \sqrt{\frac{2\pi}{3}} \left( Y_1^1 - Y_1^{-1} \right),$$

$$\frac{r_3}{r} = 2 \sqrt{\frac{\pi}{3}} Y_1^0 .$$

The angular integration is done using the orthogonality relations and gives you a matrix $$Q_{jl}$$ that you then have to sum over. Whether it saves that much work over brute force is to be determined.