This is not true--- you can make superpositions of degenerate energy eigenstates and the are still eigenstates. For example, moving at wavenumber k and -k have the same energy, but the difference is not in the reciprocal lattice.
The proper thing to say is that the Hamiltonian is translationally invariant with translation operator the lattice translations, and this means that H commutes with lattice translations, so that the eigenvalues of the lattice translations (the equivalence class of k under the reciprocal lattice) are each separately sectors where you find eigenvectors of H.
Question 1
The functions
$
\newcommand{\bfr}{\mathbf{r}}
\newcommand{\bfR}{\mathbf{R}}
\newcommand{\bfk}{\mathbf{k}}
\newcommand{\bfK}{\mathbf{K}}
u_{n,\bfk}
$ are solutions to the third equation in the question. Everything in that equation, including $V$ and $u$, is assumed to be invariant under $\bfr\to \bfr+\bfR$ for lattice vectors $\bfR$. Therefore, for each value of $\bfk$, that equation defines an eigenvalue problem on a compact manifold (specifically a torus) consisting of points $\bfr$ modulo the vectors $\bfR$. Here's the key: for any given value of $\bfk$, the differential operator in that equation is self-adjoint and elliptic, and
- The eigenvalues$^\dagger$ of a self-adjoint elliptic differential operator on a compact manifold are discrete.
$^\dagger$ In physics, the word "eigenvalue" is often used loosely, referring also to elements of the continuous part of an operator's spectrum. That's how I'm using the word here, so the statement isn't just a tautology.
This is one of those general theorems in the phycisist's standard toolkit, with many applications. For a more precise statement and an outline of the proof, see theorem 1.46 on page 47 in
The proof is not short, so I won't try to reproduce it here.
Question 2
The condition $u_{n,\bfk+\bfK}=u_{n,\bfk}$ shown in the question might be a typo (?). The correct statement is that we can take
$$
\psi_{n,\bfk+\bfK}(\bfr)=\psi_{n,\bfk}(\bfr).
\tag{1}
$$
This is equation 8.50 in Ashcroft and Mermin's Solid State Physics, 1976. To see why equation (1) is valid, let
$$
\psi_{n,\bfk}(\bfr)=\exp(i\bfk\cdot\bfr)u_{n,\bfk}(\bfr)
\tag{2}
$$
be any solution of the original Schrödinger equation (not shown in the question), where $u_{n,\bfk}(\bfr)$ is invariant under $\bfr\to\bfr+\bfR$. We can obviously rewrite (2) as
$$
\psi_{n,\bfk}(\bfr)=\exp\big(i(\bfk+\bfK)\cdot\bfr\big)
\Big(\exp(-i\bfK\cdot\bfr)u_{n,\bfk}(\bfr)\Big).
\tag{3}
$$
The combination $\exp(-i\bfK\cdot\bfr)u_{n,\bfk}(\bfr)$ is still invariant under $\bfr\to \bfr+\bfR$, so given any complete basis of solutions (2) for each $\bfk$ in the first Brillouin zone, we can simply define the Bloch wavefunctions for the other Brillouin zones by
$$
\psi_{n,\bfk+\bfK}(\bfr)\equiv
\exp\big(i(\bfk+\bfK)\cdot\bfr\big)
u_{n,\bfk+\bfK}(\bfr)
\tag{4}
$$
for all $\bfk$ in the first Brillouin zone and for all reciporocal lattice vectors $\bfK$, with
$$
u_{n,\bfk+\bfK}(\bfr)
\equiv\exp(-i\bfK\cdot\bfr)u_{n,\bfk}(\bfr).
\tag{5}
$$
Equations (3)-(5) imply equation (1). To be consistent with how we defined $\psi_{n,\bfk+\bfK}$, we can also define
$$
\epsilon_{n,\bfk+\bfK}\equiv\epsilon_{n,\bfk}
\tag{6}
$$
for all $\bfk$ in the first Brillouin zone and for all reciporocal lattice vectors $\bfK$. We can take (1) and (6) to depend smoothly on $\bfk$ within the first Brillouin zone, but this does not imply that $\epsilon_{n,\bfk}$ is a smooth function of $\bfk$ when $\bfk$ crosses a zone boundary! Here's the key point:
- The choices (1) and (6) enforce periodicity in $\bfk$ instead of enforcing smoothness in $\bfk$ at the zone boundaries.
We can think of equations (1) and (6) as defining how the wavefunctions and bands are indexed in the other Brillouin zones, given a choice of indexing in the first Brillouin zone. With this definition, the index $n$ may jump when tracing a band smoothly (differentiably) across a zone boundary, as anticipated in the question. Alternatively, we could have used a different definition in which they are taken to vary smoothly with $\bfk$ across all zones, for each $n$, but that way of indexing things might not agree with the one defined by (1) and (6). Both ways of indexing things are legitimate, but they're not necessarily equivalent.
To see that they're not always equivalent, consider the case $V=0$. This is periodic with respect to any lattice, so choose an arbitrary lattice and consider equations (1) and (6). If we had indexed the bands to be smooth for all $\bfk$, then $\epsilon_{n,\bfk}$ would be a quadratic function of $\bfk$ that increases without bound as $|\bfk|$ increases. That's not consistent with (6), so this is a case where enforcing periodicity is not the same as enforcing smoothness.
Best Answer
(1) Since $u(\textbf{r}) = u(\textbf{r}+\textbf{R})$, we can expand this part in terms of reciprocal lattice vectors, $u_k(\textbf{r}) = \sum_\textbf{G}{e^{i\textbf{G}\cdot \textbf{r}}u_\textbf{k-G}}$. We can therefore write:
\begin{equation} \psi_{\textbf k+\textbf K} = e^{i(\textbf k + \textbf K)\cdot \textbf r}\sum_\textbf{G'}{e^{i\textbf{G'}\cdot \textbf{r}}u_{\textbf k-\textbf K- \textbf G'}} = e^{i\textbf k \cdot \textbf r}\sum_\textbf{G'}{e^{i(\textbf{G'}+\textbf K)\cdot \textbf{r}}u_{\textbf k-\textbf K- \textbf G'}}=e^{i\textbf k \cdot \textbf r}\sum_\textbf{G}{e^{i\textbf{G}\cdot \textbf{r}}u_{\textbf k-\textbf G}} = \psi_\textbf k \end{equation} where $\textbf G = \textbf K+\textbf G'$.
(2) You can interpret $\textbf p'$ as being equal to $\textbf p$. This is true because the real space lattice is periodic; $\textbf k$ is always equal to $\textbf k + \textbf K$.
(3) The conserved quantity is $\textbf k$ $\textit mod$ $\textbf K$. You can see that I used this fact in the answer to (2).
You can read just about any solid state physics textbooks for complete justification though my personal favorite is Ziman's Theory of Solids.