Electric octupole moment in cartesian coordinates

I ran into the same question, googled a bit and ended up getting myself a much simpler derivation, so I leave an answer for this 4-year-old question just for the record.

(I'm not sure if any textbook refers to this method - it is as clear as it seems and at least I brought up the proof by myself.)

If you notice that this question is essentially about the formula $$ (15r_{i}r_{j}r_{k}-3r_{[i}delta_{jk]}r^{2})r'_ir'_jr'_k $$ and finding the coefficient of $r_ir_jr_k$s(in terms of $r'_i$s) of it, this is a simple arithmetic problem so you can just expand it and group them appropriately. You don't need some magic tensors out from nowhere.

Okay, the answer is $$ 15r'_ir'_jr'_k-3r'{}_{[i}delta_{jk]}r'^2 $$ as given in the answer of Igor. You can verify it by laborious expansion, or by examining each possibilities of indices (namely, $i,j,k$ are distinct, exactly two are same, and all are same.)

But wait, it's the exact same formula from the original question, only $mathbf{r}$ replaced by $mathbf{r'}$! Observing the same phenomenon in quadrupole terms, we can generalize this to hexadecapoles and any higher $2^l$-pole terms.

Proposition. $1/|mathbf{r}-mathbf{r'}|$ can be expanded in a series $$ frac{1}{|mathbf{r}-mathbf{r'}|} = sum_{l=0}^{infty} frac{F_l (mathbf{r},mathbf{r'})}{r^{2l+1}} $$ where $F_l$ is a homogeneous polynomial of degree $2l$ ($l$ for each $mathbf{r},mathbf{r'}$), and most importantly $$ F_l (mathbf{r},mathbf{r'}) = F_l (mathbf{r'},mathbf{r}). $$

Proof. Try it yourself about the degree $l+l$ part, if you're doubtful. Observe that $$ frac{1}{|mathbf{r}-mathbf{r'}|} = frac{1}{|frac{r}{r'}mathbf{r'}-frac{r'}{r}mathbf{r}|} $$ where $frac{r}{r'}mathbf{r'}$ and $frac{r'}{r}mathbf{r}$ are just vectors $mathbf{r'}$ and $mathbf{r}$, only their magnitudes interchanged.

Applying the series expression to RHS, you get $$ sum_{l=0}^{infty} frac{F_l (mathbf{r},mathbf{r'})}{r^{2l+1}} = sum_{l=0}^{infty} frac{F_l (frac{r}{r'}mathbf{r'},frac{r'}{r}mathbf{r})}{r^{2l+1}} $$ but due to the homogeneousness, $$ F_l left(frac{r}{r'}mathbf{r'},frac{r'}{r}mathbf{r}right) = left(frac{r}{r'}right)^l left(frac{r'}{r}right)^l F_l (mathbf{r'},mathbf{r}) = F_l (mathbf{r'},mathbf{r}) $$ hence we obtain $$ sum_{l=0}^{infty} frac{F_l (mathbf{r},mathbf{r'})}{r^{2l+1}} = sum_{l=0}^{infty} frac{F_l (mathbf{r'},mathbf{r})}{r^{2l+1}} $$ comparing each terms, we conclude that $F_l$ is symmetric. $Box$

Alternatively, you can figure out that $F_l (mathbf{r},mathbf{r'})$ is just the Legendre polynomial $P_l(cos gamma)$ times $r^l r'^l$, where $gamma$ is the angle between $mathbf{r}$ and $mathbf{r'}$.

This also explains why multipole tensors are traceless - because $1/|mathbf{r}-mathbf{r'}|$ is harmonic!

In particular, we know that $$ frac{F_l (mathbf{r},mathbf{r'})}{r^{2l+1}} = partial'_{i_1} cdots partial'_{i_l} left. left(frac{1}{|mathbf{r}-mathbf{r'}|} right)right|_{mathbf{r'}=0} r'_{i_1}cdots r'_{i_l} $$ But if we write $F_l(mathbf{r},mathbf{r'})= M_{i_1 cdots i_l} (mathbf{r'}) r_{i_1}cdots r_{i_l}$, we also have $F_l(mathbf{r},mathbf{r'})= M_{i_1 cdots i_l} (mathbf{r}) r'_{i_1}cdots r'_{i_l}$ so $$ M_{i_1 cdots i_l} (mathbf{r}) = partial'_{i_1} cdots partial'_{i_l} left. left(frac{1}{|mathbf{r}-mathbf{r'}|} right)right|_{mathbf{r'}=0} r^{2l+1} $$ contracting two indices $i_a = i_b = i$ of the tensor sums up to zero because $partial'_i partial'_i (1/|mathbf{r}-mathbf{r'}|)=0$.