- Yep.
- Yep, perfectly fine.
- This will work. The general method is that the expansion coefficients are given by the inner product $\int \psi^\star_{lm} \Psi_0$ where $\psi_{lm}$ are the orthonormal basis functions of your expansion and the integral is over the full space with the natural measure. In your case the basis functions are the spherical harmonics and the measure is $\mathrm{d}\Omega \equiv \sin\theta \mathrm{d}\theta \mathrm{d}\phi$. You can plug in a general spherical harmonic and by doing the integrals get the answer, though in your simple case your method is easier and equivalent. By the way, Mathematica knows about spherical harmonics already (SphericalHarmonicY, but make sure it uses the same conventions as you).
The general method is like finding the Fourier series by doing a bunch of integrals. The method you are using is like reading off the Fourier series coefficients by rearranging the function itself to look like a Fourier series. When you can do it you get the same answer as by doing the integrals, but it's not always easy to do the rearrangement.
RE edits to the question.
Q4. Yep, again. :)
Q5. It's as simple as noting that $L^2 Y_{lm} = l(l+1) Y_{lm}$ and $L_z Y_{lm}= m Y_{lm}$ (I'm dropping factors $\hbar$ which you can put back yourself if need be). Just plug your $\Psi_0$ in $H\Psi_0$ and use the fact that all the operators are linear. It's easy to do term by term and you'll see that you get just the answer you've expected. That's because all of these operators $H,L^2,L_z$ are simultaneously diagonalizable (because they commute) and we're working in the basis of their common eigenfunctions.
There is not complete uniformity in the definition of the vector spherical harmonics, so it is possible that different definitions may actually refer to expressions with different parities. However, if they are defined they way they are in Jackson's Classical Electrodynamics,
$${\bf X}_{l,m}(\theta,\phi)=\frac{1}{\sqrt{l(l+1)}}{\bf L}Y_{l,m}(\theta,\phi),$$
where ${\bf L}$ is the operator $\frac{1}{i}({\bf x}\times{\nabla})$ (which would be the angular momentum divided by $\hbar$ in quantum mechanics), then the parity of ${\bf X}_{l,m}$ is indeed equal to $(-1)^{l}$. This can be discerned fairly straightforwardly from the known parities of ${\bf L}$ and $Y_{l,m}$. Angular momentum, as an axial vector (or an element of a second-rank tensor) is even under parity, which leaves ${\bf X}_{l,m}$ with the same parity as $Y_{l,m}$.
Whenever one talks about the transformation properties of a quantity, the real physical content is associated with how the properties under active transformations. So if you want to find the behaviors of the electrostatic under parity/space inversion, you would consider a situation in which all the charges have their positions relocated—so there is a new charge distribution $\rho'(\mathbf{x})=\rho(-\mathbf{x})$. (Passive transformations, in contrast, are just relabelings of coordinates.) The electric field of $\rho'$ is clearly $\mathbf{E}'(\mathbf{x})=-\mathbf{E}(-\mathbf{x})$, exhibiting the odd parity. On the other hand, a magnetostatic field transforms as $\mathbf{B}'(\mathbf{x})=\mathbf{B}(-\mathbf{x})$, because the source $\mathbf{J}'$ has both its location and direction inverted, $\mathbf{J}'(\mathbf{x})=-\mathbf{J}(-\mathbf{x})$; that is, the current density is polar vector, like the electric field, whereas the magnetic field is an axial vector.
So when looking at the transformation of a vector like $\mathbf{E}$ or $\mathbf{X}_{l,m}$, when you decompose it into components, the components of $\mathbf{E}'$ are opposite those of $\mathbf{E}$: $E'_{x}(\mathbf{x})=-E_{x}(-\mathbf{x})$, $E_{y}'(\mathbf{x})=-E_{y}(-\mathbf{x})$, and $E_{z}'(\mathbf{x})=-E_{z}(-\mathbf{x})$, because $\mathbf{E}'$ the field of the inverted charge distribution $\rho'$ in uninverted coordinates. Since the coordinates are not inverted, the basis vectors are unaffected, and this gives rise to the correct overall transformation
$$\mathbf{E}'(\mathbf{x})=E_{x}'(\mathbf{x})\hat{\mathbf{e}}_{x}+E_{y}'(\mathbf{x})\hat{\mathbf{e}}_{y}+E_{z}'(\mathbf{x})\hat{\mathbf{e}}_{z} \\
=-E_{x}(-\mathbf{x})\hat{\mathbf{e}}_{x}-E_{y}(-\mathbf{x})\hat{\mathbf{e}}_{y}-E_{z}(-\mathbf{x})\hat{\mathbf{e}}_{z}=-\mathbf{E}(-\mathbf{x}).$$
The same thing applies to the unit vectors in the angular directions, $\hat{\mathbf{e}}_{r}$, $\hat{\mathbf{e}}_{\theta}$, and $\hat{\mathbf{e}}_{\phi}$. They are unchanged under a parity transformations, but the three corresponding components of the vector spherical harmonic $\mathbf{X}_{l,m}(\theta,\phi)$ do change their signs.
Best Answer
Using: $$ \mathbf{\hat r}\times \mathbf \Phi_{lm} = -\mathbf \Psi_{lm} $$ $$ \mathbf{\hat r}\times \mathbf \Psi_{lm} = \mathbf \Phi_{lm} $$ $$ \mathbf{\hat r}\times \mathbf Y_{lm} = 0 $$ I get from $\mathbf{\hat r}\times \mathbf E=0$ at the boundary and the independence of the harmonic functions: $$ \frac{g_l}{r}=\frac{dg_l}{dr}+\frac{g_l}{r}=0 $$ which gives $\frac{dg_l}{dr}=0$.
Note that you second boundary condition, $\mathbf{\hat r}\cdot \mathbf B = 0$ using: $$ \mathbf{\hat r}\cdot \mathbf \Phi_{lm} =\mathbf{\hat r}\cdot \mathbf \Psi_{lm} = 0 $$ $$ \mathbf{\hat r}\cdot \mathbf Y_{lm} = Y_{lm} $$ gives only $g_l=0$ at the boundary.
There is no condition on $f_l$ based on what you've given us. Are you sure you gathered all the information, and you don't have a typo? What is your reference?
Hope this helps and tell me if something's not clear.