Electric octupole moment in cartesian coordinates I'm trying to calculate the symmetric traceless tensor for the octupole moment in cartesian coordinates... I have to deal with the electrostatic potential of the form:
$\Phi^{(4)}(\textbf{r})=\int d^{3}r^{\prime}\rho(\textbf{r}^\prime)\bigg[\frac{1}{3!}\partial^{\prime}_{k}\partial^{\prime}_{j}\partial^{\prime}_{i}\bigg(\frac{1}{|\textbf{r}-\textbf{r}^\prime|}\bigg)\bigg|_{\textbf{r}^\prime=0}r_{i}^{\prime}r_{j}^{\prime}r_{k}^{\prime}\bigg]$
where the term in square brackets corresponds to the third term in the Taylor series of $1/|\textbf{r}-\textbf{r}^\prime|$. Doing the partial derivatives, the electrostatic potential may be rewritten as:
$\Phi^{(4)}(\textbf{r})=\frac{1}{3!r^{7}}\bigg[15r_{i}r_{j}r_{k}-3(r_{i}\delta_{jk}+r_{j}\delta_{ik}+r_{k}\delta_{ij})r^{2}\bigg]C_{ijk}$
where $C_{ijk}=\int d^{3}r^{\prime}\rho(\textbf{r}^\prime)r_{i}^{\prime}r_{j}^{\prime}r_{k}^{\prime}$ is the non-traceless octupole moment tensor.

The question is: How can I rearrenge the expression for the potential so that I get something of the form:
$\Phi^{(4)}(\textbf{r})=\frac{r_{i}r_{j}r_{k}}{3!r^{7}}Q^{(4)}_{ijk}$
where $Q^{(4)}_{ijk}$ is the symmetric traceless octupole moment tensor?
 A: I'm surprised this hasn't been answered for more than a year, and this is the third link that comes up on google for "octupole moment". This derivation goes as follows:
For brevity I'll use this parentheses notation for a sum with cyclic permutations of indices:
$$
A_{(i}B_{jk)} \equiv A_iB_{jk}+A_jB_{ki}+A_kB_{ij},
$$
where $A_i$ and $B_{jk}$ are some tensors of order 1 and 2 respectively. 
Let $A_{ijk}$ be a symmetric tensor of third order, and $B_{ijk}$ a symmetric tensor with an additional property that $B_{iij}=0$ (Einstein summation is implied throughout), in other words, its contraction with respect to any pair of indices is zero. We want to see how we can transform the expression 
$$A_{ijk}B_{ijk}.$$
First of all, notice that since $B_{ijk}$ has the contraction property, we can add to $A_{ijk}$ in the expression any tensor of the form $C_i\delta_{jk}$ (and therefore, of the form $C_{(i}\delta_{jk)}$), where $\delta_{jk}$ is the Kronecker delta. This is because 
$$C_i\delta_{jk}B_{ijk}=C_iB_{ijj}=0,$$ 
as $B_{ijj}=0$ by the property. Now we can choose $C_i=-\frac15A_{mmi}$ and obtain:
$$
A_{ijk}B_{ijk}=[A_{ijk}-\frac15A_{mm(i}\delta_{jk)}]B_{ijk}.
$$
The reason for that choice is that now the expression in brackets has the contraction property as you can check by contracting it with respect to any 2 indices.
This lets us add any tensor of the form $C_i\delta_{jk}$ to $B_{ijk}$ to obtain the final expression which we will use:
$$
A_{ijk}B_{ijk}=[A_{ijk}-\frac15A_{mm(i}\delta_{jk)}][B_{ijk}+C_{(i}\delta_{jk)}].
$$
Now going back to your problem, we have $A_{ijk}=r'_ir'_jr'_k$ and $B_{ijk}=15r_ir_jr_k-3r_{(i}\delta_{jk)}r^2$, and obviously they satisfy the required properties, so using the formula we've just obtained:
$$
A_{ijk}B_{ijk}=[r'_ir'_jr'_k][15r_ir_jr_k-3r_{(i}\delta_{jk)}r^2]=[r'_ir'_jr'_k-\frac15r_{(i}'\delta_{jk)}r'^2][15r_ir_jr_k-3r_{(i}\delta_{jk)}r^2+C_{(i}\delta_{jk)}]=[r'_ir'_jr'_k-\frac15r_{(i}'\delta_{jk)}r'^2][15r_ir_jr_k],
$$
where we've chosen $C_i=3r_ir^2$.
Plugging it into your expression for the potential, we now obtain:
$$
\Phi^{(4)}(\textbf{r})=\int d^{3}r^{\prime}\rho(\textbf{r}^\prime)\bigg[\frac{1}{3!}\partial^{\prime}_{k}\partial^{\prime}_{j}\partial^{\prime}_{i}\bigg(\frac{1}{|\textbf{r}-\textbf{r}^\prime|}\bigg)\bigg|_{\textbf{r}^\prime=0}r_{i}^{\prime}r_{j}^{\prime}r_{k}^{\prime}\bigg]=
\int d^{3}r^{\prime}\rho(\textbf{r}^\prime)\bigg[\frac{1}{3!r^7}[r_{i}^{\prime}r_{j}^{\prime}r_{k}^{\prime}][15r_ir_jr_k-3r_{(i}\delta_{jk)}r^2]\bigg]=\int d^{3}r^{\prime}\rho(\textbf{r}^\prime)\bigg[\frac{1}{3!r^7}[r'_ir'_jr'_k-\frac15r_{(i}'\delta_{jk)}r'^2][15r_ir_jr_k]\bigg]=\int d^{3}r^{\prime}\rho(\textbf{r}^\prime)\bigg[\frac{1}{3!r^7}[15r'_ir'_jr'_k-3r_{(i}'\delta_{jk)}r'^2][r_ir_jr_k]\bigg]=\frac{r_ir_jr_k}{3!r^7}Q^{(4)}_{ijk},
$$
where 
$$
Q^{(4)}_{ijk}=\int d^{3}r^{\prime}\rho(\textbf{r}^\prime)[15r'_ir'_jr'_k-3r_{(i}'\delta_{jk)}r'^2]
$$
is the symmetric traceless octupole moment tensor you were looking for.
A: 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$.
