I understand the proof for the expectation value of $\langle\hat{Q}\rangle$, which is shown for reference. Note, if you are already familiar with this proof then there is no need to read the contents of the quote below:
$$\langle\hat{Q}\rangle=\int_{x=-\infty}^{\infty}\psi^*\hat Q\,\psi \,dx=\int_{-\infty}^{\infty}\left(\sum\limits_m a_m\phi_m\right)^*\hat Q\color{blue}{\sum\limits_n a_n\phi_n}\,dx\tag{1}$$ where the parts marked blue hold because arbitrary $\psi$ can be written as a sum of eigenvalues $\phi_n$ with expansion coefficients $a_n$ such that
$$\psi=\color{blue}{\sum\limits_n a_n\phi_n}$$ Continuing with the proof from $(1)$:
$$\int_{-\infty}^{\infty}\left(\sum\limits_m a_m\phi_m\right)^*\hat Q\sum\limits_n a_n\phi_n\,dx=\int_{-\infty}^{\infty}\sum\limits_m {a_m}^*{\phi_m}^*\sum\limits_n a_n\,\color{red}{q_n\phi_n}\,dx\tag{2}$$
where the red part arises from the fact that $\hat Q$ satisfies the eigenvalue equation: $$\hat Q\,\phi_n=\color{red}{q_n\phi_n}$$ where $q_n$ is the associated eigenvalues with the operator $\hat Q$.Continuing the proof from $(2)$:
$$\int_{-\infty}^{\infty}\sum\limits_m {a_m}^*{\phi_m}^*\sum\limits_n a_n\,q_n\phi_n\,dx=\sum\limits_{m,n}{a_m}^*a_n\,q_n\color{#180}{\int_{-\infty}^{\infty}{\phi_m}^*\phi_n \,dx}\tag{3}$$
where the summations have been taken out the integral since they have no $x$ dependence such that the integrand comprises solely of eigenfunctions (that depend on $x$). Also, the part marked green is the Kronecker delta $\color{#180}{\large\delta_{m n}}$.So, resuming the proof from $(3)$:
$$\sum\limits_{m,n}{a_m}^*a_n\,q_n\large\delta_{m n}=\sum\limits_n {a_n}^*a_n\,q_n=\sum\limits_n|a_n|^2q_n\tag{4}$$
Since $$\large\delta_{m n}=\begin{cases}0 & m\ne n \\ 1 & m=n\end{cases}$$
Therefore we conclude from $(4)$ that
$$\langle\hat{Q}\rangle=\sum\limits_n|a_n|^2q_n$$
I fully understand the proof that $$\langle\hat{Q}\rangle=\sum\limits_n|a_n|^2q_n$$
So now I need to proof that $$\Big\langle{\hat{Q}}^2\Big\rangle=\sum\limits_n|a_n|^2{q_n}^2$$
Here is my attempt:
I know from the eigenvalue equation for $$\hat Q\,\phi_n=q_n\phi_n\tag{A}$$ So if I operate $\hat Q$ on both sides of $(\mathrm{A})$:
$$\hat Q\left(\hat Q\phi_n\right)=\hat Q \left(q_n\phi_n\right)$$
$$\implies {\hat Q}^2\phi_n=q_n\hat Q\phi_n$$
$$\implies {\hat Q}^2\phi_n={q_n}^2\phi_n\qquad\text{using the RHS of (A)}$$
From the chain rule $$\frac{du}{dx}=\frac{dy}{dx}\frac{du}{dy}$$ operating on some function $u$, I understand from this post that the chain rule can be written as $$\frac{d}{dx}=\frac{dy}{dx}\frac{d}{dy}$$ even in the absence of a function $u$ for which it can operate on.
So does this mean I can write $${\hat Q}^2\phi_n={q_n}^2\phi_n\tag{B}$$ as
$${\hat Q}^2={q_n}^2$$ where I have effectively ignored the eigenfunctions $\phi_n$.
If I now replace $\hat{Q}$ with ${\hat{Q}}^2$ in $$\langle\hat{Q}\rangle=\sum\limits_n|a_n|^2q_n$$
I will get $$\big\langle \hat{Q}^2\big\rangle=\sum\limits_n |a_n|^2{q_n}^2\tag{C}$$ as required.
Is this proof valid and if it's not valid what is the correct way to derive $(\mathrm{C})$?
EDIT:
In response to the first comment which mentions that with arbitrary operator $\hat O$ $$\langle\hat{O}\rangle=\sum\limits_n|a_n|^2o_n$$ for any operator so I can simply put $$\hat O=\hat Q^2$$ and somehow $$\langle\hat{Q}^2\rangle=\sum\limits_n|a_n|^2{q_n}^2$$ I do not understand the logic of this approach. More specifically why does letting $$\hat O=\hat Q^2$$ allow me to write $$\langle\hat{Q}^2\rangle=\sum\limits_n|a_n|^2{q_n}^2\,?$$ If anyone could please elaborate by giving some more details to this I will be most grateful, thanks.
EDIT 2:
I have been given a good answer now by 2 users especially AccidentalFourierTransform. I acknowledge now that eigenfunctions can never be ignored like I was doing before. I know that the $|a_n|^2$ are probabilities, but does the $a_n$ have some particular relationship to the eigenvalue equation $${\hat Q}^2\phi_n={q_n}^2\phi_n\,?$$ If so, this could be the source of my confusion and I would like to know how they are related?
Thanks again.
Best Answer
Let $O$ be an arbitrary operator with eigenvectors $|n\rangle$ and eigenvalues $o_n$: $$ O|n\rangle\equiv o_n|n\rangle \tag{1} $$
Then, as proven in the OP, the expectation value of $O$ is $$ \langle O\rangle=\sum_n |a_n|^2 o_n \tag{2} $$
For example, if $O=Q$ then $o_n=q_n$ and $$ \langle Q\rangle=\sum_n |a_n|^2 q_n\tag{3} $$
On the other hand, if $O=Q^2$ then $o_n=q_n^2$ and $$ \langle Q^2\rangle=\sum_n |a_n|^2 q^2_n \tag{4} $$
For generally, if $O=f(Q)$ then $o_n=f(q_n)$ and $$ \langle f(Q)\rangle=\sum_n |a_n|^2 f(q_n) \tag{5} $$