Why is quadruple moment zero for spherically symmetric charge distribution about centre? How can we show that for a spherically symmetric charge distribution, the dipole, quadrupole and all higher moments about the centre of the distribution are identically zero.

As we already know that for a spherically symmetric charge distribution potential is first term in equation 1. From that how can we conclude that higher order terms are zero, Why can't they just cancel each other?
 A: What you have done is correct. In my answer I will show why the terms with $l>0$ vanish. As you have shown,
$$\Phi = \frac{1}{4\pi\epsilon_0r}\sum^{\infty}_{l=0}\int \rho(r') P_{l}(\cos{\alpha})\bigg(\frac{r'}{r}\bigg)^{l} dV'$$
Since the distribution is spherically symmetric I have considered the origin to be the centre of the sphere. The cosine of the angle made by $r'$ and $r$ will be,
$$\cos(\alpha) = \cos(\phi' - \phi) \sin(\theta)\sin(\theta') + \cos(\theta)\cos(\theta')$$
Using the addition theorem of spherical harmonics,
$$P_{l}(\cos(\alpha)) = \frac{4\pi}{2l+1} \sum^{l}_{m=-l} Y_{ml}(\theta', \phi')Y_{ml}^{*}(\theta, \phi)$$
Substituting this in the first equation,
$$\Phi = \frac{1}{4\pi\epsilon_0r}\sum^{\infty}_{l=0}\int \rho(r') \bigg(\frac{4\pi}{2l+1}\bigg)\times \sum^{l}_{m=-l} Y_{ml}(\theta', \phi')Y_{ml}^{*}(\theta, \phi)\bigg(\frac{r'}{r}\bigg)^{l} dV'$$
If you look only at the angular part of the volume integration,
$$\int^{\pi}_{0}\int^{2\pi}_{0} Y_{ml}(\theta', \phi') d\Omega'$$
It comes out to be $\delta_{0l}\delta_{0m}$. Since there is sum over both $m$ and $l$ the only non vanishing term will be $l=0,m=0$. Hence
$$\Phi = \frac{1}{4\pi\epsilon_0r}\int^{R}_{0}\rho(r') (4\pi r'^2dr')$$
Which is the same as,
$$\Phi = \frac{Q}{4\pi\epsilon_0r}$$

Proof of Orthonormality
To solve the integral above you will require the orthonormality of spherical harmonics. Explicitly writing the spherical harmonics,
$$Y_{ml} (\theta, \phi) = A_{ml} P_{ml}(\cos(\theta))e^{im\phi}$$
Where
$$A_{ml}= \sqrt{\frac{(2l+1)(l-m)!}{(l+m)!}}$$
Since it's obvious that,
$$\int^{2\pi}_{0} e^{i\phi(m-m')} d\phi = \delta_{m m'}$$
We can write,
$$\int^{\pi}_{0}\int^{2\pi}_{0} Y_{ml}(\theta, \phi) Y_{m'l'}(\theta, \phi) d\phi d\theta = A_{ml}A_{m'l'}\delta_{mm'} \int^{1}_{-1} P_{ml}(x)P_{m'l'}(x)dx$$
$P_{m'l'}$ are known to be orthogonal for same $m$ (see: this proof). Hence,
$$\int^{\pi}_{0}\int^{2\pi}_{0} Y_{ml}(\theta, \phi) Y_{m'l'}(\theta, \phi) d\phi d\theta = \delta_{mm'}\delta_{ll'}$$
For $m'=0, l'=0$,
$$\int^{\pi}_{0}\int^{2\pi}_{0} Y_{ml}(\theta', \phi') d\Omega' = \delta_{m0}\delta_{l0}$$
A: Dipole, quadrupole, octupole etc electric moments arise from non-radial components of an electric field. By symmetry, any spherically symmetric charge distribution can have only a radial electric field. All non-radial components are zero. Hence all multi-pole electric moments are zero.
