How to find the second order perturbation to wave function?
Disclaimer This post is some kind of a legacy post. Find the notation used in the question in the other answer. I added this proof as I was not entirely certain I understood the notation correctly. As it turned out, I did not. As a result I posted the complete derivation of RSPT up to second order, trying to guide anyone through it using a different notation.$\require{cancel}$
Obviously, I will come back to the initial problem only at the end, if you are impatient, just skip ahead.
Please also note, that I removed the parts that became redundant when I posted the second answer. If you are interested in legacy versions, just look them up in the history.
Also, this page is heavily loaded with $\mathcal{M}^{ath}\mathrm{J_{a}X}$, give it time to load completely.



In Rayleigh Schrödinger perturbation theory, we set up the Hamiltonian as $$\mathbf{H}=\mathbf{H_0}+\lambda\mathbf{H'}.\tag1$$ 

The perturbed Schrödinger equation is therefore $$\mathbf{H}\Psi=W\Psi_\mathrm{ST}.\tag2$$
You can develop this into Taylor series:
\begin{align}
W                &amp;= 
     \lambda^0W_0    + \lambda^1W_1    + \lambda^2W_2    + \cdots 
     &amp;&amp;= \sum_{i=0}\lambda^i W_i\tag3\\
\Psi_\mathrm{ST} &amp;= 
     \lambda^0\Psi_0 + \lambda^1\Psi_1 + \lambda^2\Psi_2 + \cdots 
     &amp;&amp;= \sum_{i=0}\lambda^i \Psi_i\tag4\\
\end{align}

Your unperturbed reference system is a complete solution denoted as 
$$\mathbf{H}_0\Phi_i = E_0\Phi_i,\tag5$$
where $i\in\mathbb{N}$.

You chose your system to be intermediate normalised.
$$\langle\Psi_\mathrm{ST}|\Phi_0\rangle = 1\tag6$$
Since all your solutions of the unperturbed system are orthogonal,
$\langle\Phi_i|\Phi_j\rangle=\delta_{ij},$
the same holds for your perturbed system, hence
$$\langle\Psi_{i\neq0}|\Phi_0\rangle = 0.\tag7$$

Now we can combine $(1)$ through $(4)$ and gather them according to their order.
\begin{align}
\lambda^0&amp;:&amp;
  \mathbf{H}_0\Psi_0 
    &amp;= W_0\Psi_0\tag{8a}\\
\lambda^1&amp;:&amp; 
  \mathbf{H}_0\Psi_1 + \mathbf{H'}\Psi_0 
    &amp;= W_0\Psi_1 + W_1\Psi_0\tag{8b}\\
\lambda^2&amp;:&amp; 
  \mathbf{H}_0\Psi_2 + \mathbf{H'}\Psi_1 
    &amp;= W_0\Psi_2 + W_1\Psi_1 + W_2\Psi_0\tag{8c}\\
\lambda^n&amp;:&amp; 
  \mathbf{H}_0\Psi_n + \mathbf{H'}\Psi_{n-1} 
    &amp;= \sum_{i=0}^n W_i\Psi_{n-i}\tag{8d}\\
\end{align}

We do have a complete set of solutions, hence we can express perturbations as linear combinations of these.
\begin{aligned}
\Psi_0 &amp;= \Phi_0 &amp; 
\Psi_1 &amp;= \sum_i c_i \Phi_i &amp; 
\Psi_2 &amp;=\sum_i d_i \Phi_i &amp; 
\dots
\end{aligned}

Now we project equations $8$ onto $\Psi_0$. For $\mathrm{8a}$ this is trivial and it leads back to the expectation value of the energy for the unperturbed system. (In the last transformation we use the linear combinations.)
\begin{align}
\langle\Phi_0|\mathbf{H}_0|\Psi_0\rangle 
    &amp;= \langle\Phi_0|W_0|\Psi_0\rangle\\
\langle\Phi_0|\mathbf{H}_0|\Psi_0\rangle 
    &amp;= W_0\cancelto{1}{\langle\Phi_0|\Psi_0\rangle}\tag{9a}\\
W_0 &amp;= \langle\Phi_0|\mathbf{H}|\Phi_0\rangle = E_0\tag{9a'}
\end{align}
The next step is similar, and everything else will follow from there analogously. For $\mathrm{8b}$ we obtain:
\begin{align}
\langle\Phi_0|\mathbf{H}_0|\Psi_1\rangle + 
\langle\Phi_0|\mathbf{H'}|\Psi_0\rangle 
    &amp;=   \langle\Phi_0|W_0|\Psi_1\rangle 
       + \langle\Phi_0|W_1|\Psi_0\rangle\\
\langle\Phi_0|\mathbf{H}_0|\Psi_1\rangle + 
\langle\Phi_0|\mathbf{H'}|\Psi_0\rangle 
    &amp;= W_0\cancelto{c_0}{\langle\Phi_0|\Psi_1\rangle} +
       W_1\cancelto{1}{\langle\Phi_0|\Psi_0\rangle}\\
\cancelto{c_o\cdot E_0}{\sum_i c_i \langle\Phi_0|\mathbf{H}_0|\Phi_i\rangle} 
       + \langle\Phi_0|\mathbf{H'}|\Psi_0\rangle 
    &amp;= W_0\cdot c_0 + W_1\tag{9b}\\
W_1 &amp;= \langle\Phi_0|\mathbf{H'}|\Psi_0\rangle\tag{9a'}
\end{align}

Similarly we can obtain the coefficients by projecting on some other function than $\Phi_0$. I will cut this a little bit shorter.
\begin{align}
\langle\Phi_j|\mathbf{H}_0|\Psi_1\rangle 
     + \langle\Phi_j|\mathbf{H'}|\Psi_0\rangle 
    &amp;= \langle\Phi_j|W_0|\Psi_1\rangle + 
       \langle\Phi_j|W_1|\Psi_0\rangle\\
\cancelto{c_j\cdot E_j}{\langle\Phi_j|\mathbf{H}_0|\Psi_1\rangle} 
     + \langle\Phi_j|\mathbf{H'}|\Psi_0\rangle 
    &amp;= \cancelto{c_j\cdot E_0}{\langle\Phi_j|W_0|\Psi_1\rangle} 
     + \cancelto{0}{\langle\Phi_j|W_1|\Psi_0\rangle}\tag{10b}\\
c_j &amp;= \frac{\langle\Phi_j|\mathbf{H'}|\Psi_0\rangle}{E_0-E_j}\tag{10b'}
\end{align}

Now this is impossible to solve for $c_0$, but we still have the intermediate normalisation.
\begin{align}
0 &amp;= \langle\Phi_0|\Psi_1\rangle
   = \sum_i c_i \langle\Phi_0|\Phi_i\rangle
   = c_0\cancelto{1}{\langle\Phi_0|\Phi_i\rangle} + 
     \sum_{i\neq0} c_i \cancelto{0}{\langle\Phi_0|\Phi_i\rangle} = c_0\\
\end{align}

Similarly this applies to all zeroth coefficients, $c_0 = d_0 = \dots = 0$.

Let's get back to business and continue with $\mathrm{8c}$. 
\begin{align}
\langle\Phi_0|\mathbf{H}_0|\Psi_2\rangle 
  + \langle\Phi_0|\mathbf{H'}|\Psi_1\rangle 
  &amp;= \langle\Phi_0|W_0|\Psi_2\rangle 
   + \langle\Phi_0|W_1|\Psi_1\rangle 
   + \langle\Phi_0|W_2|\Psi_0\rangle\\
\scriptsize\cancelto{d_0\cdot E_0}{\langle\Phi_0|\mathbf{H}_0|\Psi_2\rangle} 
   + \sum_i c_i\langle\Phi_0|\mathbf{H'}|\Phi_i\rangle 
  &amp;= \scriptsize W_0\cancelto{d_0=0}{\langle\Phi_0|\Psi_2\rangle} 
   + W_1\cancelto{c_0=0}{\langle\Phi_0|\Psi_1\rangle} 
   + W_2\cancelto{1}{\langle\Phi_0|\Psi_0\rangle}\tag{9c}\\
W_2 &amp;= \sum_i c_i\langle\Phi_0|\mathbf{H'}|\Phi_i\rangle\tag{9c'}
\end{align}

We do now insert $\mathrm{10b'}$ and get the second order correction:
\begin{align}
W_2 &amp;= \sum_{i\neq0}\frac{
             \langle\Phi_0|\mathbf{H'}|\Phi_i\rangle
             \langle\Phi_i|\mathbf{H'}|\Psi_0\rangle
            }{E_0-E_i}\tag{9c''}
\end{align}

Now we have the expression for the second order energy, let's do the messy stuff for the second order contribution to the wave function. 

Skip here for the short version, which is still very long, but that's how we roll for the second order correction to the wave function

We start by multiplying $\mathrm{8c}$ with $\Phi_j$ from the left and integrating (projection). I will include a few more steps for this problem, hopefully that will help transferring it to your notation.
\begin{align}
\langle\Phi_j|\mathbf{H}_0|\Psi_2\rangle +
\langle\Phi_j|\mathbf{H'}|\Psi_1\rangle &amp;=
\langle\Phi_j|W_0|\Psi_2\rangle +
\langle\Phi_j|W_1|\Psi_1\rangle +
\langle\Phi_j|W_2|\Psi_0\rangle\\
\end{align}

I will go through this term by term this time, as it is quite a long equation.
Let's start with the most simplest term, the third one on the right hand side. As you noted, it cancels.
\begin{align}
\langle\Phi_j|W_2|\Psi_0\rangle 
   &amp;= W_2 \cancelto{0}{\langle\Phi_j|\Psi_0\rangle }\\
   &amp;= 0
\end{align} 

Let's continue on the right hand side:
\begin{align}
\langle\Phi_j|W_1|\Psi_1\rangle 
   &amp;= W_1 \langle\Phi_j|\Psi_1\rangle\\
   &amp;= W_1 \sum_i c_i \langle\Phi_j|\Phi_i\rangle\\
   &amp;= W_1 \left(c_j\cancelto{1}{\langle\Phi_j|\Phi_j\rangle} +  
         \sum_{i\neq j} c_i \cancelto{0}{\langle\Phi_j|\Phi_i\rangle}\right)\\
   &amp;= c_j W_1
\end{align} 

And the same for the last contribution:
\begin{align}
\langle\Phi_j|W_0|\Psi_2\rangle 
   &amp;= W_0 \langle\Phi_j|\Psi_2\rangle\\
   &amp;= W_0 \sum_i d_i \langle\Phi_j|\Phi_i\rangle\\
   &amp;= W_0 \left(d_j\cancelto{1}{\langle\Phi_j|\Phi_j\rangle} 
        + \sum_{i\neq j} d_i \cancelto{0}{\langle\Phi_j|\Phi_i\rangle}\right)\\
   &amp;= d_j E_0 &amp;\mathrm{9a':}&amp;(W_0 = E_0)
\end{align} 

Let's go to the left side:
\begin{align}
\langle\Phi_j|\mathbf{H}_0|\Psi_2\rangle 
   &amp;= \sum_i d_i \langle\Phi_j|\mathbf{H}_0|\Phi_i\rangle \\
   &amp;= \sum_i d_i E_i \langle\Phi_j|\Phi_i\rangle\\
   &amp;= d_j E_j \cancelto{1}{\langle\Phi_j|\Phi_j\rangle}
      + \sum_{i\neq j} d_i E_i \cancelto{0}{\langle\Phi_j|\Phi_i\rangle}\\
   &amp;= d_j E_j
\end{align}

And the last one can barely be simplified:
\begin{align}
\langle\Phi_j|\mathbf{H'}|\Psi_1\rangle
  &amp;= \sum_i c_i \langle\Phi_j|\mathbf{H'}|\Phi_i\rangle\\
  &amp;= \cancelto{0}{c_0 \langle\Phi_0|\mathbf{H'}|\Phi_0\rangle}
   + \sum_{i\neq 0} c_i \langle\Phi_j|\mathbf{H'}|\Phi_i\rangle\\ 
\end{align}

Now, Avengers reassemble, rearrange, and insert the expression for $W_1$:

\begin{align}
\langle\Phi_j|\mathbf{H}_0|\Psi_2\rangle +
\langle\Phi_j|\mathbf{H'}|\Psi_1\rangle &amp;=
\langle\Phi_j|W_0|\Psi_2\rangle +
\langle\Phi_j|W_1|\Psi_1\rangle +
\langle\Phi_j|W_2|\Psi_0\rangle\\
d_j E_j + \sum_{i\neq 0} c_i \langle\Phi_j|\mathbf{H'}|\Phi_i\rangle
 &amp;= c_j W_1 + d_j E_0\\
d_j (E_j - E_0)
 &amp;= c_j \langle\Phi_0|\mathbf{H'}|\Psi_0\rangle
   - \sum_{i\neq 0} c_i \langle\Phi_j|\mathbf{H'}|\Phi_i\rangle\\
d_j &amp;=
   c_j \frac{\langle\Phi_0|\mathbf{H'}|\Psi_0\rangle}{(E_j - E_0)} 
   - \sum_{i\neq 0} c_i \frac{\langle\Phi_j|\mathbf{H'}|\Phi_i\rangle}{(E_j - E_0)}
\end{align}

Now introduce the first order base, $c_i$ and $c_j$ from $\mathrm{10b'}$, there will be some minus shuffling as well.
\begin{align}
d_j 
  &amp;=
   \frac{\langle\Phi_j|\mathbf{H'}|\Psi_0\rangle}{(E_0 - E_j)} 
   \frac{\langle\Phi_0|\mathbf{H'}|\Psi_0\rangle}{(E_j - E_0)} 
   - \sum_{i\neq 0} \frac{\langle\Phi_i|\mathbf{H'}|\Psi_0\rangle}{(E_0 - E_i)}
     \frac{\langle\Phi_j|\mathbf{H'}|\Phi_i\rangle}{(E_j - E_0)}\\
  &amp;=
   \frac{\langle\Phi_j|\mathbf{H'}|\Psi_0\rangle}{(E_0 - E_j)} \cdot
   \frac{\langle\Phi_0|\mathbf{H'}|\Psi_0\rangle}{(-1)(E_0 - E_j)} 
   - \sum_{i\neq 0} \frac{\langle\Phi_i|\mathbf{H'}|\Psi_0\rangle}{(E_0 - E_i)}
     \frac{\langle\Phi_j|\mathbf{H'}|\Phi_i\rangle}{(-1)(E_0 - E_j)}\\
 d_j &amp;=
   \sum_{i\neq 0} 
   \frac{\langle\Phi_j|\mathbf{H'}|\Phi_i\rangle%
         \langle\Phi_i|\mathbf{H'}|\Psi_0\rangle%
        }{(E_0 - E_j)(E_0 - E_i)} 
   -\frac{\langle\Phi_j|\mathbf{H'}|\Psi_0\rangle%
          \langle\Phi_0|\mathbf{H'}|\Psi_0\rangle}{(E_0 - E_j)^2} \tag{10c'}
\end{align}

Well, I am officially exhausted right now. 



And lastly, a small remark on your notation (I keep this paragraph for nostalgic reasons). It took me quite some time understanding it, and as we saw (in some earlier version of this answer), I was mistaken.
While I do not mind the hat notation for operators, indicating order on top of the quantities is problematic in my opinion. To me this place is reserved for powers. I suggest you change it to one of the following: $E_{n,0},\ {}^{(0)}\!E_n,\ E^{(0)}_n$. When you have a look at the other answer, you will notice, I have decided to stick to the last option. But as I do not know which book your are using it is also not easy to recommend. I believe there are as many different notations as there are books out there. I personally preferred the book "Introduction to Computational Chemistry, 2nd Edition"by Frank Jensen.
As we settled in the comments, I appreciate that you stick to Dirac's original notation. And also, notation is like Marvel and DC, you pick one and stick with it, you can still appreciate the other, but just not as much. Or a football (Americans pronounce this as "soccer") club, or ...
As long as you understand it and can explain it to someone else, you are fine.
And now we all deserve some chocolate and some coffee and cake. Thank you for reading this post until the end.