Bloch's theorem tells us that the energy eigenvectors of a Hamiltonian with a periodic potential can be written
$$\psi_{n\mathbf k}(\mathbf x) = e^{i\mathbf k \cdot \mathbf x} u_{n\mathbf k}(\mathbf x)$$
where $n\in \mathbb Z$, $\mathbf k\in \mathrm{BZ}$ (the first Brillouin zone), and $u_{n\mathbf k}(\mathbf x)$ is periodic with the same periodicity as the lattice. Applying the Hamiltonian operator yields
$$\big(H \psi_{n\mathbf k}\big)(\mathbf x)= \left[-\frac{\hbar^2}{2m} \nabla^2 + V(\mathbf x)\right]e^{i\mathbf k \cdot \mathbf x}u_{n\mathbf k}(\mathbf x) $$
$$= e^{i\mathbf k\cdot \mathbf x}\left[-\frac{\hbar^2}{2m}(\nabla + i\mathbf k)^2 + V(\mathbf x) \right]u_{n\mathbf k}(\mathbf x) = E_{n\mathbf k} e^{i\mathbf k\cdot \mathbf x} u_{n\mathbf k}(\mathbf x) \tag{$\star$}$$
Cancelling the factor $e^{i\mathbf k \cdot \mathbf x}$ from both terms in $(\star)$ yields that $u_{n\mathbf k}$ is a solution of the equation $H_\mathbf k u_{n\mathbf k} = E_{n\mathbf k} u_{n\mathbf k}$, where $H_{\mathbf k} \equiv -\frac{\hbar^2}{2m} (\nabla +i\mathbf k)^2 + V(\mathbf x)$, defined on the unit cell with periodic boundary conditions.
More concretely, let $\mathscr u$ denote the unit cell. Consider the Hilbert space $L^2(\mathscr u)$ of square-integrable functions on the unit cell equipped with the standard inner product
$$\langle \psi,\phi\rangle := \int_{\mathscr u} \mathrm d^n x \ \overline{\psi(\mathbf x)} \phi(\mathbf x)$$
Further define the Bloch Hamiltonian $H_\mathbf k$ to act on the twice-weakly differentiable elements of $L^2(\mathscr u)$ with periodic boundary conditions. One can show that $H_{\mathbf k}$ is self-adjoint with discrete spectrum, and therefore that one can construct an orthonormal basis $\{u_{n\mathbf k}\}$ of solutions to the eigenvalue equation $H_\mathbf k u_{n\mathbf k} = E_{n\mathbf k}u_{n\mathbf k}$.
From some comments,
Why do we go and change the relevant Hilbert space here? Put differently, I think that during this procedure, in eq. (2), $H_k$ is an operator of $L^2(V)$ (sq. int. function on $V$, the sample volume, respecting the BvK b.c.) and moreover $u_{nk}\in L^2(V)$. Can't we say directly that $H_k$ is self-adjoint on this Hilbert space (but ofc. defined only on a dense subspace) and that then we can choose the $u_{nk}$ as orthonormal (or orthogonal and normalize it on the unit cell)? Why do you (and all textbooks I've seen) go first to $L^2(u)$ for this?
[...]
Why can't we leave the BvK boundary conditions? Through the calculation you've done, we simply find that $u_{nk}$ are eigenfunctions of $H_k$, which are lattice periodic due to the Bloch theorem anyway and hence also fulfill the BvK boundary conditions. [From another point of view: One can see $V$ as a single unit cell of a crystal with volume $NV$ and the math would be the same as you did, namely treating $H_k$ as an operator over the unit cell $V$ of this new crystal.] I don't see a reason why we would need to change boundary conditions to be lattice periodic.
It's worth taking a moment to think about what Bloch's theorem actually says. We begin with the Hamiltonian
$$H = -\frac{\hbar^2}{2m} \nabla^2 + V$$
where $V$ is a periodic function. We take the Hilbert space to be $L^2(V)$, where $V$ is the sample volume, and impose periodic boundary conditions on $V$. These are the so-called Born-von Karman boundary conditions.
The fact that $V$ is periodic means that $H$ commutes with a set of translation operators $T_i$. As per the usual recipe, since the $T_i$'s commute among themselves, we seek to simplify the problem by finding a simultaneous eigenbasis of $H$ and the $T_i$'s.
One can show without excessive difficulty that if a function $\psi$ is an eigenfunction of the $T_i$'s, then it must be of the form $\psi(x) = e^{i\mathbf k \cdot \mathbf x}u(\mathbf x)$ where $u$ is periodic with the same periodicity as the lattice. This identification is not unique, however. If $\mathbf b$ is any arbitrary inverse lattice vector, then the function $\exp[i\mathbf b \cdot \mathbf x]$ is also periodic on the lattice, and so
$$e^{i\mathbf k\cdot \mathbf x}u(\mathbf x) = e^{i\mathbf k \cdot \mathbf x} e^{i \mathbf b \cdot \mathbf x} e^{-i\mathbf b \cdot \mathbf x} u(\mathbf x) = e^{i\tilde{\mathbf k} \cdot \mathbf x} \tilde u(\mathbf x)$$
where $\tilde{\mathbf k}= \mathbf k + \mathbf b$. In order to resolve this ambiguity, we restrict $\mathbf k$ to the first Brillouin zone, which makes the decomposition of $\psi$ into $e^{i\mathbf k \cdot \mathbf x}u(\mathbf x)$ unique.
To recap, a simultaneous eigenfunction of all of the $T_i$'s can always be uniquely expressed in the form $e^{i\mathbf k \cdot \mathbf x} u(\mathbf x)$, where $\mathbf k$ lies in the first Brillouin zone and $u$ is periodic with the same periodicity as the lattice. Until this point, we've said nothing of the Hamiltonian, but we can plug such a function in to the eigenvalue equation to obtain
$$H \psi = E \psi$$
$$\underbrace{\left[-\frac{\hbar^2}{2m}\left(\nabla + i\mathbf k\right)^2 + V(\mathbf x)\right]}_{\equiv H_\mathbf k} u(\mathbf x) = E u(\mathbf x)$$
where we have done some algebra to cancel out the $e^{i\mathbf k \cdot \mathbf x}$ from each side.
As a result, we find that a simultaneous eigenvector of the $T_i$'s and the Hamiltonian takes the form $e^{i\mathbf k \cdot \mathbf x} u(\mathbf x)$ where $\mathbf k$ lies in the first Brillouin zone, $u$ is periodic with the same periodicity as the lattice, and $u$ is an eigenfunction of $H_\mathbf k$ as defined above. This is the Bloch theorem.
The point of all this is that $H$ has a complicated, continuous spectrum, which is difficult to deal with. Bloch's theorem allows us to choose a $\mathbf k$ from the first Brillouin zone and then seek the lattice-periodic eigenfunctions $u_{n\mathbf k}$ of $H_\mathbf k$, which form a nice discrete set; from there, the corresponding eigenfunction of $H$ is $e^{i\mathbf k\cdot \mathbf x}u_{n\mathbf k}(\mathbf x)$.
