How to justify unit cell normalization of Bloch functions? As is widely known, non-square-integrable wavefunctions don't belong to the Hilbert space, and therefore cannot represent physical states.
This is the case for e.g. oscillating wavefunctions and generally wavefunctions that do not vanish sufficiently rapidly at infinity (see however Normalizable wavefunction that does not vanish at infinity for counterexamples).
A particular example of such functions are Bloch waves, arising in solid state physics whenever we have a crystal lattice (and therefore a periodic potential). These are wavefunctions of the form (in one dimension)
$$\psi_k(x) = u_k(x) e^{ikx}$$
where $u_k(x)$ is a periodic complex-valued function with the same periodicity as the lattice. Such wavefunctions are clearly non-normalizable over the whole space since
$$\int_{- \infty}^{+ \infty} |\psi_k(x)|^2 dx = \int_{- \infty}^{+ \infty} |u_k(x)|^2 dx $$
and $|u_k(x)|^2$ is periodical.
Since we want to be able to calculate expectation values, a typical approach in solid state physics is to normalize this kind of functions not over the whole space, but only over the unit cell of the lattice. This in one dimension corresponds to writing the normalization condition as
$$ \int_{0}^{a} |\psi_k(x)|^2 dx = 1 $$
where $a$ is the cell parameter.
This surely allows to normalize the wavefunction, and to carry on calculations, but isn't equivalent to saying that we are sure to find the electron inside a particular unit cell, at every time?
How is such a procedure justified from a physical point of view?
 A: 
This surely allows to normalize the wavefunction, and to carry on calculations, but isn't equivalent to saying that we are sure to find the electron inside a particular unit cell, at every time?

This doesn't normalize the wavefunction; $\int_{-\infty}^\infty |\psi_{nk}(x)|^2 \mathrm dx$ is infinite, as always.  However, $\int_0^a |u_{nk}(x)|^2\mathrm dx$ is always going to be finite, and sometimes it happens to be convenient to re-scale $u_{nk}$ to make that quantity equal to $1$.  Alternatively, sometimes it is helpful to choose the normalization $\int_0^a |u_{nk}(x)|^2 \mathrm dx = a$.
For example, the "orthogonality" of differing eigenstates suggests that $\int \psi^*_{nk}(x) \psi_{mq}(x) \mathrm dx \propto \delta(k-q)\delta_{nm}$.  Indeed, if $u_{nk}$ is periodic with period $a$, then it can be expanded as a Fourier series:
$$u_{nk}(x) = \frac{1}{\sqrt{a}}\sum_{r\in \mathbb Z} c^{(nk)}_r e^{-2\pi i rx/a} \iff c^{(nk)}_r = \frac{1}{\sqrt{a}}\int_0^a u_{nk} (x) e^{2\pi irx/a} \mathrm dx$$
Therefore,
$$\psi_{nk}(x) = e^{ikx} u_{nk}(x) = \frac{1}{\sqrt{a}}\sum_{r\in \mathbb Z}c^{(nk)}_r e^{i(k-2\pi r/a)x}$$
As a result,
$$\int \mathrm dx \ \psi^*_{nk}(x) \psi_{mq}(x) = \frac{1}{a}\sum_{r,s\in \mathbb Z} c_r^{(nk)*} c_{s}^{(mq)} \int \mathrm dx\ e^{-i[k-q + 2\pi (r-s)/a]x}$$
$$= \frac{2\pi}{a} \sum_{r,s\in \mathbb Z} c_r^{(nk)*} c_{s}^{(mq)} \delta\big(k-q + \frac{2\pi}{a} (r-s)\big)$$
However, $k$ and $q$ are both taken from the first Brillouin zone, and so they differ by less than $2\pi/a$.  As a result, if $r-s\neq 0$ then the delta function will never "fire" and the result is zero.  Therefore, we can replace it with $\delta(k-q) \delta_{rs}$.  From there, we obtain
$$\int \mathrm dx \psi^*_{nk}(x) \psi_{mq}(x)  = \frac{2\pi}{a} \sum_{r\in \mathbb Z} c^{(nk)*}_r c^{(mk)}_r \delta(k-q)$$
We now notice that
$$\int_0^a u_{nk}^*(x) u_{mk}(x) \mathrm dx = \frac{1}{a} \sum_{r,s\in \mathbb Z} c_r^{(nk)*}c_s^{(mk)} \int_0^a \mathrm dx \ e^{2\pi i(r-s)x/a}$$
$$= \frac{1}{a} \sum_{r,s\in \mathbb Z} c_r^{(nk)*}c_s^{(mk)} a \delta_{rs} = \sum_{r\in \mathbb Z}c^{(nk)*}_rc^{(mk)}_r$$
Defining the Bloch function inner product
$$\langle u_{nk}, u_{mk}\rangle \equiv \int_0^a u_{nk}^*(x) u_{mk}(x) \mathrm dx$$
yields
$$\int \psi^*_{nk}(x)\psi_{mk}(x) \mathrm dx = \frac{2\pi}{a} \langle u_{nk},u_{mk}\rangle\delta(k-q)$$
Finally, we note that $u_{nk}$ and $u_{mk}$ are eigenfunctions of the Bloch Hamiltonian $H(k)$ with eigenvalues $E_n(k)$ and $E_{m}(k)$ respectively; since $H(k)$ is self-adjoint, it follows that we may choose $u_{nk}$ and $u_{mk}$ to be orthogonal (of course, if $E_n(k) \neq E_m(k)$ then they are automatically orthogonal).  Choosing them to be normalized such that $\langle u_{nk},u_{mk}\rangle = a\delta_{nm}$ yields the final result
$$\int_{-\infty}^\infty \psi^*_{nk}(x) \psi_{mq}(x) \mathrm dx = 2\pi  \delta_{nm} \delta(k-q)$$
as expected by analogy with the plane waves of a free particle.  On the other hand, choosing the OP's normalization $\langle u_{nk},u_{mk}\rangle = \delta_{nm}$ yields an extra factor of $a$ in the above; however, if you only care about the $u_{nk}$'s (which is common when examining the topology of the band structure, for example) then this may be more convenient.
A: I think that normalizing on a single cell does not mean that

we are sure to find the electron inside a particular unit cell

I think is just a way of simplifying calculation exploiting the fact that for certain periodic potential it doesn't matter on how many cells our wave function is defined on, any domain you choose yield the same expectation value.
Here I give a plausibility argument of this.
Think of the Bloch wave defined on the whole space, of the shape $\psi_k(r)=\exp(ikr)u_k(r)$, with u_k periodic with periodicity $a$, the same of the lattice.
Now slice it, so that you have one piece per cell.
Take one cell, normalize the piece of function to that cell:
$$\begin{eqnarray} 
1 &=& \int_{0}^{a} \psi_k(r)^{*}\psi_k(r) \\
&=& \int_{0}^{a} \exp(-ikr)u_k(r)^{*}\exp(ikr)u_k(r) \\
&=&\int_{0}^{a} u_k(r)^{*}u_k(r) 
\end{eqnarray}$$
Now let V a periodic potential with the same periodicity of the lattice and that
$$V(r)\exp(\pm ikr)u_k(r)=\exp(\pm ikr)V(r) u_k(r)$$
The expectation value of V over the single cell is:
$$\begin{eqnarray} 
\int_{0}^{a} \psi_k(r)^{*}V(r)\psi_k(r)&=&
\int_{0}^{a} \exp(-ikr)u_k(r)^{*}V(r)\exp(ikr)u_k(r) \\
&=&\int_{0}^{a} u_k(r)^{*}V(r)u_k(r)
\end{eqnarray}$$
Now take, say 5, different cells at random, and consider the wave function on this domain, we normalize it so that
$$ \sum_{i=1}^{5} \int_{x_i}^{x_i+a} \psi^{i}_{k}(r)^{*}\psi^{i}_{k}(r) = 1 $$
where $x_i$ are lattice points.
Due to the periodicity we are summing 5 copies of the same value so that each integral is:
$$\begin{eqnarray} 
\int_{x_i}^{x_i+a} \psi^{i}_{k}(r)^{*}\psi^{i}_{k}(r)(r) &=& \frac{1}{5} \\ &=&\int_{x_i}^{x_i+a} u^{i}_k(r)^{*}u^{i}_k(r)
\end{eqnarray}$$
then $u_k(r+x_i)=\sqrt(5)*u^{i}_k(r)$ (since $u_k$ and $u^{i}_k$ are slices of the same function (which is periodic) apart from a normalization constant).
when we calculate the expected value of V we get:
$$\begin{eqnarray} 
 \sum_{i=1}^{5} \int_{x_i}^{x_i+a} \psi^{i}_{k}(r)^{*}V(r)\psi^{i}_{k}(r) &=&\\
&=& \sum_{i=1}^{5} \int_{x_i}^{x_i+a} \exp(-ikr)u^{i}_k(r)^{*}V(r)\exp(ikr)u^{i}_k(r) \\
&=&   5*\int_{x_i}^{x_i+a} u^{i}_k(r)^{*}V(r)u^{i}_k(r) \\
&=& \int_{x_i}^{x_i+a} u_k(r+x_i)^{*}V(r)u_k(r+x_i)dr   \\
&=& \int_{0}^{a} u_k(y)^{*}V(r)u_k(y)dy
\end{eqnarray}$$
Which is the same expectation value we had obtained when normalizing on a single cell.
