As part of a larger problem, I am trying to find the average squared Hamiltonian of a system with eigenfunctions $\psi_{1,1}$, $\psi_{1,2}$, $\psi_{2,1}$, $\psi_{2,2}$ and the following wave function:
$$ \Psi\left(\mathbf{r};t=0\right)=c\sum_{j=1}^2 \psi_{ij}\left(\mathbf{r}\right) $$
The problem defines the following operators:
\begin{align*} \hat{H}\psi_{ij} &= iE\psi_{ij} \\ \hat{Q}\psi_{ij} &= jQ\psi_{ij} \end{align*}
where $ \{i,j\}\in\mathbb{R} $. I have already calculated that
\begin{align*} p\left(\mathbf{r}\right) &= c\,\left(\langle\psi_{i1}|\psi_{i1}\rangle + \langle\psi_{i1}|\psi_{i2}\rangle + \langle\psi_{i2}|\psi_{i1}\rangle + \langle\psi_{i2}|\psi_{i2}\rangle\right) \\ 1 &= c\,\left(1 + 0 + 0 + 1\right) \\ c &= \frac{1}{2} \end{align*}
and
\begin{align*} \langle\Psi|\hat{H}|\Psi\rangle &= \langle\psi_{i1}|\hat{H}|\psi_{i1}\rangle + \langle\psi_{i1}|\hat{H}|\psi_{i2}\rangle + \langle\psi_{i2}|\hat{H}|\psi_{i1}\rangle + \langle\psi_{i2}|\hat{H}|\psi_{i2}\rangle \\ &= \frac{i}{2}\left(E_{i1}\langle\psi_{i1}|\psi_{i1}\rangle + E_{i2}\langle\psi_{i1}|\psi_{i2}\rangle + E_{i2}\langle\psi_{i2}|\psi_{i1}\rangle + E_{i2}\langle\psi_{i2}|\psi_{i2}\rangle\right) \\ &= \frac{i}{2}\left(E_{i1}\langle\psi_{i1}|\psi_{i1}\rangle + 0 + 0 + E_{i2}\langle\psi_{i2}|\psi_{i2}\rangle\right) \\ &= \frac{i}{2}\left(E_{i1}+E_{i2}\right) \end{align*}
However, I'm not quite sure how to scale it up to $ \langle\Psi|\hat{H}^2|\Psi\rangle $. I am putting forward the assumption that $ \langle\hat{H}^2\rangle - \langle\hat{H}\rangle^2 $ should = 0 since I am working with eigenfunctions, and that therefore
\begin{align*} \langle\Psi|\hat{H}^2|\Psi\rangle &= \,?!? \\ &= \left(\frac{i}{2}\left(E_{i1}+E_{i2}\right)\right)^2 \\ &= \frac{i^2}{4}\left(E_{i1}^2 + E_{i1}E_{i2} + E_{i2}^2\right) \end{align*}
But I am unsure how to prove it, hence the ?!?.