Dissatisfied with textbook explanations for why $\vec k$ in Bloch's theorem can be restricted to thefirst Brillouin Zone (BZ) By Bloch's theorem, all the eigenfunctions of a Hamiltonian with a periodic potential $$U({\vec r}+{\vec R})=U({\vec r})$$  can be chosen to have the form $$\psi_{n{\vec k}}({\vec r})=e^{i{\vec k}\cdot{\vec r}}u_{n{\vec k}}({\vec r})$$ where ${\vec k}$ can be any reciprocal lattice vector. However, the statement is that one can always confine the range of $\vec k$ to be inside the first Brillouin Zone (BZ). At this point, let me explain my confusion.
Suppose I first consider a wavevector $\vec k_{\rm BZ}\in{\rm BZ}$, for which Bloch theorem above reads $$\psi_{n\vec k_{\rm BZ}}(\vec r)=e^{i\vec k_{\rm BZ}\cdot\vec{r}}u_{n\vec k_{\rm BZ}}(\vec r).$$
Next, consider another wavevector $\vec k$ that lies outside the BZ. But any $\vec k\notin {\rm BZ}$ can always be written as $$\vec k=\vec k_{\rm BZ}+\vec G$$ where  $\vec k_{\rm BZ}\in {\rm BZ}$ and $\vec G$ is an appropriate reciprocal lattice translation vector.  
Substituting $\vec k=\vec k_{\rm BZ}+\vec G$ in the Bloch's theorem above, we find $$\psi_{n\vec k}(\vec r)=e^{i\vec{k}_{\rm BZ}\cdot\vec{r}}\underbrace{e^{i\vec G\cdot\vec r}u_{n\vec{k}}(\vec{r})}\equiv e^{i\vec{k}_{\rm BZ}\cdot\vec{r}}u^\prime_{n\vec k
}(\vec r).$$ In the last step, we noted that the function inside the underbrace is also periodic with the periodicity $\vec{R}$ on account of the fact $e^{i\vec G\cdot\vec R}=1$ and defined it to be $$u^\prime_{n\vec{k}}(\vec r)\equiv e^{i\vec G\cdot\vec r}u_{n\vec{k}}(\vec{r}).$$


*

*My impression was/is that the reason we can remain restricted to BZ because $\psi_{n\vec b_{\rm BZ}}$, in some sense, is not a different solution from $\psi_{n\vec k}$. But it is not clear to me how or in what way, $\psi_{n\vec k}=\psi_{n,\vec k_{\rm BZ}+\vec G}$ is not a different solution from $\psi_{n\vec b_{\rm BZ}}$. Please help.

 A: Consider the 1D Schrödinger's equation with $U(x+a)=U(x)$:
$$-f''(x)+U(x)f(x)-Ef(x)=0,\tag1$$
where $f$ is required (as a boundary condition) to be bounded at infinity: $|f(x)|<\infty$ as $|x|\to\infty$.
Substituting
$$f(x)=u_k(x)\exp(ikx),\tag2$$
we get
$$-u_k''(x)-2iku_k'(x)+(k^2-E+U(x))u_k(x)=0.\tag3$$
By Bloch's theorem, boundary conditions for $u_k$ are
$$
\left\{
\begin{align}
u_k(x+a)&=u_k(x),\\
u_k'(x+a)&=u_k'(x).
\end{align}
\right.\tag4
$$
Now substituting
$$k=k_1+G,\tag5$$
where $G$ is a multiple of $2\pi/a$, we transform $(3)$ into
$$-u_{k_1+G}''(x)-2i(k_1+G)u_{k_1+G}'(x)+((k_1+G)^2-E+U(x))u_{k_1+G}(x)=0.\tag6$$
The boundary conditions for this equation are unchanged, since it's just a replacement of a parameter.
But if now we substitute
$$u_{k_1+G}(x)={u_1}_{k_1}(x)\exp(-iGx),\tag7$$
we'll get an equation
$$-{u_1}_{k_1}''(x)-2ik_1{u_1}_{k_1}'(x)+(k_1^2-E+U(x)){u_1}_{k_1}(x)=0,\tag8$$
which is isomorphic to $(3)$. Moreover, since $\exp(-iGx)$ is periodic with period of $2\pi/a$, boundary conditions are also the same as for $(3)$, i.e. $(4)$. This means that ${u_1}_{k_1}$ and $u_k$ span the same set of solutions.
Now, combining $(2)$, $(5)$ and $(7)$, we get
$$\begin{align}
f(x)=u_k(x)\exp(ikx)&=\big[{u_1}_{k_1}(x)\exp(-iGx)\big]\exp(i(k_1+G)x)=\\
                    &={u_1}_{k_1}(x)\exp(ik_1x),
\end{align}\tag9
$$
which is here expressed both in terms of $k$ and in terms of $k_1$, is one and the same solution—for wavenumbers that differ by a whole number of reciprocal lattice constants.
A: I have ultimately found a simple derivation of $\psi_{\vec k}=\psi_{\vec k+\vec G}$ in the text by Ibach and Luth. 
Expanding the periodic potential $U(\vec r)$ and the solution of time-independent Schrdinger equation (TISE) $\psi(\vec r)$ in Fourier series $$U(\vec r)=\sum\limits_{\vec G}U_{\vec G}e^{i\vec G\cdot\vec r},~\psi(\vec r)= \sum\limits_{\vec q}c_{\vec q}e^{i\vec q\cdot\vec r}$$ and plugging into the TISE $$\Big(-\frac{\hbar^2}{2m}\nabla^2+U(\vec r)\Big)\psi=E\psi,$$ it can be shown that the solutions $\psi(\vec r)$ are of the form [Aschroft and Mermin, Solid State Physics, Page 137] $$\psi_{\vec k}(\vec r)=\sum\limits_{\vec G}c_{\vec k-\vec G'}e^{i(\vec k-\vec G')\cdot\vec r}.$$ Now one can follow the a litte piece of algebra as in in Ibach and Luth:
Changing $\vec k\to \vec k+\vec G$, $$\psi_{\vec k+\vec G}=\sum\limits_{\vec G'}c_{k+\vec G-\vec G'}e^{i(\vec k+\vec G-\vec G')\cdot r}.$$ Renaming $\vec G'-\vec G=G^{''}$, $$\psi_{\vec k+\vec G}=\sum\limits_{\vec G^{''}}c_{\vec k+\vec G^{''}}e^{i(\vec k-\vec G^{''})\cdot\vec r}=\psi_{\vec k}(\vec r).$$ Since, $\psi_{\vec k}$ is shown to be periodic, now it is easy to show that the energy is also periodic $E_{\vec k+\vec G}=E_{\vec k}$ which allows us to remain restricted to first Brillouin zone.
A: We have a periodic lattice. Bloch’s theorem ensures that all we have to worry about is the eigenstates at one site that is multiplied by the translator. In other words the overall wavefunction is essentially a periodic function with a period of $2\pi/a$ where $a$ is the lattice constant. 

Looking at spatial variations in the system, the smallest resolution for our translator is the lattice constant $a$. Because the variations within the lattice is captured by our on-site wavefunction $u(r)$. This means that the frequencies we need to look at lie between $0$ and $2\pi/a$. Larger frequencies don’t give us any additional information. 
