Bloch functions vs. Bloch state vectors Let
$$\hat{H}=\frac{\hat{\mathbf{p}}^2}{2m} + V_L(\mathbf{r})$$
where the lattice potential $V_L(\mathbf{r})=V_L(\mathbf{r}+\mathbf{R})$ for any lattice vector $\mathbf{R}$, and let $\hat{T}_{\mathbf{R}}$ denote translation by $\mathbf{R}$. Skipping the details for now, we can derive the following eigenvalue equations:
\begin{align}
    \hat{T}_{\mathbf{R}}|n,\mathbf{k}\rangle&=\exp\left(-i\mathbf{k}\cdot\mathbf{R}\right)|n,\mathbf{k}\rangle\\
    \hat{H}|n,\mathbf{k}\rangle&=E_n\left(\mathbf{k}\right)|n,\mathbf{k}\rangle
\end{align}
The first equation implies that the wave function $\psi_{n\mathbf{k}}\left(\mathbf{r}\right)$ associated with $|n,\mathbf{k}\rangle$ obeys
$$\psi_{n\mathbf{k}}\left(\mathbf{r}+\mathbf{R}\right)=\exp\left(i\mathbf{k}\cdot\mathbf{R}\right)\psi_{n\mathbf{k}}\left(\mathbf{r}\right)$$
and letting $$u_{n\mathbf{k}}\left(\mathbf{r}\right)=\exp\left(-i\mathbf{k}\cdot\mathbf{r}\right)\psi_{n\mathbf{k}}\left(\mathbf{r}\right)$$
we see that $u_{n\mathbf{k}}\left(\mathbf{r}+\mathbf{R}\right)=u_{n\mathbf{k}}\left(\mathbf{r}\right)$, which is Bloch's theorem.
Here's my question. What is $|u_{n\mathbf{k}}\rangle$? It is tempting to claim that $|u_{n\mathbf{k}}\rangle=\exp\left(-i\mathbf{k}\cdot\mathbf{r}\right)|n,\mathbf{k}\rangle$, but this does not seem to be the case:
\begin{align*}
|u_{n\mathbf{k}}\rangle&=\int d^3r|\mathbf{r}\rangle u_{n\mathbf{k}}\left(\mathbf{r}\right)\\
&= \int d^3r|\mathbf{r}\rangle\exp\left(-i\mathbf{k}\cdot\mathbf{r}\right)\psi_{n\mathbf{k}}\left(\mathbf{r}\right)
\end{align*}
and the $\exp\left(-i\mathbf{k}\cdot\mathbf{r}\right)$ factor cannot just be pulled out of the integral.
 A: The trick is to replace the $\mathbb C$-valued function $\exp(i\mathbf k\cdot\mathbf r)$ by the operator $\exp(i\mathbf k\cdot\hat{\mathbf r})$ (which no longers depend on $\mathbf r$.
More explicitely, you are defining $u$ by :
$$ u_{n\mathbf{k}}\left(\mathbf{r}\right)=\exp\left(-i\mathbf{k}\cdot\mathbf{r}\right)\psi_{n\mathbf{k}}\left(\mathbf{r}\right)$$
ie :
$$\langle \mathbf r|u_{n\mathbf k}\rangle  =\exp(-i\mathbf k\cdot\mathbf r)\langle \mathbf r|\psi_{n\mathbf k}\rangle$$
Then, you can compute :
\begin{align}
|u_{n\mathbf k}\rangle  &= \int \text d\mathbf r|\mathbf r\rangle \langle\mathbf r|u_{n\mathbf k}\rangle  \\
&= \int \text d\mathbf r|\mathbf r\rangle\exp(-i\mathbf k\cdot\mathbf r)\langle \mathbf r|\psi_{n\mathbf k}\rangle \\
&= \int \text d\mathbf r\exp(-i\mathbf k\cdot\hat{\mathbf r})|\mathbf r\rangle\langle \mathbf r|\psi_{n\mathbf k}\rangle \\
&= \exp(-i\mathbf k\cdot\hat{\mathbf r})\int \text d\mathbf r|\mathbf r\rangle\langle \mathbf r|\psi_{n\mathbf k}\rangle \\
&= \exp(-i\mathbf k\cdot\hat{\mathbf r})|\psi_{n\mathbf k}\rangle
\end{align}
