Why are bloch factor orthogonal? The Bloch wave can be expressed as:
$$
\psi_{n\mathbf{k}}(\mathbf{r}) = u_{n\mathbf{k}}(\mathbf{r})\,e^{i\mathbf{k}\cdot \mathbf{r}} \tag{A1}
$$
In this problem Bloch wave they say that $u_{n\mathbf{k}}(r)$ is orthogonal. I would like to ask whether $u_{n\mathbf{k}}(r)$ itself can be non-orthogonal, but if the Bloch wave is a set of orthonormal basis, the premise is that $u_{n\mathbf{k}}(r)$ must be orthogonal, so we mandate:
$$
\int_{\mathrm{unit \,cell}} u_{n\mathbf{k}}(\mathbf{r})\,u_{m\mathbf{k}}(\mathbf{r})d\mathbf{r} = \delta_{nm} \tag{A2}
$$
$\delta$ is the Dirac Function.

Thanks to a commenter for the reminder that in both answers one and two they give the origin of the $u_{n\mathbf{k}}(\mathbf{r})$ quadrature and state that this is derived from such an equation:
$$
[\dfrac{(i\hbar\nabla + \hbar\mathbf{k})^2}{2m} + V(\mathbf{r})] u_{n\mathbf{k}}(\mathbf{r}) = E_{n\mathbf{k}}u_{n\mathbf{k}}(\mathbf{r}) \tag{B1}
$$
I know where this wave equation came from, First use the momentum operator $\hat{p}=-i\hbar \nabla$:
$$
\begin{align}
\hat{p}\psi_{n\mathbf{k}}(\mathbf{r}) =& e^{i\mathbf{k}\cdot \mathbf{r}}(\hat{p} + \hbar \mathbf{k})u_{n\mathbf{k}}(\mathbf{r}) \\
\hat{p}^2\psi_{n\mathbf{k}}(\mathbf{r}) =& e^{i\mathbf{k}\cdot \mathbf{r}}(\hat{p} + \hbar \mathbf{k})^2u_{n\mathbf{k}}(\mathbf{r})
\end{align}
$$
Substituting this into the Schrodinger equation gives eq(B1), but I don't know how to derive eq(A2) from eq(B1)
 A: 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}$.
