Chemistry - How to find the second order perturbation to wave function?
Solution 1:
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 &= \lambda^0W_0 + \lambda^1W_1 + \lambda^2W_2 + \cdots &&= \sum_{i=0}\lambda^i W_i\tag3\\ \Psi_\mathrm{ST} &= \lambda^0\Psi_0 + \lambda^1\Psi_1 + \lambda^2\Psi_2 + \cdots &&= \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&:& \mathbf{H}_0\Psi_0 &= W_0\Psi_0\tag{8a}\\ \lambda^1&:& \mathbf{H}_0\Psi_1 + \mathbf{H'}\Psi_0 &= W_0\Psi_1 + W_1\Psi_0\tag{8b}\\ \lambda^2&:& \mathbf{H}_0\Psi_2 + \mathbf{H'}\Psi_1 &= W_0\Psi_2 + W_1\Psi_1 + W_2\Psi_0\tag{8c}\\ \lambda^n&:& \mathbf{H}_0\Psi_n + \mathbf{H'}\Psi_{n-1} &= \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 &= \Phi_0 & \Psi_1 &= \sum_i c_i \Phi_i & \Psi_2 &=\sum_i d_i \Phi_i & \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 &= \langle\Phi_0|W_0|\Psi_0\rangle\\ \langle\Phi_0|\mathbf{H}_0|\Psi_0\rangle &= W_0\cancelto{1}{\langle\Phi_0|\Psi_0\rangle}\tag{9a}\\ W_0 &= \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 &= \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 &= 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 &= W_0\cdot c_0 + W_1\tag{9b}\\ W_1 &= \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 &= \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 &= \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 &= \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 &= \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 &= \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 &= \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 &= \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 &= \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 &= \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 &= W_2 \cancelto{0}{\langle\Phi_j|\Psi_0\rangle }\\ &= 0 \end{align}
Let's continue on the right hand side: \begin{align} \langle\Phi_j|W_1|\Psi_1\rangle &= W_1 \langle\Phi_j|\Psi_1\rangle\\ &= W_1 \sum_i c_i \langle\Phi_j|\Phi_i\rangle\\ &= 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)\\ &= c_j W_1 \end{align}
And the same for the last contribution: \begin{align} \langle\Phi_j|W_0|\Psi_2\rangle &= W_0 \langle\Phi_j|\Psi_2\rangle\\ &= W_0 \sum_i d_i \langle\Phi_j|\Phi_i\rangle\\ &= 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)\\ &= d_j E_0 &\mathrm{9a':}&(W_0 = E_0) \end{align}
Let's go to the left side: \begin{align} \langle\Phi_j|\mathbf{H}_0|\Psi_2\rangle &= \sum_i d_i \langle\Phi_j|\mathbf{H}_0|\Phi_i\rangle \\ &= \sum_i d_i E_i \langle\Phi_j|\Phi_i\rangle\\ &= 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}\\ &= d_j E_j \end{align}
And the last one can barely be simplified: \begin{align} \langle\Phi_j|\mathbf{H'}|\Psi_1\rangle &= \sum_i c_i \langle\Phi_j|\mathbf{H'}|\Phi_i\rangle\\ &= \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 &= \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 &= c_j W_1 + d_j E_0\\ d_j (E_j - E_0) &= 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 &= 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 &= \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)}\\ &= \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 &= \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.
Solution 2:
Since I had troubles with the notation the first time around, I did some real quantum chemistry (with a pencil and some paper) and finally was able to derive the whole RSPT ansatz. Since this is equivalent to the first approach, but yet also standalone, I decided to add it as a separate answer. In a few key steps, this approach is slightly different, so a one to one transformation will not be possible. $$%Introducing some shortcuts \require{cancel} \newcommand{\op}[1]{\hat{\mathrm{#1}}} %\op{H} \newcommand{\order}[1]{^{(#1)}} %E_n\order{1} \newcommand{\overlap}[3]{\mathcal{#1}_{#2}\order{#3}} %\overlap{S}{m}{1} \newcommand{\integral}[3]{\mathcal{#1}_{#2,#3}} %\integral{V}{i}{j} \newcommand{\tagref}[1]{\mathrm{#1}} $$
We start again by setting up the Hamiltonian, as a reference operator with the perturbation operator. $$\op{H} = \op{H}\order{0} +\lambda\op{V} \tag1$$
The perturbed Schrödinger equation therefore becomes $$(\op{H}\order{0} +\lambda\op{V})|n\rangle = E_n|n\rangle, \tag2$$ and we again develop the terms of the energy and wave function into Taylor power series. \begin{align} E_n &= E_n\order{0} + \lambda E_n\order{1} + \lambda^2 E_n\order{2} + \dots & &= \sum_i \lambda^i E_n\order{i} \tag3\\ |n\rangle &= |n\order{0}\rangle + \lambda |n\order{1}\rangle + \lambda^2 |n\order{2}\rangle + \dots & &= \sum_i \lambda^i |n\order{i}\rangle \tag4\\ \end{align}
We can write the unperturbed system, that provides us with a complete solution and therefore with a complete set of energies and wave functions, i.e. $n\in\mathbb{N}$. $$\op{H}\order{0} |n\order{0}\rangle = E_n\order{0} |n\order{0}\rangle \tag5$$
We also know from this, that all solutions are orthonormal, i.e. \begin{align} \langle n\order{0} |n\order{0}\rangle &= 1, \langle n\order{0} |m\order{0}\rangle &= 0. \tag6 \end{align}
Further we choose our perturbed wave function to be (intermediate) normalised, i.e. \begin{align} \langle n\order0|n\rangle &= 1, & \langle n|n\rangle &= 1, & \langle n\order{i}|n\order{j}\rangle &= \delta_{ij} \begin{cases} 1, & i=j\\ 0, & i\neq j\\ \end{cases}.\tag7 \end{align}
We continue with gathering the equations for all our orders. \begin{align} \lambda^0 &:& \op{H}\order0 |n\order0\rangle &= E_n\order0 |n\order0\rangle\tag{8a}\\ \lambda^1 &:& \op{H}\order0 |n\order1\rangle + \op{V} |n\order0\rangle &= E_n\order0 |n\order1\rangle + E_n\order1 |n\order0\rangle\tag{8b}\\ \lambda^2 &:& \op{H}\order0 |n\order2\rangle + \op{V} |n\order1\rangle &= E_n\order0 |n\order2\rangle + E_n\order1 |n\order1\rangle E_n\order1 |n\order0\rangle\tag{8c}\\ \vdots\\ \lambda^n &:& \op{H}\order0 |n\order{n}\rangle + \op{V} |n\order{n-1}\rangle &= \sum_{i=0}^n E_n\order{i} |n\order{n-i}\rangle\tag{8d}\\ \end{align}
The two fundamental tricks of algebra are adding zero and multiplying by one.
We use this to set up our bases, we basically expand our bases into an auxiliary base (also known as 'resolution-of-the-identity'), which is possible, since the unperturbed system provides the complete set. This is what makes this approach somewhat different. While in my other answer we are just expanding the wafer function into an ordinary boring linear combination, we are using here the complete Hilbert space to do the same thing, but a little more elegant. However, this is also the key step: Multiplying by one.
At the same time we introduce an overlap integral $\integral{S}{m}{i}$, which makes writing a little easier further on. For any order $i\neq0$ we can therefore write:
\begin{align}
|n\order{i}\rangle &=
\sum_{m} |m\order0\rangle\langle m\order0|n\order{i}\rangle\\
&=
\sum_{m\neq n} |m\order0\rangle
\underbrace{\langle m\order0|n\order{i}\rangle}_{\overlap{S}{m}{i}}
+ |n\order0\rangle\cancelto{0}{\langle n\order0|n\order{i}\rangle}\\
|n\order{i}\rangle &= \sum_{m\neq n} |m\order0\rangle \overlap{S}{m}{i}\\
\end{align}
The bases we need are for first and second order, therefore we write \begin{align} |n\order{1}\rangle &= \sum_{m\neq n} |m\order0\rangle \overlap{S}{m}{1}, & |n\order{2}\rangle &= \sum_{m\neq n} |m\order0\rangle \overlap{S}{m}{2}, \end{align} which seems rather trivial.
The next step is getting the expressions for the energies and corrections to the wave function.
For $\tagref{8a}$ this is trivial and it just yields the expectation value of the unperturbed Hamiltonian. We project on (multiply from the left with) $\langle n\order0|$, hence
\begin{align}
\langle n\order0|E_n\order0 |n\order0\rangle &=
\langle n\order0|\op{H}\order0 |n\order0\rangle \\
E_n\order0 \cancelto{1}{\langle n\order0|n\order0\rangle} &=
\langle n\order0|\op{H}\order0 |n\order0\rangle \\
E_n\order0 &= \langle n\order0|\op{H}\order0 |n\order0\rangle. \tag{9a}\\
\end{align}
For the first order correction to the energy we use the hermiticity of the Hamiltonian, i.e. we rewrite $\tagref{8a}$ in bra instead of the ket notation we used all way through. $$\langle n\order0| \op{H}\order0 = E_n\order0 \langle n\order0|$$ We then further go about and project $\tagref{8b}$ on $\langle n\order0|$, hence \begin{align} \langle n\order0|\left(\op{H}\order0|n\order1\rangle + \op{V}|n\order0\rangle \right) &= \langle n\order0|\left(E_n\order0|n\order1\rangle + E_n\order1|n\order0\rangle \right)\\ \langle n\order0|\op{H}\order0|n\order1\rangle + \langle n\order0| \op{V}|n\order0\rangle &= \langle n\order0|E_n\order0|n\order1\rangle + \langle n\order0| E_n\order1|n\order0\rangle\\ E_n\order0 \cancelto{0}{\langle n\order0|n\order1\rangle} + \langle n\order0| \op{V}|n\order0\rangle &= E_n\order0\cancelto{0}{\langle n\order0|n\order1\rangle} + E_n\order1 \cancelto{1}{\langle n\order0| n\order0\rangle}\\ E_n\order1 &= \langle n\order0| \op{V}|n\order0\rangle \tag{9b} \end{align}
Now for the first order correction to the waffle function. I have to admit, this was the part where I struggled most. We will project on $\langle k\order0|\neq\langle n\order0|$, therfore obtaining: \begin{align} \langle k\order0|\left(\op{H}\order0|n\order1\rangle + \op{V}|n\order0\rangle \right) &= \langle k\order0|\left(E_n\order0|n\order1\rangle + E_n\order1|n\order0\rangle \right)\\ \langle k\order0|\op{H}\order0|n\order1\rangle + \langle k\order0| \op{V}|n\order0\rangle &= \langle k\order0|E_n\order0|n\order1\rangle + \langle k\order0| E_n\order1|n\order0\rangle\\ E_k\order0 \langle k\order0|n\order1\rangle + \langle k\order0| \op{V}|n\order0\rangle &= E_n\order0\langle k\order0|n\order1\rangle + E_n\order1 \cancelto{0}{\langle k\order0| n\order0\rangle}\\ E_k\order0 \langle k\order0|\sum_{m\neq n} |m\order0\rangle \overlap{S}{m}{1} + \langle k\order0| \op{V}|n\order0\rangle &= E_n\order0\langle k\order0|\sum_{m\neq n} |m\order0\rangle \overlap{S}{m}{1}\\ \end{align} Before it gets to messy, we are going to reduce each of the terms individually and then reinsert, to obtain our expression. \begin{align} E_k\order0 \langle k\order0|\sum_{m\neq n} |m\order0\rangle \overlap{S}{m}{1} &= \sum_{m\neq n,k}E_k\order0 \cancelto{0}{\langle k\order0|m\order0\rangle} \overlap{S}{m}{1} + E_k\order0 \cancelto{1}{\langle k\order0|k\order0\rangle} \overlap{S}{k}{1}\\ &= E_k\order0 \overlap{S}{k}{1}\\ \end{align} The second term stays the same, and the last one simplifies analogously. \begin{align} E_n\order0 \langle k\order0|\sum_{m\neq n} |m\order0\rangle \overlap{S}{m}{1} &= \sum_{m\neq n,k}E_n\order0 \cancelto{0}{\langle k\order0|m\order0\rangle} \overlap{S}{m}{1} + E_n\order0 \cancelto{1}{\langle k\order0|k\order0\rangle} \overlap{S}{k}{1}\\ &= E_n\order0 \integral{S}{k}{1}\\ \end{align} And reinsterting we can write and resubstitute for the overlap integral $\overlap{S}{k}{1}$, with $m=k$: \begin{align} E_k\order0 \overlap{S}{k}{1} + \langle k\order0| \op{V}|n\order0\rangle &= E_n\order0 \overlap{S}{k}{1}\\ \integral{S}{k}{1} &= \frac{\langle k\order0| \op{V}|n\order0\rangle}{E_n\order0 -E_k\order0}\\ |n\order1\rangle &= \sum_{k\neq n} |k\order0\rangle \frac{\langle k\order0| \op{V} |n\order0\rangle}{E_n\order0 -E_k\order0} \tag{10b}\\ \end{align}
To make matters a little bit easier further on, we will introduce another short notation for the perturbed integrals: $$\integral{V}{i}{j} := \langle i\order0|\op{V}|j\order0\rangle$$ And therefore we can write \begin{align} |n\order1\rangle &= \sum_{k\neq n} |k\order0\rangle \frac{\integral{V}{k}{n}}{E_n\order0 -E_k\order0} \tag{10b}\\ \end{align}
For the second order correction to the energy we project $\tagref{8c}$ on $\langle n\order0|$. This is again fairly straight forward, the key is to insert $\tagref{10b}$ at the end to obtain the expression. \begin{align} \langle n\order0|\left(\op{H}\order0|n\order2\rangle + \op{V}|n\order1\rangle\right) &= \langle n\order0|\left( E_n\order0|n\order2\rangle + E_n\order1|n\order1\rangle + E_n\order2|n\order0\rangle\right)\\ \langle n\order0|\op{H}\order0|n\order2\rangle + \langle n\order0|\op{V}|n\order1\rangle &= E_n\order0 \langle n\order0|n\order2\rangle + E_n\order1 \langle n\order0|n\order1\rangle + E_n\order2 \langle n\order0|n\order0\rangle\\ \small E_n\order0 \cancelto{0}{\langle n\order0|n\order2\rangle} + \langle n\order0|\op{V}|n\order1\rangle &= \small E_n\order0 \cancelto{0}{\langle n\order0|n\order2\rangle} + E_n\order1 \cancelto{0}{\langle n\order0|n\order1\rangle} + E_n\order2 \cancelto{1}{\langle n\order0|n\order0\rangle}\\ E_n\order2 &= \langle n\order0|\op{V}|n\order1\rangle\\ % &= \langle n\order0|\op{V} \sum_{k\neq n} |k\order0\rangle % \frac{\langle k\order0| \op{V} % |n\order0\rangle}{E_n\order0 -E_k\order0}\\ &= \sum_{k\neq n} \frac{\langle n\order0|\op{V} |k\order0\rangle \langle k\order0| \op{V} |n\order0\rangle}{E_n\order0 -E_k\order0}\tag{9c}\\ &= \sum_{k\neq n} \frac{\integral{V}{n}{k} \integral{V}{k}{n}}{E_n\order0 -E_k\order0} \end{align}
Now let's get to the real fancy part: The second order correction to the basil. We project once again $\tagref{8c}$ on $\langle k\order| \neq \langle n\order 0|$. We do of course need our expression for the first order base and our expansion for the second order. But let us take one step at the time. \begin{align} \langle k\order0|\left( \op{H}\order0|n\order2\rangle + \op{V}|n\order1\rangle\right) &= \langle k\order0 \left( E_n\order0|n\order2\rangle + E_n\order1|n\order1\rangle + E_n\order2|n\order0\rangle\right)\\ \small \langle k\order0| \op{H}\order0|n\order2\rangle + \langle k\order0| \op{V}|n\order1\rangle &= \small \langle k\order0| E_n\order0|n\order2\rangle + \langle k\order0| E_n\order1|n\order1\rangle + \langle k\order0| E_n\order2|n\order0\rangle\\ \end{align} And we proceed again term by term, because it is pretty messy. We are starting with the right hand side, and there with the third term, as this one cancels, because of our normalisation. \begin{align} \langle k\order0| E_n\order2|n\order0\rangle &= E_n\order2 \cancelto{0}{\langle k\order0| n\order0\rangle} =0\\ \end{align} For the second term on the right we need our definition of the base. Since we are projecting on $k$ we choose the base to be represented by $m$. We also need the expression for the first order correction to the energy $\tagref{9b}$. Also note, that we are using the contracted notation. \begin{align} \langle k\order0| E_n\order1|n\order1\rangle &= E_n\order1 \langle k\order0|n\order1\rangle\\ &= E_n\order1 \langle k\order0| \sum_{m\neq n} |m\order0\rangle \frac{\integral{V}{m}{n}}{E_n\order0 - E_m\order0}\\ &= \sum_{m\neq n,k} E_n\order1 \cancelto{0}{\langle k\order0| m\order0\rangle} \frac{\integral{V}{m}{n}}{E_n\order0 - E_m\order0} + E_n\order1 \cancelto{1}{\langle k\order0| k\order0\rangle} \frac{\integral{V}{k}{n}}{E_n\order0 - E_k\order0}\\ &= \frac{\integral{V}{n}{n}\integral{V}{k}{n}}{E_n\order0 - E_k\order0}\\ \end{align}
For the first term on the right hand side we just need our expansion of the second order correction of the base. \begin{align} \langle k\order0| E_n\order0|n\order2\rangle &= E_n\order0 \langle k\order0|n\order2\rangle\\ &= E_n\order0 \langle k\order0| \sum_{m\neq n} |m\order0\rangle \overlap{S}{m}{2}\\ &= \sum_{m\neq n,k} E_n\order0 \cancelto{0}{\langle k\order0|m\order0\rangle} \overlap{S}{m}{2} + E_n\order0 \cancelto{1}{\langle k\order0|k\order0\rangle} \overlap{S}{k}{2}\\ &= E_n\order0 \overlap{S}{k}{2}\\ \end{align}
We have to use the same for the first term on the left hand side. In addition we are using the bra notation for the energy to simplify. \begin{align} \langle k\order0| \op{H}\order0|n\order2\rangle &= E_k\order0 \langle k\order0|n\order2\rangle\\ &= E_k\order0 \langle k\order0|\sum_{m\neq n} |m\order0\rangle \overlap{S}{m}{2} \\ &= \sum_{m\neq n,k} E_k\order0 \cancelto{0}{\langle k\order0|m\order0\rangle} \overlap{S}{m}{2} + E_k\order0 \cancelto{1}{\langle k\order0|k\order0\rangle} \overlap{S}{k}{2}\\ &= E_k\order0 \overlap{S}{k}{2}\\ \end{align}
For the second term on the left side we once again need the expression for the first order correction to bass, choosing $m$, as we are already projecting on $k$. \begin{align} \langle k\order0| \op{V}|n\order1\rangle &= \langle k\order0| \op{V} \sum_{m\neq n}|m\order0\rangle \frac{\integral{V}{m}{n}}{E_n\order0 -E_m\order0}\\ &= \sum_{m\neq n}\frac{\integral{V}{k}{m} \integral{V}{m}{n}}{E_n\order0 -E_m\order0}\\ \end{align}
Now, let's put all the stuffing back into the turkey. Then we rearrange to complete bird and solve for our desired expression for the base. \begin{align} E_k\order0 \overlap{S}{k}{2} + \sum_{m\neq n}\frac{\integral{V}{k}{m} \integral{V}{m}{n}}{E_n\order0 -E_m\order0} &= E_n\order0 \overlap{S}{k}{2} + \frac{\integral{V}{n}{n}\integral{V}{k}{n}}{E_n\order0 - E_k\order0}\\ (E_n\order0 - E_k\order0)\overlap{S}{k}{2} &= \sum_{m\neq n}\frac{\integral{V}{k}{m} \integral{V}{m}{n}}{E_n\order0 -E_m\order0} - \frac{\integral{V}{n}{n}\integral{V}{k}{n}}{E_n\order0 - E_k\order0}\\ \overlap{S}{k}{2} &= \sum_{m\neq n}\frac{\integral{V}{k}{m} \integral{V}{m}{n}}{(E_n\order0 -E_m\order0)(E_n\order0 - E_k\order0)} - \frac{\integral{V}{n}{n}\integral{V}{k}{n}}{(E_n\order0 - E_k\order0)^2}\\ |n\order2\rangle &= \sum_{k\neq n} \sum_{m\neq n} |k\order0\rangle \frac{\integral{V}{k}{m} \integral{V}{m}{n}}{(E_n\order0 -E_m\order0)(E_n\order0 - E_k\order0)}\\ %split into two lines &\phantom{=~} - \sum_{k\neq n} |k\order0\rangle \frac{\integral{V}{k}{n}\integral{V}{n}{n}}{(E_n\order0 - E_k\order0)^2}\\ \end{align}
And last, but not least, let's remove that nasty shorthand and make everything pretty again. \begin{align} |n\order2\rangle &= \sum_{k\neq n} \sum_{m\neq n} |k\order0\rangle \frac{ \langle k\order0|\op{V}|m\order0\rangle \langle m\order0|\op{V}|n\order0\rangle }{(E_n\order0 -E_m\order0)(E_n\order0 - E_k\order0)} - \sum_{k\neq n} |k\order0\rangle \frac{ \langle k\order0|\op{V}|n\order0\rangle \langle n\order0|\op{V}|n\order0\rangle }{(E_n\order0 - E_k\order0)^2}\\ \end{align}
And there we are. Pretty.
Now you will probably have noticed, that you were one step away from the ultron goal. In your last equation you simply had to reduce to the non-zero parts.
\begin{equation} ( E_n^0 - E_k^0 ) \langle k^0 |n^2\rangle = \sum_{m \neq n} \frac{ \langle k^0| \hat{V} |m^0 \rangle \langle m^0 | \hat{V} | n^0 \rangle}{E_n^0 - E_m^0} - \langle n^0 | \hat{V} | n^0 \rangle \sum_{m \neq n} \frac{\langle k^0 |m^0 \rangle \langle m^0 | \hat{V} | n^0 \rangle}{E_n^0 - E_m^0} \end{equation}
\begin{align} ( E_n^0 - E_k^0 ) \langle k^0 |n^2\rangle &= \sum_{m \neq n} \frac{ \langle k^0| \hat{V} |m^0 \rangle \langle m^0| \hat{V} |n^0 \rangle }{E_n^0 - E_m^0}\\ &\phantom{=~} - \sum_{m \neq n,k} \frac{ \cancelto{0}{\langle k^0 |m^0 \rangle} \langle n^0 | \hat{V} | n^0 \rangle \langle m^0 | \hat{V} | n^0 \rangle }{E_n^0 - E_m^0}\\ &\phantom{=~} - \frac{ \cancelto{1}{\langle k^0 |k^0 \rangle} \langle k^0 | \hat{V} | n^0 \rangle \langle n^0 | \hat{V} | n^0 \rangle }{E_n^0 - E_k^0}\\ ( E_n^0 - E_k^0 ) \langle k^0 |n^2\rangle &= \sum_{m \neq n} \frac{ \langle k^0| \hat{V} |m^0 \rangle \langle m^0| \hat{V} |n^0 \rangle }{E_n^0 - E_m^0} - \frac{ \langle k^0 | \hat{V} | n^0 \rangle \langle n^0 | \hat{V} | n^0 \rangle }{E_n^0 - E_k^0}\\ \langle k^0 |n^2\rangle &= \sum_{m \neq n} \frac{ \langle k^0| \hat{V} |m^0 \rangle \langle m^0| \hat{V} |n^0 \rangle }{(E_n^0 - E_m^0)(E_n^0 - E_k^0)} - \frac{ \langle k^0 | \hat{V} | n^0 \rangle \langle n^0 | \hat{V} | n^0 \rangle }{\left(E_n^0 - E_k^0\right)^2}\\ \end{align}
Now you just have to expand into the complete base by multiplying with $\sum_{k\neq n}|k^0\rangle$ and you get to the same expression as above.
And as a theoretician a proof is only complete with a coffee stain, and trust me, there have been many.
Image is courtesy of Roger Karlsson (http://www.free-photo-gallery.org/photos/coffee-stain/) obtained from flicker.