Solution (source code)

= Solution

Put $\gamma=2m\hbar\omega$. If $\psi_n$ is a normalized <energy eigenstate> at $E_n$, the factorization from part (a) gives
$$
\|a_+\psi_n\|^2=\langle\psi_n,a_-a_+\psi_n\rangle=\gamma(n+1).
$$
It follows inductively that $v_n=(a_+)^n\psi_0$ has squared <norm> $\gamma^n n!$. The <commutator> with $H$ establishes its <energy eigenvalue> at each step. Hence
$$
\boxed{\psi_n=\frac{(a_+)^n\psi_0}{(2m\hbar\omega)^{n/2}\sqrt{n!}}}
$$
is normalized and has $E_n=(n+\tfrac12)\hbar\omega$. Distinct levels are <orthogonal> because $H$ is a <self-adjoint operator>. In particular,
$$
a_+\psi_n=\sqrt{\gamma(n+1)}\,\psi_{n+1},\qquad a_-\psi_n=\sqrt{\gamma n}\,\psi_{n-1}.
$$

To decompose the specified initial <wavefunction>, retain this exact phase convention for the <creation and annihilation operators>. Their difference gives
$$
X=\frac{a_+-a_-}{2im\omega}.
$$
Since $a_-\psi_0=0$ and $a_-a_+\psi_0=\gamma\psi_0$,
$$
\begin{aligned}
X^2\psi_0&=-\frac1{4m^2\omega^2}(a_+^2-a_+a_--a_-a_++a_-^2)\psi_0\\
&=\frac{\hbar}{2m\omega}\left(\psi_0-\sqrt2\,\psi_2\right).
\end{aligned}
$$
This <quadratic excitation of an oscillator ground state>, $\Psi(x,0)=C_0^{-1}X^2\psi_0$, contains only the $n=0$ and $n=2$ <energy eigenstates>. The printed initial function has an unspecified overall normalization, so the <Born rule> requires division by its squared <norm>. After <wavefunction normalization> the state is
$$
\frac1{\sqrt3}\psi_0-\sqrt{\frac23}\psi_2.
$$
The minus sign is consistent with the paper's convention: its $a_+$ is $i$ times the usual real-coordinate raising operator multiplied by $\sqrt\gamma$, so its $\psi_2$ has the opposite sign to the usual real second <harmonic oscillator eigenfunction>. The <energy measurement> probabilities are independent of this choice of phase.

Therefore the complete answer is
$$
\boxed{\Pr\left(E=\frac12\hbar\omega\right)=\frac13,\qquad \Pr\left(E=\frac52\hbar\omega\right)=\frac23.}
$$
\b[Every other energy has probability zero.] At any later time the two coefficients only acquire their separate <dynamical phases>, so these <probabilities> remain the same. A direct check of the normalization is
$$
\int_{-\infty}^{\infty}x^4e^{-\alpha x^2}\,dx=\frac{d^2}{d\alpha^2}\sqrt{\frac\pi\alpha}=\frac{3\sqrt\pi}{4\alpha^{5/2}},\qquad \alpha=\frac{m\omega}{\hbar}.
$$
This also confirms the factor of three in the squared <norm> of the two-component expansion.