Perturbation theory with degeneracy even after 1st order Most textbooks on basic quantum mechanics tell you that when your initial Hamiltonian $H_0$ has degenerate states, then before you can do (time independent) perturbation theory with a perturbation matrix $V$ on it, you have to first diagonalize $H_0 + V$ in the subspace of the degenerate states. 
Sounds nice enough, but what if $V$ is already diagonal in that subspace?
The example I have in mind is a hydrogen molecule in the tight-binding + Heitler-London limit, e.g. $U \gg t$. The Hamiltonian is 
$$\left(\begin{array}{cccc} U & 0 & t & t\\ 0 & U & t & t \\ t & t & 0 & 0 \\ t & t & 0 & 0\end{array}\right)$$
where the first two basis states are the ionized states (both electrons on the same atom) and the other two basis states are the covalent states.
Now, this Hamiltonian is easy enough to diagonalize directly, and if the new ground-state energy is expanded up to second order in $t$, we get $-4t^2/U$ as the drop in energy.
I now try to picture this in simple terms of second order perturbation theory, and there I only see a drop of $-2t^2/U$: Because from either of the unperturbed ground states (the covalent ones), I can have two second-order processes, with amplitudes $t^2$ each and with energy denominator $U$ each. Since I basically have a two-fold degenerate ground state, is this where the missing factor $2$ comes from? If so, what's the mechanism behind this?
 A: This answer is basically what I learned from the Lubos Motl's answer(hopefully I understand it right), with a little bit generalization.
Question raised
We want to solve a Hamiltonian $H=H_0+\lambda V$ perturbatively . What we know is the solution of $H_0$, whose eigenstates are $|n^0\rangle$ with eigenenergies $E_n^0$ has no degeneracy; it also has a degenerate state $|d_i^0\rangle$ with eigenenergy $E_d^0$, where $i=1,2,3\cdots$g, the degenerate vector space $D$ is thus spanned by $|d_i^0\rangle$ 

Solution
In general the degeneracy is lifted when we turn on the perturbation. 
$$
H|d_i\rangle=E_{di}|d_i\rangle
$$
However, we must keep in mind that when we switch off the perturbation, $d_i$ is not necessarily tend to $d_i^0$ as $\lambda\to 0$. Rather, it will tend to a state $|l_i^0\rangle$ which is a linear combination of $|d_i^0\rangle$. Clearly, the vector space spanned by $|l_i^0\rangle$ is also $D$.
$$
|d_i\rangle=|l_i^0\rangle+\lambda|d_i^1\rangle+\lambda ^2|d_i^2\rangle+\cdots\\
E_{di}=E_d^0+\lambda E_d^1+\lambda^2E_d^2+\cdots
$$
Substitute to the Shordinger eqution, we get:
First order:
$$
(H_0-E_d^0)|d_i^1\rangle=(E_d^1-V)|l_i^0\rangle
$$
Using $\langle d_j^0|$ to inner product, we get $\langle d_j^0|(E_d^1-V)|l_i^0\rangle=0$, i.e. vector $(E_d^1-V)|l_i^0\rangle$ has no component in the vector space $D$. We can then safely write :
$$\tag{1}
\begin{aligned}
|d_i^1\rangle&=\frac{1}{H_0-E_d^0}(E_d^1-V)|l_i^0\rangle\\
                     &=P\frac{1}{H_0-E_d^0}P(E_d^1-V)|l_i^0\rangle\\
                     &=-P\frac{1}{H_0-E_d^0}PV|l_i^0\rangle
\end{aligned}
$$
where $P=\sum_{n\notin D}|n^0\rangle\langle n^0|$ is the projection operator. Clearly, $|d_i^1\rangle$ also has no component in vector space $D$.
Also, let $|l_i^0\rangle=\sum_j c_{ij}|d_i^0\rangle$, we get the secular equation: $(E_d^1\delta_{ij}-V_{ij})c_{ij}=0$, solving this, we can get $E_d^1$ and $c_{ij}$ thus the vector $|l_i^0\rangle$ that diagnolize $V$.
Second order
$$
H_0|d_i^2\rangle+V|d_i^1\rangle=E_d^0|d_i^2\rangle+E_d^1|d_i^1\rangle+E_d^2|l_i^0\rangle
$$
Using bra $\langle l_i^0|$ to do the inner product, and provided $|l_i^0\rangle$ is normalized, we get:
$$\tag{2}
\begin{aligned}
E_{di}^2&=\langle l_i^0|V|d_i^1\rangle \\
         &=-\langle l_i^0|VP\frac{1}{H_0-E_d^0}PV|l_i^0\rangle \\
         &=-\sum_{n\notin D}\frac{|\langle l_i^0|V|n^0\rangle|^2}{E_n^0-E_d^0}
\end{aligned}
$$
A very important remark
Note that in eqn(2)  $|l_i^0\rangle$ is the state that  $| d_i(\lambda)\rangle$ evolve to when we turn off the perturbation gradually ($\lambda\to 0)$. As seen in first order perturbation , it is indeed the vectors that diagonalize $V$, and the corresponding $V_{ii}$ is just the first order energy. 
However, when the degeneracy is not lifted in the first order. The "good" ket  $|l_i^0\rangle$ does not seems so good, because $|d_i(\lambda)\rangle$ will still have ambiguity which final state to evolve as $\lambda\to 0$.
In this case, we need second order to lift the degeneracy, and  $\Lambda=VP\frac{1}{H_0-E_d^0}PV$ plays the role of $V$ as in the first order. the states that diagonalize $\Lambda$ should be considered as "good", you can construct a matrix of $g\times g$ dimension using $|d_i^0\rangle$ and obtain the good states $|k_i^0\rangle$ as a linear combination of them, the corresponding $\Lambda_{ii}=\langle k_i^0|\Lambda|k_i^0\rangle$ is thus the second order perturbation energy.
A: First, just to be sure about the answers to this particular problem: the eigenvalues of the $4\times 4$ matrix are 
$$0,\quad U\quad {\rm and}\quad U/2\pm \sqrt{(U/2)^2+4t^2}$$
When expanded to the first nontrivial order, the last two eigenvalues are 
$$ 0 - \frac{4t^2}U \quad {\rm and} \quad U+\frac{4t^2}U. $$
Note that the corrections to the energy arise at order $t^2$ so the first-order perturbation theory is not enough in this case.
Second, the problem that prevents us from choosing the right eigenstates by the simple method is - as you correctly notice - that the matrix $V$, i.e. the matrix multiplied by $t$, has a vanishing upper left $2\times 2$ block as well as the right lower $2\times 2$ block - both of these blocks vanish. 
So $V$ doesn't lift the degeneracy "inside the degenerate subspaces" only. This is, of course, related to the fact that the first-order $O(t)$ corrections to the energy eigenvalues vanish.
The standard formula of perturbation theory for the second-order corrections to energy is
$$ E_n = E_n^{(0)} + t \langle n^{(0)} |V|n^{(0)}\rangle + t^2 \sum_{k\neq n}\frac{\left| \langle k^{(0)} |V| n^{(0)} \right|^2}{E_n^{(0)}-E_k^{(0)}} +O(t^3)  $$
Now, the $t^2$ term should give us $\pm 4t^2/U$ if it works. And of course, it does as long as we choose the right superpositions as the zeroth-order eigenvectors.
In particular, let's just reveal the eigenvectors. The eigenvalue $U$ comes with the eigenvector $(1,-1,0,0)^T / \sqrt{2}$ and similarly the eigenvalue $0$ comes with the eigenvector $(0,0,1,-1)^T / \sqrt{2}$. The $U+4t^2/U$ comes from $(1,1,0,0)^T/\sqrt{2}$ and similarly $0-4t^2/U$ comes from $(0,0,1,1)^T/\sqrt{2}$. The transposition means that the vectors should be written in the conventional column form. I added the $1/\sqrt{2}$ factor to make all of them normalized - and they're orthogonal, too.
For each calculated state (and its second order correction to energy), there's just one nonzero term in the $\sum_{k\neq n}$ sum. It has the denominator - the energy difference - $U$ if we calculate $E_n$ near $U$ or $-U$ if we calculate $E_n$ near $0$. And it connects the states $(1,1,0,0)^T/\sqrt{2}$ with $(0,0,1,1)^T/\sqrt{2}$ or vice versa.
Note that the matrix element of the $V$ matrix with the two $((1,1),(1,1))$ off-diagonal blocks between the two states I mentioned at the end of the previous paragraph is $1/2$ (from the two $1/\sqrt{2}$ normalization factors) times $8$ (because there are eight nonzero entries $t$ in your matrix, or $1$ in mine, and each of them contributes the same one to the matrix element). 
So the matrix element is simply $4$ and the perturbation theory gives you the right $\pm 4t^2/E$ correction to the energy, from a single term, with the proper sign.
Now, the only step I haven't quite justified was the right choice of the eigenvectors - such as $(1,1,0,0)^T  / \sqrt{2}$. How could I have seen this was the right one? Well, in this case, it was the intuitively right choice. While the $2\times 2$ blocks on the diagonal preserved the degeneracy, the energy shift came from the $2\times 2$ off-block-diagonal blocks and those made it natural to use this basis.
More generally, if we're in a similar situation that the degeneracy isn't lifted by first-order corrections in $V$ in the subspaces, what we really have to diagonalize is $V (H_0-E_0)^{-1\prime} V$ where the prime indicates that one has to omit the (divergent) terms from the vanishing energy differences (of course, we need a sensible result). This is the operator whose expectation value de facto gives us the second-order energy correction.
In this particular case, this operator will have a block diagonal form, $U^{-1} {\rm diag} [ ((+1,+1),(+1,+1)),  ((-1,-1),(-1,-1)) ]$, and by diagonalizing, you get the right initial eigenvectors to deal with. Once you have the right eigenvectors to start with, their perturbations are infinitesimal at each order of the perturbation theory and the standard formulae of perturbation theory work without any extra subtleties, as the example above showed.
Again, the only thing one has to be careful about are the right zeroth-order initial eigenvectors. In a more generic case, they're given as eigenvectors of $V$ because $V$ lifts the degeneracy in each subspace. If it doesn't lift the degeneracy, only higher-order terms do so. But the operator $V (H-E_0)^{-1\prime} V$ plays the same role as $V$, and by diagonalizing it, we get the right initial eigenvectors.
A: Some remarks:
1) A 3-dimensional example is the matrix
 2 Εo         ΑΕ                     0

  Α*Ε          Εo                    2ΑΕ

    0            2Α*Ε                2 Εo

This also has two degenerate states that are not directly connected, only indirectly via the nondegenerate second state. In fact one of the degenerate states
is unshifted even with the full diagonalization
2)Diagonalising to find the eigenvalues: This is ok as a post-mortem, but kind of cheating when one wants to do the problem in perturbation theory: If one has done the
analysis to get the eigenvalues and eigenvectors of the full matrix, why do perturbation theory? Of course in the examples shown, it's actually easier  and faster to do the problem exactly than to do "post-mortem" analysis on what did not work and why
In other words the problem is: How do you go about it if you do absolutely want to use pert theory without first finding eigenvalues and eigenvectors of the full  matrix.
3) What one can easily see is that because the degenerate states are not connected by V or H' matrix elements, one needs to  involve states outside the degenerate block. I am not sure what the second order propagator term  solves in this case-the issue is not the operator, but the states used. Anyway, the point is that perturbation theory will not do  here. In my example, if one writes
|1'>=c11|1>+c21|2>+c31|3> and correspondingly for the rest
(Ε1 –2 Ε0 ) c11  = ΑΕ c21     (1)
(Ε3 –2 Ε0 ) c33  = 2ΑΕ* c23    (2)
 (Ε1 -ε2)c21=Η21’c11 +Η23’c33+  Η22’ c21 =>(Ε1 -Ε0 )c21=Η21’c11 +Η23’c33=>
c21= [ Η23’c33 +Η21’c11 ] /(Ε1 –Ε0 )  (3)
(Ε3 -ε2)c23=Η23’c33 +Η21’c11+  Η22’ c23 =>(Ε3 -Ε0)c23=Η23’c33 +Η21’c11=>
c23= [Η23’c33 +Η21’c11 ] /(Ε3 –Ε0 )  (4)
and eliminating c23 and c21 in (1) and (2) we get a homogeneous system
(Ε1 –2 Ε0 ) c11  = ΑΕ[ AEc33 + Α*Εc11 ] /(Ε1 –Ε0 )]
(Ε3 –2 Ε0 ) c33  = 2Α*Ε[ AEc33 + Α*Εc11 ] /(Ε3 –Ε0 ) 
Setting the determinant to 0 gets an equation connecting E1 and E3, but these are two unknowns. So bottom line: in such a case where degenerate states are not connected by the perturbation,  diagonalize the part consisting of these states plus the states that are "most connected" to them, i.e. have the strongest matrix elements
