How degenerate time independent perturbation theory works? Let's consider the usual setup for time independent perturbation theory:
$$H=H_0+\varepsilon H'$$
and we can then set up the usual expansion:
$$(H_0+\varepsilon H')[|n_0\rangle+\varepsilon |n_1\rangle+\varepsilon ^2 |n_2\rangle+...]=(E_n^{(0)}+\varepsilon E_n^{(1)}+\varepsilon ^2 E_n^{(2)}+...)[n_0\rangle+\varepsilon |n_1\rangle+\varepsilon ^2 |n_2\rangle+...]$$
Long story short: when we have to solve a problem using perturbation theory the only thing we are interested about is how to calculate the corrections to the eigenstates and the eigenvalues.
If we are in the case of time independent non degenerate perturbation theory then this task is preatty straight forward once you know the formulas for the corrections:
$$E^{(k)}_n=\langle n_0|H'|n_{k-1}\rangle$$
$$|n_k\rangle=\frac{1}{H_0+E^{(0)}_n}|_{|n_0\rangle}[(E_n^{(1)}-H')|n_{k-1}\rangle+E_n^{(2)}|n_{k-2}\rangle+.....+E_n^{(k)}|n_0\rangle]$$
Done! Wonderful! But of course what if our Hamiltonian is degenerate? On textbooks I have found reasons why the old formulas do not work. I also understood that in some cases the perturbation cancels the degeneracy and in some other cases it doesn't. And there are also talks about the need to diagonalize the matrix in the degenerate space (this last point is not clear to me at the moment). Ok. But in practice: How can I set up and solve the perturbative expansion in the degenerate case? What are the formulas for the correction? (Knowing why the formulas work would also be nice but it is not the main point of this question)
Those are simple questions but I cannot seem to find any direct answer in my books or lecture notes. I would like a nice and concise answer. This topic seems really complicated to me as a beginner and I would like a summary of what is going on here. Especially from a practical point of view, on how can we solve exercises and expansions in the degenerate case.
 A: The main idea behind perturbation theory for degenerate states is to find not only corrections but also the states that are being corrected. Only specific states would acquire small corrections, others will be corrected by $O(1)$ terms. Let's consider as simple example. Consider a two level system given by the following Hamiltonian
\begin{equation}
H = 
\left (
\begin{array}{ccc}
m & \varepsilon \\
\varepsilon & m
\end{array}
\right ),
\end{equation}
with $\varepsilon \ll m$. The system can be solved exactly giving
\begin{equation}
E_\pm=m\pm\varepsilon ~~ \text{and} ~~ 
| \psi_\pm \rangle= \left (
\begin{array}{ccc}
1 \\
\pm 1
\end{array}
\right ).
\end{equation}
Now imagine we tried getting this result using perturbation theory. The unperturbed Hamiltonian is
\begin{equation}
H = 
\left (
\begin{array}{ccc}
m & 0 \\
0 & m
\end{array}
\right ),
\end{equation}
has degenerate eigenstates
\begin{equation}
| \psi^{(0)} \rangle = 
c_1\left (
\begin{array}{c}
1 \\
0
\end{array}
\right )
+
c_2\left (
\begin{array}{c}
0 \\
1
\end{array}
\right ),
\end{equation}
all with energy $E^{(0)}=m$. It is clear that only if you choose your unperturbed states to be
\begin{equation}
| \psi^{(0)}_{1,2} \rangle = \left (
\begin{array}{ccc}
1 \\
\pm 1
\end{array}
\right )
\end{equation}
corrections due to the perturbation is small (in this case it vanishes). How could we obtain that result without solving the system exactly? For that you are choosing an arbitrary basis for the unperturbed system $| \varphi_i \rangle$ and express the "true" unperturbed (and perturbed) eigenstates as linear combinations of those:
\begin{equation}
| \psi^{(0)}_i \rangle = c^{(0)}_{ij} | \varphi_j \rangle,
~~ \text{and} ~~
| \psi^{(1)}_i \rangle = c^{(1)}_{ij} | \varphi_j \rangle.
\end{equation}
Then multiplying the Schrödinger equation
\begin{equation}
(H_0+\varepsilon V) \left ( | \psi^{(0)}_i \rangle + \varepsilon | \psi^{(1)}_i \rangle \right )=
(E^{(0)}+\varepsilon E^{(1)}_i) \left ( | \psi^{(0)}_i \rangle + \varepsilon | \psi^{(1)}_i \rangle \right )
\end{equation}
by $\langle \phi_k |$ one gets
\begin{equation}
\sum_{j}\langle \varphi_k | V | \varphi_j \rangle c_{ij}^{(0)} = E_i^{(1)} c_{ik}^{(0)}.
\end{equation}
Omitting the index $i$ we see that these equations are nothing else but equations for eigenstates
\begin{equation}
\sum_j V_{k j} c_j = E^{(1)}c_k, 
\end{equation}
which implies that $\det (V-E^{(1)})=0$. From this equation $E_i^{(1)}$ and $c_{ij}^{(0)}$ are derived simultaneously.
Back to our example, we can chose
\begin{equation}
| \varphi_1 \rangle =
\left (
\begin{array}{ccc}
1 \\
0
\end{array}
\right ),
~~ \text{and} ~~
| \varphi_2 \rangle =
\left (
\begin{array}{ccc}
0 \\
1
\end{array}
\right ).
\end{equation}
The Schrödinger equation becomes
\begin{equation}
\left (
\begin{array}{cc}
m & \varepsilon \\
\varepsilon & m
\end{array}
\right )
\left (
\begin{array}{ccc}
c_{i1}^{(0)}+\varepsilon c_{i1}^{(1)} \\
c_{i2}^{(0)}+\varepsilon c_{i2}^{(1)}
\end{array}
\right ) 
= 
\left ( m+\varepsilon E_i^{(1)} \right )
\left (
\begin{array}{ccc}
c_{i1}^{(0)}+\varepsilon c_{i1}^{(1)} \\
c_{i2}^{(0)}+\varepsilon c_{i2}^{(1)}
\end{array}
\right ),
\end{equation}
or after simplification
\begin{equation}
\varepsilon\left (
\begin{array}{ccc}
c_{i2}^{(0)} \\
c_{i1}^{(0)}
\end{array}
\right ) 
= 
\varepsilon E_i^{(1)}
\left (
\begin{array}{ccc}
c_{i1}^{(0)} \\
c_{i2}^{(0)}
\end{array}
\right ),
\end{equation}
whose solution is
\begin{equation}
E^{(1)}=\pm 1, ~~ \text{for} ~~ 
\left (
\begin{array}{ccc}
1 \\
\pm 1
\end{array}
\right ),
\end{equation}
which is exactly what we had before.
A: What you are interested in is called secular equation.
The classical source is the second volume of Landau&Lifshitz https://books.google.ru/books?id=neBbAwAAQBAJ&pg=PA110&hl=ru&source=gbs_selected_pages&cad=2#v=onepage&q&f=false
Let $\psi_{n}^{(0)}, \psi_{n^{'}}^{(0)}$ be the eigenfunctions, belonging to the same eigenvalue $E_n^{(0)}$. By $\psi_{n}^{(0)}, \psi_{n^{'}}^{(0)}$ we assume unperturbed functions, selected in some arbitrary way. The correct eigenfunction in zeroth order are linear combinations of form:
$$
c_{n}^{(0)} \psi_{n}^{(0)} + c_{n^{'}}^{(0)} \psi_{n^{'}}^{(0)} + \ldots
$$
The substitution in the first order of perturbation for the energy $E_n^{(0)} + E^{(1)}$ into the second equation in your post gives:
$$
E^{(1)} c_{n}^{(0)} = \sum_{n^{'}} H_{n n^{'}} c_{n^{'}}^{(0)} 
$$
Or rewrite it in a following manner:
$$
\sum_{n^{'}} (H_{n n^{'}} - E^{(1)} \delta_{n n^{'}})c_{n^{'}}^{(0)} = 0
$$
This equation has solutions, as a system with zero right hand side, only if the matrix, defining the system is degenerate. For the square matrix it is equivalent to the vanishing of determinant:
$$
\boxed{\det(H_{n n^{'}} - E^{(1)} \delta_{n n^{'}}) = 0}
$$
This equation is the aforementioned secular equation. And the eigenvalue $E^{(1)}$ of the perturbation determines the energy correction, and the solutions of the equation the coefficients  $c_{n^{'}}^{(0)}$.
A: It is possible to set up an expansion for the degenerate case but only if you use the “right” basis.  The “right” basis is this basis that diagonalises the perturbation in the degenerate subspace of interest.  Then by construction there will be no off-diagonal terms in this subspace, i.e. in this new basis
with basis vectors $\vert\alpha_i\rangle$ so that $\hat V\vert\alpha_i\rangle=\lambda_i\vert\alpha_i\rangle$, you have
$\langle \alpha _k\vert \hat V\vert \alpha_j\rangle=\delta_{kj}$ so you never divide by $0$ since the expansion does not include terms where $k=j$.
If you use this new basis then you may proceed as if the problem was not degenerate.  The procedure can still fail if the perturbation $\hat V$ has repeated eigenvalues in the degenerate subspace of interest; in this case there’s nothing to be done, i.e no obvious perturbative expansion will exists for those remaining degenerate states.
