Traceless multipole moments vs non-traceless moments There are two different possibilities to define the electric quadrupole tensor:
On the one hand, one can define
\begin{align}Q_{kl} = \int \rho(\mathbf r') \cdot r'_k \, r'_l  d^3r',\end{align}
while on the other hand most of the textbooks (for example Jackson) define
\begin{align}
Q_{kl} = \int \rho(\mathbf r') \cdot (3r'_k \, r'_l - (r')^2 \, \delta_{kl}) \cdot d^3r'.\end{align}
For all higher moments, these two definitions (traceless vs. nontraceless) are also possible. Is there any physical reason why one of the definitions should be preferred?
 A: The reason the traceless moments are chosen is that they form an irreducible representation of the rotation group. 
In contrast, the non-traceless moments contain a nontrivial subspace which is invariant under rotations, which is obviously the span of the diagonal tensor
$$
Q_{kl}=\delta_{kl}.
$$
(That is, the invariant subspace is the space of all diagonal tensors.)
The existence of this invariant subspace is, for various reasons, undesirable. In essence, doing a multipole expansion is about estimating the potential of the system on a sphere of large radius $r$. Because of this, you want the functions in your asymptotic expansion to play well with rotations; more specifically, you want them to embody representations of the rotation group.
As it turn out, all the potentials in a given irreducible representation must have the same dependence on $r$. This encodes the fact that if the angular dependence is of similar character, you can rotate any function within the representation subspace into any other, in the same way you can rotate a $z$-pointing dipole into an $x$-pointing one. (In reality, this is a little more complex, and if you rotate an $xz$ quadrupole you won't get a single quadrupole term, but you will be able to express the result as a sum of pure quadrupoles. This is also true for dipoles if the rotation is arbitrary.)
Thus, irreducible representations are forced by the spherical symmetry to have the same radial dependence throughout, and one usually wants to make this apparent. The presence of invariant subspaces would break this property, which is why we use notation to explicitly take those out.

More practically, though, the reason the traceless moments are the ones to use is that they actually do describe the far-field potential of the system. To see this, consider a charge density $\rho$ which is contained inside a bounded volume $\Omega$ around the origin. Its electrostatic potential for $\mathbf r$ outside this volume can then be approximated as
\begin{align}
\Phi(\mathbf r) 
& = \int_\Omega\mathrm d\mathbf r'\, 
\frac{\rho(\mathbf r')}{|\mathbf r-\mathbf r'|}
 = \frac1r \int_\Omega\mathrm d\mathbf r'\, 
\frac{\rho(\mathbf r')}{\left(1-2\mathbf r\cdot \mathbf r'/r^2+r'^2/r^2\right)^{1/2}}
\\&=
 \frac1r \int_\Omega\mathrm d\mathbf r'\,\rho(\mathbf r')
\sum_{n=0}^\infty \begin{pmatrix}-1/2\\n\end{pmatrix}
\left(-2\mathbf r\cdot \mathbf r'/r^2+r'^2/r^2\right)^{n}
\\&\approx
 \frac1r \int_\Omega\mathrm d\mathbf r'\,\rho(\mathbf r')
\left[
1
-\frac12 \left(-2\mathbf r\cdot \mathbf r'/r^2+r'^2/r^2\right)
+\frac38 \left(-2\mathbf r\cdot \mathbf r'/r^2+r'^2/r^2\right)^{2}
+\cdots
\right].
\end{align}
If you want to approximate this potential up to quadrupole order, then you can just truncate every term which scales worse than $r^{-3}$. This means, for starters, that you can neglect every term in the last square except the square of $-2\mathbf r\cdot \mathbf r'/r^2$. However, you will have another term that contributes to second order in $r'$, which comes from the contribution at $n=1$, and which has a distinct, isotropic character.
If you do this truncation procedure, you get
\begin{align}
\Phi(\mathbf r)
 &\approx
 \frac1r \int_\Omega\mathrm d\mathbf r'\,\rho(\mathbf r')
\left[
1
+\frac{\mathbf r\cdot \mathbf r'}{r^2}
+\frac32\frac{(\mathbf r\cdot \mathbf r')^2}{r^4}-\frac12\frac{r'^2}{r^2}
\right]
\\ &=
 \frac1r \int_\Omega\mathrm d\mathbf r'\,\rho(\mathbf r')
\left[
1
+\frac{x_k x_k'}{r^2}
+\frac12\frac{3x_k x_l x_k'x_l'-r^2 r'^2}{r^4}
\right],
\end{align}
where I'm using Einstein summation notation to keep things simple.
This is almost ready to perform the essential step of the multipole approximation, which is to separate out the different terms corresponding to the different orders,
\begin{align}
\Phi(\mathbf r)
 &\approx
 \frac1r \int_\Omega\rho(\mathbf r')\mathrm d\mathbf r'
+
\frac1r \int_\Omega\frac{x_k x_k'}{r^2} \rho(\mathbf r')\mathrm d\mathbf r'
+
\frac1r \int_\Omega\frac12\frac{3x_k x_l x_k'x_l'-r^2 r'^2}{r^4}\rho(\mathbf r')\mathrm d\mathbf r'
,
\end{align}
and then to factor these out into potentials (which depend on $\mathbf r$) and moments (which depend on $\mathbf r'$). To get it into this form, you need to rewrite the term $r^2 r'^2$ as $x_k x_l \delta_{kl} r'^2$, so that the $x_kx_l$ term will factorize, and when you do that you get
\begin{align}
\Phi(\mathbf r)
 &\approx
 \frac1r \int_\Omega\rho(\mathbf r')\mathrm d\mathbf r'
+
\frac{x_k}{r^3} \int_\Omega x_k' \rho(\mathbf r')\mathrm d\mathbf r'
+
\frac{x_k x_l}{r^5} \int_\Omega\frac{3 x_k'x_l'-\delta_{kl} r'^2}{2}\rho(\mathbf r')\mathrm d\mathbf r'
.
\end{align}
As you can see, the traceless moments emerge naturally, from a careful and consistent treatment of all the terms in the expansion. In particular, the trace-cancelling term $-\delta_{kl} r'^2$ comes from a term one order down - from the squared $r'^2/r^2$ term that was excluded from the dipole term. These are hard to keep track of, but they are essential for a correct description, because they contribute at the same order as the quadrupole ones you want to describe.
You could, of course, define the quadrupole moments as the non-traceless ones, but then you would have to put in extra work into the form of the potential, each term of which which would then depend on more than one moment. In practice no one does this, because it takes you even further away from the spherical-harmonics description of the moments which is their natural representation.
