General Question(s) on Bloch's Theorem (Apologies for the long set-up, I wanted to make sure the question was properly motivated. Skip to the bottom for the actual question.)
Bloch's theorem is a foundational theorem in solid-state physics. It states, given a single-particle Hamiltonian $H = T + V$ with a periodic potential,
$$
V(\vec{r}+\vec{R}) = V(\vec{r}),
$$
the eigenstates of $H$ can always be chosen to satisfy
$$
\psi_{n\vec{k}}(\vec{r} + \vec{R}) = e^{i \vec{k} \cdot \vec{R}} \psi_{n\vec{k}}(\vec{r}) .
$$
Here, $\vec{R}$ refers to a vector in the Bravais lattice, which defines the periodicity of $V$. The theorem is easy to prove: it essentially follows from choosing eigenstates which are simultaneous eigenstates of the discrete translation operators. Alternatively, it follows by simply rewriting the Hamiltonian in Fourier space, where it is immediately seen that $V$ connects only states differing in momentum by a reciprocal lattice vector $\vec{K}$; therefore $H$ can be block diagonalized into sectors labeled by momenta $\vec{k}$ in the first Brillouin zone, each $\vec{k}$th sector containing all momentum eigenstates differing from $\vec{k}$ by a reciprocal lattice vector.
There are a number of statements which commonly follow the proof of Bloch's theorem (for example, see the discussion in Ashcroft and Mermin). First, it is claimed that the additional index $n$ used to resolve any degeneracies in the crystal momentum $\vec{k}$ is necessarily discrete. This is argued by considering the differential equation obeyed by the cell-periodic Bloch function $u_{n\vec{k}}(\vec{r}) = e^{-i\vec{k} \cdot \vec{r}} \psi_{n\vec{k}}(\vec{r})$:
$$
\left[ \frac{1}{2m} \left( -i \hbar \vec{\nabla} + \hbar \vec{k} \right)^2 + V(\vec{r}) \right] u_{n\vec{k}}(\vec{r}) = \varepsilon_{n\vec{k}} u_{n\vec{k}}(\vec{r})
$$
together with the boundary condition $u_{n\vec{k}}(\vec{r} + \vec{R}) = u_{n\vec{k}}(\vec{r})$. It is then claimed that since $u_{n\vec{k}}(\vec{r})$ satisfies this eigenvalue problem in a single unit cell, the spectrum of solutions for fixed $\vec{k}$ must be discrete, and the additional index $n$ used to label this spectrum can be chosen as a discrete index. The second statement commonly made is that for given $n$, the Bloch wavefunctions can always be chosen to be periodic in the Brillouin zone: $u_{n,\vec{k} + \vec{K}} = u_{n\vec{k}}$.
Both of these points are central to the theory of band structure. Together with the assumption that the eigenvalues $\varepsilon_{n\vec{k}}$ vary continuously with $\vec{k}$, the first point states that the spectrum of $H$ is organized into bands. The second point states that these bands are always periodic in the Brillouin zone.
Here are the questions, which I hope are sufficiently closely related to comprise one overall question:

*

*How can I see more precisely that the spectrum in $n$ is discrete for a fixed crystal momentum $\vec{k}$? In other words, given the above differential equation for the cell-periodic Bloch function, is there a way to show that the eigenvalue spectrum of the differential equation is discrete?

*While it's perfectly clear that I can take the crystal momenta to lie within the Brillouin zone, how can I see that each band must strictly be periodic in the Brillouin zone? For example, consider a one-dimensional problem for ease of visualizing the bands. Why can't the $n$th band start at some fixed energy $\varepsilon_{n, -\pi / a}$, increase monotonically with $k$, and end at $\varepsilon_{n,\pi / a} = \varepsilon_{n+1, -\pi / a}$, so that each band continuously wraps into the next one?

 A:  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

*

*Kottke (2016), "Linear Analysis on Manifolds: Notes for Math 7376, Spring 2016" (https://ckottke.ncf.edu/neu/7376_sp16/lacm.pdf).

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.
