Bloch theorem and Bloch states In deriving the Bloch theorem, they first give a Schrodinger equation,
$$\hat H\psi_n(r)=E_n\psi_n(r)$$  $n=0,1,2…$
$\hat T_R$ is a translation operator commuting with $\hat H$.
Its eigenvalue is
$$\hat T_R\psi_{k}(r)=e^{ikR}\psi_{k}(r)$$
$k=…$
For the reason that $\hat T_R$ commutes with $\hat H$, and they should have the same complete set of eigenvectors $\lbrace\psi_{nk}(r)\rbrace$. We should have
$$\hat H\psi_{nk}(r)=E_n\psi_{nk}(r)$$
$n=0,1,2…$
$k=…$
But this puzzles me, as we have $n$ eigenvalues for the Schrodinger equation at the beginning. For the Schrodinger equation itself I can't see where does these $k$ quantum number comes from. The set $\lbrace\psi_n(r)\rbrace$, with the number of $n$ eigenvectors, should itself be a complete set according to the Schrodinger equation, is it contradictory that we need $n*k$ number of eigenvectors $\lbrace\psi_{nk}(r)\rbrace$ to form a complete set?
At first, I thought that for a given $n$ value, $\psi_{nk}(r)$ might be degenerate for different $k$ value. But later, I realized that this is not true as a specific $n$ the eigenvalue (Energy) $\epsilon$ still change for different $k$. 

I think for the SEQ. $\hat H\psi_n(r)=E_n\psi_n(r)$   $n=0,1,2…$,
I agree with you that '$k$ labels the degeneracy of the energy eigenvalues.' as we have
$$\hat T_R\psi_{n}(r)=e^{ikR}\psi_{n}(r)$$ $k=k_1,k_2,...$.
This means that for each $k$ in $k=k_1,k_2,...$, $e^{ikR}\psi_{n}(r)$ corresponds to the same energy level as $\psi_{n}(r)$. But $e^{ikR}\psi_{n}(r)$ ($k=k_1,k_2,...$) is not the same as $\lbrace\psi_{nk}(r)\rbrace$, they are quite different. 
 A: I do not know where you found the description of  Bloch wavefunctions you presented, but your picture is misleading in my opinion. Referring to the standard notation $E_{n}(\vec{k})$,  $n$ is the band index and just indicates the band the energy level belongs to. 
To understand the rationale behind the picture you can initially omit $n$. In that case $\vec{k}$ completely fixes a simultaneous eigenfunction $\psi_k$ of $H$ and all $T_{\vec{R}}$. 
We have in particular
$$(T_{\vec{R}}\psi_{\vec k})(\vec{x})= e^{-i \vec{R}\cdot \vec{k}}\psi_{\vec k}(\vec{x})$$
where you see that the eigenvalue $e^{-i \vec{R}\cdot \vec{k}}$ of $T_{\vec{R}}$ is a function of $\vec{k}$.
The energy eigenvalue itself $E(\vec{k})$ is function of $\vec{k}$, too. 
The degeneracy of energy levels is due to the fact that 
just in view of the symmetry of the Hamiltonian you may have $$E(\vec k)= E(\vec h)\quad \mbox{for some $\vec k\neq \vec h$}\:.$$
The set of eigenvectors $\psi_{\vec k}$ such that $E(\vec{k}) = E$ span a closed subspace which is invariant under all $T_{\vec{R}}$ and therefore supports a unitary representation of the translational group of the $T_{\vec {R}}$. 
In more physical models, where $\vec{k}$ varies almost continuously, you see that the energies are grouped into bands. It happens that fixing  $\vec k$ is not enough to determine a unique simultaneous eigenstate of $H$ and all $T_{\vec{R}}$. To do it a further index $n$ is necessary and a Hilbert basis of simultaneous eigenfunctions is indicated by $\{\psi_{n \vec{k}}\}$.  This further index does not affect the eigenvalues of the translation operators,
$$(T_{\vec{R}}\psi_{n\vec k})(\vec{x})= e^{-i \vec{R}\cdot \vec{k}}\psi_{n \vec k}(\vec{x})$$
but it affects the energy eigenvalues
$$H \psi_{n \vec k} = E_n(\vec k) \psi_{n \vec k}$$
 and generally you have
$$E_{n}(\vec{k}) \neq E_{n'}(\vec{k}) \quad \mbox{if} \quad n \neq n'$$
Moreover for fixed $n$, the map $\vec{k} \mapsto E_n(\vec{k})$ is approximatively continuous and defines a so-called band of energy.
Two different bands are separated by a gap $\Delta E(\vec{k}) = E_{n}(\vec{k}) - E_{n'}(\vec{k})$ which may depend on $\vec{k}$.
