Hamiltonian with periodic potential in second quantization I'm working with the following Hamiltonian 
$$\hat{H}=\int\mathrm{d}\mathbf{x}\sum_{\sigma\in\left\lbrace\uparrow,\downarrow\right\rbrace}\hat{\psi}_\sigma^\dagger(\mathbf{x})\left[-\frac{\hbar^2\nabla^2}{2m}+V(\mathbf{x})\right]\hat{\psi}_\sigma(\mathbf{x})$$
(plus a term concerning an interaction potential, which I can do myself when I understand this part.)
The $\hat\psi^\dagger(\mathbf{x})$ and $\hat\psi(\mathbf{x})$ are the fermionic field creation and annihilation operators. $$V(\mathbf{x})=V(\mathbf{x}+a\mathbf{\hat{e}}_i)$$ is a periodic potential with cubic symmetry, where $a$ is the lattice constant. The minima of the potential are denoted by $\mathbf{x}_\mathbf{j}$, with $$\mathbf{j}=j_x \mathbf{\hat{e}}_x+j_y \mathbf{\hat{e}}_y+j_z \mathbf{\hat{e}}_z$$ the vector that labels the minima. The coordinates are chosen in such a way that $\mathbf{x}=0$ coincides with a minimum.
The single particle Schrödinger eqn. is solved according to Bloch's theorem by eigenfunctions $\phi_{\mathbf{k},m}(\mathbf{x})$ with energies $\epsilon_{\mathbf{k},m}$. The Wannier functions $W_m(\mathbf{x}-\mathbf{x_j})$ that describe a localized particle at site $\mathbf{x_j}$ are determined by
$$\phi_{\mathbf{k},m}(\mathbf{x})=\frac{1}{\sqrt{N}}\sum_\mathbf{j}\mathrm{e}^{i\mathbf{k}\cdot\mathbf{x_j}}W_m(\mathbf{x}-\mathbf{x_j})$$
The aim is to write the Hamiltonian in terms of Wannier functions and the creation/annihilation operators, $\hat{c}_{\mathbf{j},m,\sigma}^\dagger$ and $\hat{c}_{\mathbf{j},m,\sigma}$, that create/annihilate an electron with spin state $\left|{\sigma}\right\rangle$ at lattice site $\mathbf{x_j}$ in the band with index $m$.
To do so, I started by expanding the field operators in terms of the Bloch functions and the operators $\hat{c}_{\mathbf{j},m,\sigma}^\dagger$ and $\hat{c}_{\mathbf{j},m,\sigma}$, i.e. I wrote
$$\hat{\psi}(\mathbf{x})=\sum_{\mathbf{j},\mathbf{k},m}\hat{c}_{\mathbf{j},m,\sigma}\phi_{\mathbf{k},m}(\mathbf{x})$$
Plugging this into the Hamiltonian gives me
$$\hat{H}=\sum_{\mathbf{j},\mathbf{j'}}\sum_{\mathbf{k},\mathbf{k'}}\sum_{m,m'}\sum_\sigma\hat{c}_{\mathbf{j},m,\sigma}^\dagger\hat{c}_{\mathbf{j'},m',\sigma}\int\mathrm{d}\mathbf{x}\;\phi_{\mathbf{k},m}^*(\mathbf{x})\left[-\frac{\hbar^2\nabla^2}{2m}+V(\mathbf{x})\right]\phi_{\mathbf{k'},m'}(\mathbf{x})$$
Next, I write this in terms of the Wannier functions:
$$\hat{H}=\sum_{\mathbf{j},\mathbf{j'}}\sum_{\mathbf{k},\mathbf{k'}}\sum_{m,m'}\sum_\sigma\hat{c}_{\mathbf{j},m,\sigma}^\dagger\hat{c}_{\mathbf{j'},m',\sigma}\int\mathrm{d}\mathbf{x}\;W_{m}^*(\mathbf{x}-\mathbf{x_j})\mathrm{e}^{-i\mathbf{k}\cdot\mathbf{x_j}}\left[-\frac{\hbar^2\nabla^2}{2m}+V(\mathbf{x})\right]W_{m}(\mathbf{x}-\mathbf{x_{j'}})\mathrm{e}^{-i\mathbf{k'}\cdot\mathbf{x_{j'}}}$$
However, this is the point where I'm stuck. The final result should be something like
$$\hat{H}=\sum_{\mathbf{j},\mathbf{j'}}\sum_{m,m'}\sum_\sigma t_{\mathbf{j},\mathbf{j'};m,m'}\hat{c}_{\mathbf{j},m,\sigma}^\dagger\hat{c}_{\mathbf{j'},m',\sigma}$$
where
$$t_{\mathbf{j},\mathbf{j'};m,m'}=\int\mathrm{d}\mathbf{x}\;W_{m}^*(\mathbf{x}-\mathbf{x_j})\left[-\frac{\hbar^2\nabla^2}{2m}+V(\mathbf{x})\right]W_{m}(\mathbf{x}-\mathbf{x_{j'}})$$
I need to get rid of the two exponentials somehow, but I can't get it. I think I might have done something wrong in the expansion of $\hat{\psi}_\sigma$, but I don't know what.
Thanks in advance.
 A: I would procede in a different way. I will consider the 1D case for simplicity.
At first you write the field as a superposition of Bloch modes:
\begin{equation*}
\hat{\psi}(\vec{r})= \sum_{s=\{\uparrow,\downarrow\}} \sum_n \sum_{k=1}^M b_{s\,n\,k} \,  e^{i\tilde{k}x} f_{nk}(x)
\end{equation*}
where $n$ labels the bands, and $k$ runs over the first Brillouin zone. $\tilde{k}=\frac{2\pi}{L}k$ is the quantized wave number and $f_{nk}$ is a function with the same periodicity of the lattice. 
Performing the substitution the Hamiltonian reads: 
\begin{equation*}
\hat{H}=\sum_{s=\{\uparrow,\downarrow\}}\sum_{n} \sum_{k=1}^M \epsilon_{nk}\left(b_{s\,n\,k}^\dagger b_{s\,n\,k} \right)
\end{equation*}
At this point you can jump from the momentum modes picture to the spatial modes picture, just by substituing:
\begin{equation*}
b_{s\,n\,k}= \sum_{j=1}^M \frac{1}{\sqrt{M}} a_{s\,n\,j} e^{-i\tilde{k}ja}
\end{equation*}
where $j$ labels the $M$ sites. 
You then need to approximate the dispersion reation with 
\begin{equation*}
\epsilon_n (k) \approx  2\epsilon_n[1-\cos(\tilde{k})] 
\end{equation*}
do the substitution, write $\cos$ in exponential form, recognize the definition of Kronecker delta. At the end you will end up with:
\begin{equation*}
\hat{H}=2\sum_{s=\{\uparrow,\downarrow\}}\sum_n \sum_{j=1}^M \epsilon_n a_{s\,n\,j}^\dagger a_{s\,n\,j} -\sum_{s=\{\uparrow,\downarrow\}}\sum_n \sum_{j=1}^M \left[ a_{s\,n\,j+1}^\dagger a_{s\,n\,j} + a_{s\,n\,j}^\dagger a_{s\,n\,j-1}\right]
\end{equation*}
Let me know if you need the full computation.
A: I eventually solved it using bracket notation.
Define $\left|m,\mathbf{k}\right\rangle$ as the state of an electron with momentum $\mathbf{k}$ in band $m$, i.e. in real space $\phi_{m,\mathbf{k}}(\mathbf{x})=\left\langle\mathbf{x}\,\right|\left.m,\mathbf{k}\right\rangle$. We can then write the Wannier state as
$$
\left|m,\mathbf{j}\right\rangle=\frac{1}{\sqrt{N}}\sum_{\mathbf{k}}\mathrm{e}^{-i\mathbf{k}\cdot\mathbf{x_j}}\left|m,\mathbf{k}\right\rangle
$$
such that $W_m(\mathbf{x}-\mathbf{x_j})=\left\langle \mathbf{x}\,\right|\left.m,\mathbf{j}\right\rangle$. Working in this basis, we have
$$
\begin{align}
\hat{\psi}_\sigma^\dagger(\mathbf{x})\left|0\right\rangle\equiv\left|\mathbf{x},\sigma\right\rangle &=\sum_m\sum_\mathbf{j}\sum_{\sigma'}\left|m,\mathbf{j},\sigma'\right\rangle\left\langle m,\mathbf{j},\sigma'\right|\left.\mathbf{x},\sigma\right\rangle\\ &=\sum_m\sum_\mathbf{j}\left|m,\mathbf{j},\sigma\right\rangle\left\langle m,\mathbf{j}\right|\left.\mathbf{x}\right\rangle\\ &=\sum_m\sum_\mathbf{j}W_m^*(\mathbf{x}-\mathbf{x_j})\,\hat{c}_{m,\mathbf{j},\sigma}^\dagger\left|0\right\rangle
\end{align}
$$
If we plug this into the original Hamiltonian, we will finally get
$$
\begin{align}
\hat{H}&=\sum_{m,m'}\sum_{\mathbf{j},\mathbf{j'}}\sum_\sigma\hat{c}_{m,\mathbf{j},\sigma}^\dagger\hat{c}_{m',\mathbf{j'},\sigma}\int\mathrm{d}\mathbf{x}\,W_m^*(\mathbf{x}-\mathbf{x_j})\left[-\frac{\hbar^2\nabla^2}{2m}+V(\mathbf{x})\right]W_{m'}(\mathbf{x}-\mathbf{x_{j'}})\\ &=\sum_{m,m'}\sum_{\mathbf{j},\mathbf{j'}}\sum_\sigma t_{m,m';\mathbf{j},\mathbf{j'}}\,\hat{c}_{m,\mathbf{j},\sigma}^\dagger\hat{c}_{m',\mathbf{j'},\sigma}
\end{align}
$$
