Wave function decomposition Problem:
Given the wave function $\Psi_0=A\sin^2(\theta)$ along with the Hamiltonian operator of a physical system: $H=\frac{L^2}{2I}+g B L_z$,


*

*find the eigenvalues and eigenfunctions of $\hat{H}$ and 

*derive the form of $\Psi(t)$.


Solution:
I apply the $\hat{U}(t)$ time-evolution operator on $\Psi_0$ to get $\Psi(t)$: $\Psi(t)=\hat{U}(t)\Psi_0=\exp\left(-\frac{i}{\hbar} H t\right) \Psi_0$.
In order to calculate the effect of $\hat{H}$ on $\Psi_0$, I need to write $\Psi_0$ as a linear combination of the spherical harmonics $Y_{\ell m}$.
I have a list of the first few (up to $\ell=2$) spherical harmonics in my textbook, but I am not so sure on how to proceed. The best guess I can think of is: $Y_{20}=\sqrt{\frac{5}{16\pi}}(3\cos^2(\theta)-1)$ if I consider that $\cos^2(\theta)=1-\sin^2(\theta)$.
Then $\Psi_0 = \frac{2}{3}A - \left(\frac{A}{3}\sqrt{\frac{16\pi}{5}} \right)Y_{20}$.
My questions are these (edited more than once):


*

*How do I tackle the constant term ? Do I use $Y_{00}=\frac{1}{2}\sqrt{\frac{1}{\pi}}$ for that?

*Is it physically meaningful to have $Y_{\ell m}$ terms with different $\ell$ in the same wave function, like e.g. $Y_{00} + Y_{20}$ ? I am asking because until now I was encountering only terms with the same $\ell$ but different $m$. I presume it is.

*Is there anything fundamentally wrong with my strategy?

*How do I calculate constant A? I used the normalization condition $|\Psi_0|^2=1$, and since $\hat{U}(t)$ is unitary, $\Psi(t)$ will also be normalized.
Finally, $\Psi_0 = \sqrt{\frac{5}{6}}Y_{00}-\sqrt{\frac{1}{6}}Y_{20}$. At this point I have the following question:

*How do I calculate $\exp\left(-\frac{i}{\hbar}Ht\right)$ say on $Y_{20}$? I recollect a theorem which states that if operator $\hat{A}$ has eigenvalues $\alpha$, then operator $f(\hat{A})$ has the same eigenfunctions with eigenvalues $f(\alpha)$ ($f$ needs to fullfil some criteria that escape me). Anyway, the "problem" is that in this case $\hat{H}$ is a function of both $L^2$ and $L_z$.
Given that $Y_{\ell m}$ constitute a common set of eigenfunctions for both $L^2$ and $L_z$, could I claim that $Y_{20}$ is also an eigenfunction of $H(L,L_z)$ with eigenvalues $\frac{\hbar^2 \ell (\ell+1)}{2I} + gB(m\hbar), m \in Z$ ?
EDIT: I use this Mathematica snippet to verify my results (up until now):
$\text{In:=} \text{FullSimplify}\left[\frac{4}{3} \sqrt{\pi } A Y_{00}-\frac{A}{3} \sqrt{\frac{16 \pi }{5}} Y_{20}\text{ /.}\,\\ \left\{Y_{20}\to \sqrt{\frac{5}{16 \pi }} \left(3 \text{Cos}^2(t)-1\right),Y_{00}\to 
\frac{1}{2}\sqrt{\frac{1}{\pi }}\right\}\right]\\
\text{Out=}A-A \text{Cos}^2[t]$
 A: *

*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.
A: Let
$$P_{l}^{m}=|Y_{l}^{m}><Y_{l}^{m}|.$$
Then
$$H=\sum_{l,m}\{\frac{l(l+1)}{2I}+gBm\}P_{l}^{m},\;\hat{U}(t)=\sum_{l,m}%
\exp[-\frac{it}{\hslash}\{\frac{l(l+1)}{2I}+gBm\}]P_{l}^{m},$$
and
$$\Psi(t)=\sum_{l,m}\exp[-\frac{it}{\hslash}\{\frac{l(l+1)}{2I}+gBm\}]P_{l}%
^{m}\Psi_{0}
\\=\sum_{l,m}\exp[-\frac{it}{\hslash}\{\frac{l(l+1)}{2I}%
+gBm\}]<Y_{l}^{m}|\Psi_{0}>|Y_{l}^{m}>$$
Your remark:  " I recollect a theorem which states that if operator $\hat{A}$
has eigenvalues $\alpha$, then operator $f(\hat{A})$ has the same
eigenfunctions with eigenvalues $f(\alpha)$ " is completely correct  and for
general $\Psi_{0}$ you only need to calculate the overlap (inner product)
$<Y_{l}^{m}|\Psi_{0}>.$ 
$\Psi_{0}$ which are combinations of terms with
different $l$'s and $m$'s, say
$$\Psi_{0}=\sum_{l_{1},m_{1}}c_{l_{1}m_{1}}Y_{l_{1}}^{m_{1}},$$
result in
$$\Psi(t)=\sum_{l_{1},m_{1}}c_{l_{1}m_{1}}\exp[-\frac{it}{\hslash}\{\frac
{l_{1}(l_{1}+1)}{2I}+gBm_{1}\}]Y_{l_{1}}^{m_{1}}$$
since $<Y_{l}^{m}|Y_{l_{1}}^{m_{1}}>=\delta_{ll_{1}}\delta_{mm_{1}}$.
