I'm frustrated that I can't understand something that must be simple and fundamental. I'd appreciate any answer to the question, but also any clarifications of how my presentation of the theorem/proof, or general quantum theory, is incorrect. Apologies for the long question; I could not figure out how to present the question completely without all of these details.
Bloch's Theorem states that, for a system with a periodic Hamiltonian — $H({\bf r} + {\bf R}) = H({\bf r})$ for all ${\bf R} \in A$, where $A$ is the set of all lattice symmetry translations, i.e., all integer linear combinations of the primitive lattice vectors $({\bf a_1}, {\bf a_2}, {\bf a_3})$ — any energy eigenstate satisfying $H | \psi \rangle = E | \psi \rangle$ "can be chosen so that associated with each $\psi$ is a wavevector ${\bf k}$ such that$^1$:"
$$
\psi({\bf r} + {\bf R}) = {\rm e}^{i {\bf k} \cdot {\bf R}} \, \psi( {\bf r})
$$
for each ${\bf R}\in A$. If we write:
$$
\psi({\bf r}) = {\rm e}^{i {\bf k} \cdot {\bf r}} u_{ {\bf k}}({\bf r})\,,
$$
then we see that
$$
\psi({\bf r} + {\bf R}) = {\rm e}^{i {\bf k} \cdot {\bf r}} {\rm e}^{i {\bf k} \cdot {\bf R}} \, u_{ {\bf k}}( {\bf r} + {\bf R})\, ,
$$
which implies that $u_{\bf k}({\bf r} + {\bf R}) = u_{\bf k}({\bf r})$ for all ${\bf R} \in A$. This is an alternative statement of the theorem: any energy eigenstate can be written as a plane wave times a function, labeled by the wavevector ${\bf k}$, that has the periodicity of the lattice.
Proofs of the theorem, e.g., from Ashcroft and Mermin or the wikipedia page, rely on choosing the wavevector to be of the form:
$$
{\bf k} = {\bf x} \cdot {\bf b} = x_1 {\bf b_1} + x_2 {\bf b_2} + x_3 {\bf b_3} \quad ({\bf x} \in \mathbb{R}^3)
$$
where the ${\bf b_i}$ are the primitive reciprocal lattice vectors satisfying ${\bf a_i} \cdot {\bf b_j} = 2 \pi \delta_{ij}$. The fact that the Hamiltonian commutes with each unitary translation operator $T_{\bf R}$ implies that energy eigenstates are also eigenstates of these translation operators. Letting
$$
{\bf R} = {\bf n} \cdot {\bf a} = n_1 {\bf a_1} + n_2 {\bf a_2} + n_3 {\bf a_3} \quad ({\bf n} \in \mathbb{Z}^3)
$$
we can write the eigenvalues of the translation operators as
$$
T_{\bf R} \, \psi({\bf r}) = \exp[i\, 2\pi \, (n_1 \tilde{x}_1 + n_2 \tilde{x}_2 + n_3 \tilde{x}_3)] \, \psi({\bf r})
$$
where the $\tilde{x}_i$ values specify the eigenvalues of the ${\bf a_i}$ primitive translations, i.e.,
$$ T_{\bf a_i} \, \psi({\bf r}) = \exp[i\, 2 \pi \, \tilde{x}_i] \, \psi({\bf r})\,.
$$
But, if we choose ${\bf k} = {\bf \tilde{x}} \cdot {\bf b}$, then we can write
$$
T_{\bf R} \, \psi({\bf r}) = \exp[i\, {\bf k} \cdot {\bf R}] \, \psi({\bf r})
$$
for all ${\bf R}\in A$. Therefore:
$$
\psi({\bf r} + {\bf R}) = T_{\bf R} \, \psi({\bf r}) = {\rm e}^{i\,{\bf k} \cdot {\bf R}} \, \psi({\bf r})
$$
for all ${\bf R} \in A$, showing that we can find a ${\bf k}$ such that the theorem is satisfied.
My confusion here is that it seems that the vector ${\bf k}$ is uniquely specified given an energy eigenstate $\psi$, i.e., it is given by the values $(\tilde{x}_1, \tilde{x}_2, \tilde{x}_3)$ that yield the eigenvalues of the primitive translations, ${\bf a_i}$, acting on that eigenstate $\psi$. That interpretation of the proof must be incorrect, however, because it implies a single ${\bf k}$ vector for each energy eigenstate. And, it is my undestanding, from the application of Bloch theory, is that the wavevector ${\bf k}$ should be thought of as a free parameter, an arbitrary real-valued 3-vector label (though its values can be restricted to the "first Brillion zone"). Indeed it is this continous freedom that gives the bands of energies. What am I missing here?
Possible solution: If I write the simultaneous energy and translation eigenstates as $| \psi_{\alpha, {\bf k}}\rangle$ then I must write the action of the translation operator as:
$$
T_{\bf R} \, |\psi_{\alpha, {\bf k}} \rangle = \exp[i\, 2\pi \, (n_1 \tilde{x}_{\alpha, {\bf k},1} + n_2 \tilde{x}_{\alpha, {\bf k},2} + n_3 \tilde{x}_{\alpha, {\bf k},3})] \, | \psi_{\alpha, {\bf k}} \rangle
$$
and then I choose ${\bf k} = {\bf \tilde{x}_{\alpha, {\bf k}}} \cdot {\bf b}$. This means that the $\tilde{x}_i$ values are specific to an eigenstate labeled by $(\alpha, {\bf k})$, not just by $\alpha$. Therefore we actually have complete freedom to choose ${\bf k} \in \mathbb{R}^3$. This argument feels a bit circular to me.
1 – From Ashcroft and Mermin, Chapter 8.
Best Answer
In what follows, operators have hats, un-hatted things are real or complex numbers
Suppose (in full generality, not necessarily assuming that we have a single particle state) that we have a Hamiltonian operator $\hat H$ that is periodic in the sense that $[\hat T_{\mathbf{R}}, \hat H] = 0$ for all $\mathbf{R} = n_1 {\bf a}_1 + n_2 {\bf a}_2 + n_3 {\bf a}_3$, where $\hat T_{\mathbf{R}} = \exp(i {\bf \hat p \cdot R })$ is the translation operator and $\underline{n} = (n_1,n_2,n_3) \in \mathbb{Z}^3$.
The game we are playing here is simply Hamiltonian diagonalisation. The facts we are using are
Using results (1-3), we find that there exists a spanning set for the Hilbert space $\{ | E, \lambda(\underline{n}) \rangle \}$, defined by $$\hat H | E, \lambda(\underline{n}) \rangle = E | E, \lambda(\underline{n}) \rangle\\ \hat T_{\underline{n}} | E, \lambda(\underline{n}) \rangle = \lambda(\underline{n}) | E, \lambda(\underline{n}) \rangle$$
I called this a spanning set because it is not yet clear whether these vectors are linearly independent. In fact, they are not: result (4) implies that $\lambda(\underline{n}) = \lambda(1,0,0)^{n_1} \lambda(0,1,0)^{n_2} \lambda(0,0,1)^{n_3}$, i.e. we can throw away all but three $\hat{T}_{\underline{n}}$ eigenvectors and span exactly the same space. In other words, we gain no additional information about a state from translations by more than one unit cell.
In this sense the only independent, commuting operators are
$$\hat H, \hat T_{\mathbf{a_1}}, \hat T_{\mathbf{a_2}}, \hat T_{\mathbf{a_3}}.$$
Unitarity allows us to rename $\lambda(1,0,0) = e^{ik_1'}, \lambda(0,1,0) = e^{ik_2'}, \lambda(0,0,1) = e^{ik_3'}$, or more concisely
$$ \lambda(\underline n) = \exp(i\underline k \cdot \underline n) $$
for $\underline{k}'$ chosen from $(\mathbb{R}/0\sim 2\pi)^3$. Reintroducing the units (via $\sum_i k_i' n_i = k_i \mathbf{b}_i \cdot \mathbf{R}_{\underline n} \equiv \mathbf{k} \cdot \mathbf{R}_{\underline n}$) proves the statement,
The allowed values of $k$
Recall that the defining eigenvalues of a particular quantum state are
$$ E, e^{ik'_x},e^{ik'_y},e^{ik'_z} $$
It remains to see how the $k'_{x,y,z}$ are arranged. Firstly, note that the physical thing is $e^{ik_j}$ - because the complex exponential is periodic, the natural way to think of the $k$ values is as real numbers modulo $2\pi$ - $k$ and $k + 2\pi$ label precisely the same state.
Suppose that instead of being infinite, you put your lattice on a torus such that there are finitely many ($L$, say) unit cells in every direction. The translation operator becomes compactified due to $\hat{T}_{L\mathbf{a_j}} = \hat{1} \forall i=1,2,3$, and further, this implies that $\lambda(1,0,0)^L = 1$, i.e. $k'_i = \frac{2\pi n}{L}, n \in \mathbb{Z}$. Note that in 1D, this gives you exactly $L$ distinct eigenvalues modulo $2\pi$, and so in 3D, you end up with $L^3$ possible $\lambda(\underline n)$ eigenvalues / $k_{\underline n}$ (projective) eigenvalues. This argument established that the spectrum of the operator $\hat{T}_{\mathbf{a}_1}$ consists of $L$ evenly spaced complex phases on the unit circle.
Note the subtlety here - there are $L^3 + 1$ commuting operators ($\hat{T}_{\underline{n}}, \hat{H}$). If there is one local degree of freedom per unit cell (for a total of $N$ degrees of freedom), the system is actually overconstrained by the $k$ eigenvalues - this effectively makes $E$ a function of $k$. If there are two degrees of freedom per unit cell (e.g. graphene), then there are two bands, and $E$ is needed to distinguish them.
This counting argument suggests that it is helpful to express the full Hilbert space $\mathcal{H}$ as the product space $\otimes_{k \in BZ} \mathcal{H}_k$.
where the dimension of $\mathcal{H}_k$ is independent of the lattice size. As one takes $L \to \infty$, you find that $k$ can be any rational multiple of $2\pi$. Since the rationals are dense, one can argue for replacing the sum by an integral, but these technicalities are not really needed for solving most problems.