I have a very simple problem: There is a charge $-q$ at $(0, 0, d)$ and $(0, 0, -d)$ as well a charge $2q$ at $(0, 0, 0)$. I have to calculate the quadrupole moment using spherical coordinates. I use $x=\cos(\phi)$.

The formula I build up is this: $$q_{2,m} = \frac{(2-m)!}{(2+m)!} \frac{(-1)^{m+1}}8 \int_{-1}^1 \mathrm dx \, \left(1-x^2\right)^{m/2} \frac{\mathrm d^{2+m}}{\mathrm dx^{2+m}} \left(x^2-1\right)^2$$ $$\times\int_0^{2\pi} \mathrm d \phi \, \mathrm e^{\mathrm im \phi} \int_0^\infty \mathrm dr \, r^2 \rho(r, x, \phi).$$

Is this even correct?

Then I tried to figure out what the $\rho$ has to be, this is what I currently believe to be correct: $$\rho(r, x, \phi) = q \delta(\phi) \left[2 \delta(r) \delta(x-1) - \delta(r-d) \delta(x-1) - \delta(r-d) \delta(x+1)\right].$$

So now I would just need to calculate the integrals, right? The problem for positive $m$ would be that $\int \mathrm dx \, \delta(x-1) \left(1-x^2\right)^{m/2}$ would be zero. And $\int \mathrm dr \, r^2 \delta(r)$ would also be zero. So the $q_{2,1}$ and $q_{2,2}$ are zero.

I previously used the following charge density (although there was a negative radius): $$\rho(r, x, \phi) = q \delta(x) \delta(\phi) \left[2\delta(r) - \delta(d-r) - \delta(d+r)\right].$$

And I got: $$q_{2,-2} = 3! q d^2 ,\quad q_{2,-1} = 0 ,\quad q_{2,0} = -q d^2 ,\quad q_{2,1} = 0 ,\quad q_{2,2} = \frac 14 qd^2.$$

In one of the books, it says that only the $Q_{3,3}$ element of the tensor would be non-zero.

So how can I obtain the correct answer here?

-
By symmetry, only the $q_{2,0}$ will contribute to the quadrupole moments, and it will be the leading term to the multipole expansion. The $m=\pm 2$ components require that the charge change sign after a $90^\circ$ rotation, like you would get with unit charges at $(\pm1,0,0)$, positive, and $(0,\pm1,0)$, negative. – Emilio Pisanty Nov 7 '12 at 20:21
If I have the $\delta(x-1)$ and $\delta(x+1)$, I will still have 0 for $m=0$. If I take $\delta(x)$, i should have a nonzero answer. Which $x$ do I have to choose in the origin? – Martin Ueding Nov 7 '12 at 20:29
There's some additional obvious problems with your answer. For one, the $q_{2,\pm2}$ components must always be complex conjugates (i.e. equal if they're real) for the total field to be real-valued. This comes about because $Y_{l,-m}=Y_{lm}^\ast$ as $Y_{lm}\propto e^{im\phi}$. – Emilio Pisanty Nov 7 '12 at 20:52
Ah! yes, you're using $x=\cos\theta$. Those multiple delta functions looked funny but they are in fact OK. That use for $x$ in this context is quite misleading; more traditional notations are e.g. $u$ or $v$. – Emilio Pisanty Nov 7 '12 at 20:54

I would say the quickest route would be to use the multipole expansion for the potentials of the different charges and then add those up. Thus you'd use $$\frac{1}{|\mathbf{r}-\mathbf{r}'|}=\sum_{l=0}^\infty\sum_{m=-l}^l\frac{ 4\pi}{2l+1} \frac{r'^l}{r^{l+1}}Y_{lm}^\ast(\theta',\phi')Y_{lm}(\theta,\phi),$$ valid for $r>r'$ (see Jackson, 2nd ed., p.111, eq 3.70). To implement this, take $\mathbf{r}'=(0,0,\pm d)$, so that $\theta'=0,\pi$ and $\phi'$ is irrelevant. Then $Y_{lm}^\ast(\theta',\phi')$ is zero unless $m=0$, so you have $$Y_{lm}^\ast(\theta',\phi')=\delta_{m,0}N_{lm}P_l(\cos\theta)=\delta_{m,0}\sqrt{\frac{2l+1}{4\pi}}P_l(\pm1)=\delta_{m,0}\sqrt{\frac{2l+1}{4\pi}}(\pm1)^l.$$ Thus, in this case, $$\frac{1}{|\mathbf{r}-\mathbf{r}'|}=\sum_{l=0}^\infty (\pm d)^l \sqrt{\frac{ 4\pi}{2l+1}} \frac{Y_{l0}(\theta,\phi)}{r^{l+1}}=\sum_{l=0}^\infty (\pm d)^l \frac{P_{l}(\cos\theta)}{r^{l+1}}.$$
If you put it all together, then, the $l=0$ components of the charges at $(0,0,\pm d)$ will cancel with the one at the origin (all the other multipole terms of which vanish, of course), and the $l=1$ terms will kill each other due to that sign. The full potential will therefore be $$\Phi=q\sum_{l=2}^\infty \left[1+(-1)^l \right] d^l \frac{P_{l}(\cos\theta)}{r^{l+1}}=\sum_{k=0}^\infty 2q d^{2+2k} \frac{P_{2+2k}(\cos\theta)}{r^{3+2k}}.$$
Now, as to how to extract the multipole moments, it depends on what normalization you use. I would use a normalization as $$\Phi=\sum_{l=0}^\infty\sum_{m=-l}^l q_{lm}\frac{Y_{lm}(\theta,\phi)}{r^{l+1}},$$ but it depends on what the multipolar fields are defined like in your book. Since the nonzero-$m$ terms have already vanished, we must have $q_{2,\pm1}=q_{2,\pm2}=0$ as the problem's symmetry demands (as I indicated already in the comments). The $m=0$ component is then simply $$q_{2,0}=2q d^2,$$ but, again, this depends on the normalization you're using.