Non-diagonal matrix element in the Wannier functions basis

Question

I'm going through Marzari's paper about maximally localized wannier functions. There is a passage I'm trying to understand but can't seem to get. On the third page, between equation 5 and 6, they have the following:

$$ \left\langle u_{nk}|u_{mk+b}\right\rangle = \left\langle\psi_{nk}|e^{-ib\cdot r}|\psi_{mk+b}\right\rangle $$ which is trivial to prove using Bloch's theorem: $\left.|\psi_{nk}\right\rangle = e^{ik\cdot r}\left.|u_{nk}\right\rangle$. However, what I don't get is the following part where they state that the above equation is equal to:

$$ \left\langle\psi_{nk}|e^{-ib\cdot r}|\psi_{mk+b}\right\rangle = \sum_R e^{-ik\cdot R}\left\langle Rn|e^{-ib\cdot r}|0m\right\rangle $$

where $\left.|Rn\right\rangle=\int dk e^{-ik\cdot R}\left.|\psi_{nk}\right\rangle$ is a Wannier function centered in cell R (the inverse transform is $\left.|\psi_{nk}\right\rangle = \sum_R e^{ik\cdot R}\left.|Rn\right\rangle$). I'm trying to prove this to myself but can't figure out how they just "note" this. As per the definition of Wannier functions, I would believe there would be a double sum on Wannier centers like so:

$$ \left\langle\psi_{nk}|e^{-ib\cdot r}|\psi_{mk+b}\right\rangle = \sum_{RR'} e^{-ik\cdot R}\left\langle Rn|e^{-ib\cdot r}|R'm\right\rangle e^{i(k+b)\cdot R'} $$

and I believe it would be possible to use translation/symmetry arguments to rearrange this equation such that the rightmost ket is expressed from the wannier function at the origin $\left.|0m\right\rangle$. However this would not kill one of the sum as the statement suggests.

So my question is: How did they arrive at this equation? My guess is that one could use the properties of the Bloch function to deduce this but I'm not sure how.

As a follow up: this result is somewhat applied to any operator matrix element (like the Hamiltonian). Therefore I believe the proof of this statement should hold for any operator.

N.B.: In all the above, R, r, k and b are all 3D vectors.

Answer 1

Ok I believe I got the answer right after posting the question. The trick is to revert back to real space (using the identity $1=\int dr\left.|r\right\rangle\left\langle r|\right.$:

$$ \sum_{RR'}e^{-ik\cdot R}e^{i(k+b)\cdot R'}\int dr\left\langle Rn|r\right\rangle e^{-ib\cdot r}\left\langle r|R'm\right\rangle $$

but wannier functions are related by translation operation $\left\langle r|R'm\right\rangle = \left\langle r - R'|0m\right\rangle$. Then, doing a change of variable $r\rightarrow r+R'$, you get $$ \sum_{RR'}e^{-ik\cdot R}e^{i(k+b)\cdot R'}\int dr\left\langle Rn|r+R'\right\rangle e^{-ib\cdot (r+R')}\left\langle r|0m\right\rangle = \sum_{RR'}e^{-ik\cdot (R - R')}\int dr\left\langle R-R'n|r\right\rangle e^{-ib\cdot r}\left\langle r|0m\right\rangle$$

where an exponential got cancelled out. From then, we can get rid of the temporary real space integral: $$ = \sum_{RR'}e^{-ik\cdot (R - R')}\left\langle R-R'n|e^{-ib\cdot r}|0m\right\rangle $$

and now we can use a simple symmetry argument. Since we sum on all possible values of $R-R'$ and since that $R-R'$ is itself a wannier center as well (R and R' are part of a Bravais lattice), this double sum is equivalent to a single sum on all bravais lattice vectors. Therefore, we can get rid of the sum on $R'$ for instance and replace any occurance of $R-R'$ for $R$ only. Therefore,

$$ = N\sum_{R}e^{-ik\cdot R}\left\langle Rn|e^{-ib\cdot r}|0m\right\rangle. $$ where $N$ is the number of Wannier functions. Also, in the OP I did not include a factor of $1/\sqrt{N}$ in the definition of the wannier functions but usually it is included and these would cancel out the factor of $N$ appearing here.

As for a more general proof related to any operator, as long as it is diagonal in the real space and have a translational symmetry (which is true for the Hamiltonian) and that we compare matrix elements of the form $\left\langle \psi_{nk}|O(r)|\psi_{mk}\right\rangle$ (at the same k instead of k+b like above), the same argument holds.