Why are the zeroth order terms in degenerate perturbation theory the eigenstates of the perturbing Hamiltonian? I have for quite some time now tried to find a satisfactory answer to this, but I haven't yet. In perturbation theory, with small parameter $\lambda$, we expand the eigenstate as
$$| E \rangle=| E^{(0)} \rangle + \lambda | b \rangle + ...$$
Where $| E^{(0)} \rangle$ is an eigenstate of the unperturbed Hamiltonian. The problem is that when there is degeneracy, there is a choice of these eigenstates. I know the answer is that they are chosen to be eigenvectors of the perturbing Hamiltonian, but my question is why.
 A: This is really a great question.
Take a look at the figure attached here. The circle indicates a 2-dimensional degenerate subspace at $\lambda=0$.
In red we indicate two possible basis states for the subspace. In blue we show another possible choice of basis states.
Now look at the curvy green lines, these are the states as they evolve for $\lambda > 0$. The green curves do not connect with the red states at $\lambda=0$. The perturbation series is a Taylor series, which is a continuous function (it's a polynomial), so there is no way to make the series go from the red states at $\lambda=0$ to the green states for $\lambda > 0$.
Clearly, we have to start the perturbation series at $\lambda=0$ with the blue states, because those connect up with the green ones when $\lambda > 0$. The green states are, by definition, eigenstates of the perturbed Hamiltonian, so the blue states at $\lambda=0$ must also be eigenstates of the perturbed Hamiltonian.

Degenerate perturbation theory states. We show the states as functions of the perturbation parameter $\lambda$ (green), and various choices of the unperturbed states (blue, red).
A: As I cannot comment (reputation>50 needed), I put my comment here, as this might be considered an answer, too.
The answer and picture of DanielSank is very nice, so just to add what I think is a satisfactory answer as well, just in words.
On page 139 of my favourite textbook:-), the explanation is right in the beginning of the "Secular equation" paragraph.
"We denote by $\psi_n^{(0)}$, $\psi_{n'}^{(0)},\dots$ the eigenfunctions belonging to the same eigenvalue $E_n^{(0)}$ of the energy. The choice of these functions is, as we know, not unique; instead of them we can choose any $s$ (where $s$ is the degree of degeneracy of the level $E_n^{(0)}$) independent linear combinations of these functions. The choice ceases to be arbitrary, however, if we subject the wave functions to the requirement that the change in them under the action of the small applied perturbation should be small.
...
The coefficients in these combinations are determined, together with the
corrections in the first approximation to the eigenvalues, as follows."
Then it turns out that the coefficients are obtained as a solution of the characteristic (secular) equation and that such eigenvectors are the eigenvectors of the perturbed Hamiltonian.
