First understand the method of image charges. The idea behind the method is to bypass actually solving the differential equation with boundary conditions, and instead "cheat" by guessing the correct solution. To this end, we find a configuration of imaginary charges that together with the real ones will make the potential on all surfaces be what is given.
In your case, you have two surface, each with constant potential zero. For a plane, if you have a charge $q$ at $(x,y,z)$ and another one at $(x,-y,z)$, clearly the potential at $y=0$ will be zero. 586
For the sphere, if a charge is at radius $R$ from the sphere with radius $a$, an image charge with charge $-qa/R$ should be placed at radius $a^2/R$, but the same idea remains.
Now, first add the image charge for the sphere. Now the potential on the sphere is zero, but on the plane, we still have a gradient. However, if you mirror both charges off the plane, you should see that the potential on both the sphere and the plane is zero.
Writing this all down, we have: $$\begin{eqnarray} \phi({\bf x}) &=& \frac{q}{\sqrt{(x-R\cos\alpha)^2+(y-R\sin\alpha)^2+z^2}} \ &-& \frac{q}{\sqrt{(x-R\cos\alpha)^2+(y+R\sin\alpha)^2+z^2}} \ &-&\frac{qa}{R\sqrt{(x-(a^2/R)\cos\alpha)^2+(y-(a^2/R)\sin\alpha)^2+z^2}} \ &+&\frac{qa}{R\sqrt{(x-(a^2/R)\cos\alpha)^2+(y+(a^2/R)\sin\alpha)^2+z^2}} \end{eqnarray}$$ Clearly, at $y=0$ we have $\phi=0$. Now, what about on the hemisphere? points on the hemisphere satisfy $z^2 = a^2 - x^2-y^2$, so that substituting leads us to the potential on the hemisphere:
$$\begin{eqnarray} \phi({\bf x}_{sphere}) &=& \frac{q}{\sqrt{-2xR\cos\alpha -2yR\sin\alpha + R^2 +a^2}} \ &-& \frac{q}{\sqrt{-2xR\cos\alpha +2yR\sin\alpha + R^2 +a^2}} \ &-&\frac{qa}{R\sqrt{-2x(a^2/R)\cos\alpha -2y(a^2/R)\sin\alpha + (a^2/R)^2 + a^2}} \ &+&\frac{qa}{R\sqrt{-2x(a^2/R)\cos\alpha +2y(a^2/R)\sin\alpha + (a^2/R)^2 + a^2}} \ &=&0 \end{eqnarray}$$ By putting the $a/R$ into the square root.
To write down the multipole expansion, you just need to write down the <PERSON> expansion of the potential around $1/r$, with $r = \sqrt{x^2+y^2+z^2}$. This gives that any expression of the form: $$\lim_{r\to\infty}\frac{1}{\sqrt{r^2+b}}= \frac{1}{r} - \frac{b}{2r^3} + \frac{3b^2}{8r^5} + \cdots$$ You can then use this expansion on the components of the potential to get your result. 521