Perturbed wavefunctions I've been studying a number of TISE perturbation problems, where the Hamiltonian is $H = H_{0} + \epsilon H^{\prime}$, the wave function for bound state $n$ is $|n\rangle = \sum_{m=0}^{\infty} \epsilon^{m} |n^{(m)}\rangle$, the energy corresponding to $|n\rangle$ is $E_{n} = \sum_{m=0}^{\infty} \epsilon^{m} E_{n}^{(m)}$, $\epsilon$ is the expansion/bookkeeping parameter, and $m$ is the expansion order.  The perturbation procedure is straightforward.  For each $n$: 1) substitute $|n\rangle$ and $E_{n}$ into the perturbed TISE; 2) group together all terms for each $m$; 3) solve for the $|n^{(m)}\rangle$ and $E_{n}^{(m)}$ in terms of all $|p^{(0)}\rangle \neq |n^{(0)}\rangle$ unperturbed wave functions.
The perturbed wave functions are not guaranteed to be normalized, which is easily remedied: $|n\rangle \, \rightarrow |n\rangle / \, C_{n}$, where $C_{n}^{2} = \langle n|n\rangle$.  I noticed that the perturbed wave functions are not orthonormal, i.e., $\langle n|p\rangle \neq 0$ in general, suggesting that the matrix elements of the perturbation component of the Hamiltonian are not symmetric.
What does the non-orthonormality mean physically?  Do the perturbed wave functions span the complete space of wave functions, i.e., can a weighted sum represent any wave function in the space?  If I do something to make them orthonormal, e.g., Gram-Schmidt, will there still be a 1-to-1 correspondence between wave functions and energies (adiabatic changes)? I suspect that the answer is yes for small perturbations, but I am unsure for large perturbations.  What does it mean when the 1-to-1 correspondence disappears?
EDIT: OK, Let's try to address Bob's comments.  Hopefully, it will lead to increased understanding on my part ...
I was originally considering all orders, but let's do what Bob suggests and keep only the zeroth and first orders.  For $\epsilon = 1$, the wave function is $|n\rangle = |n^{(0)}\rangle + |n^{(1)}\rangle$.  Taking the inner product of two states yields $\langle m|n\rangle = \langle m^{(0)}|n^{(0)}\rangle + \langle m^{(0)}|n^{(1)}\rangle + \langle m^{(1)}|n^{(0)}\rangle + \langle m^{(1)}|n^{(1)}\rangle$.
When $m = n$, the first term on the right is unity (assuming normalization of unperturbed states), and the second and third terms are zero by definition due to the perturbation/normalization process).  The last term is small and real.  Normalization of $|n\rangle$ is trivial.
When $m \neq n$, the first term of the right is identically zero, since the unperturbed wave functions are defined as orthogonal.  I have not been able to prove that the last three terms on the right are zero or cancel.  The same statement holds for higher orders only more so.  In other words, just because $\langle n^{(1)}|n^{(0)}\rangle$ is zero doesn't necessarily mean that $\langle n^{(1)}|m^{(0)}\rangle$ or $\langle m^{(1)}|n^{(0)}\rangle$ are zero.
The $m = n$ and $m \neq n$ cases above can be derived using the expansion
$|n^{(1)}\rangle = \sum_{m\neq n} |m^{(0)}\rangle \langle m^{(0)}|n^{(1)}\rangle$ $=$ $\sum_{m\neq n} |m^{(0)}\rangle \frac{\langle m^{(0)}|H_{i}|n^{(0)}\rangle}{E_{n}^{(0)}-E_{m}^{(0)}}$.
Am I missing something?  It certainly is possible.
 A: OK, for $m \neq n$,
$$\langle m|n\rangle = \langle m^{(0)}|n^{(0)}\rangle + \langle m^{(0)}|n^{(1)}\rangle + \langle m^{(1)}|n^{(0)}\rangle + \langle m^{(1)}|n^{(1)}\rangle, $$
but as you recognized the first term  is identically zero; moreover we ignore the last term to this order, since it is of second order, and should only be considered in conjunctions with 2nd order wavefunction corrections.
How do you see the 2nd term cancels the third?
Since 
$$|n^{(1)}\rangle =  \sum_{k\neq n} |k^{(0)}\rangle \frac{\langle k^{(0)}|H_{i}|n^{(0)}\rangle}{E_{n}^{(0)}-E_{k}^{(0)}},$$
(Note it is rotten karma to use an active index as a dummy variable), we dot on the left by $\langle m^{(0)}|$, so the sum collapses to just one term, by the original orthonormality at zero order,
$$
\langle m^{(0)}|n^{(1)}\rangle= \frac{\langle m^{(0)}|H_{i}|n^{(0)}\rangle}{E_{n}^{(0)}-E_{m}^{(0)}}=-\langle n^{(0)}|m^{(1)}\rangle ~. 
$$
Observe how this expression managed to be antisymmetric in m,n, by dint of the denominator. So orthonormality to first order is established.
But this is a general result.
In fact, you can't avoid it. At some point in your course, you must have learned that for hermitian hamiltonians, 
$$
\langle m|H|n\rangle= E_n \langle m|n\rangle= E_m\langle m|n\rangle .  
$$ 
the middle and rightmost side are gotten by action on the ket and bra, respectively. But, for non degenerate energies $E_m\neq E_n$, this can only hold for orthogonal states, $\langle m|n\rangle=0$. So, if this failed in perturbation theory, there should be a consistency mistake around...
