Dipole matrix elements for Bloch wavefunctions? Suppose we have a one-dimensional periodic system with lattice constant $a_0$. From Bloch's theorem, we can express the wavefunction for an electron in band $m$ with crystal momentum $k$ $\left\langle x \middle| \psi_{m,k} \right\rangle$ as follows:
$$
\left\langle x \middle| \psi_{m,k} \right\rangle = e^{i k x} u_{m,k}(x),
$$
where $u_{m,k}(x + a_0) = u_{m,k}(x)$. I don't understand the following expression for the matrix elements of the position operator:
$$
\langle \psi_{m,k} | x | \psi_{m',k'}\rangle = i \delta_{m,m'} \delta_{k,k'} \frac{\partial}{\partial k} + i \delta_{k,k'} X_{m,m'},
$$
where
$$
X_{m,m'}= i N \int_0^{a_0} e^{i (k - k') x} u^*_{m,k}(x) \frac{\partial}{\partial k} u_{m',k'}(x) dx.
$$
The second term is easy enough to understand. The first term, however... how can the expectation value between two eigenstates be a derivative? Am I missing something obvious here?
 A: I think that the correct expression involves $\delta'(k-k')$ rather than a partial derivative.
Suppose  we construct a wave-packet with states taken from a single band with index  $n$:
$$
\varphi({\bf r}) = \int_{Z^*} \frac{d^3k} {(2\pi)^3}a({\bf k}) \psi_{n,{\bf k}}({\bf r})=
\int_{Z^*} \frac{d^3k} {(2\pi)^3}a({\bf k}) e^{i{\bf k}\cdot {\bf r}} u_{n,{\bf k}}({\bf r}).
$$
Here $Z^*$ is the Brillouin zone.
The state  is normalized if
$$
\int_{Z^*}   \frac{d^3k} {(2\pi)^3} | a({\bf k})|^2 =1.
$$
The   spatially localized wave-packet  is now in the domain of the unbounded operator  $\bf r$ and
$$
\langle{\psi_{n',{\bf k}}}|{r_\mu}|{\varphi}\rangle= i\left\{\frac {\partial}{\partial k_\mu} a({\bf k})\delta_{n'n}  +\langle{u_{n'{\bf k}}}|{\partial_\mu  u_{n,{\bf k}}}\rangle  a({\bf k})\right\}.
$$
A: 
The second term is easy enough to understand. The first term, however... how can the expectation value between two eigenstates be a derivative? Am I missing something obvious here?

There is nothing special going on here $-$ this is a pretty common feature.
For contrast, consider the matrix element of the position operator between two (free-particle) momentum eigenstates:
$$
⟨p|\hat x|p'⟩ = i\hbar \frac{\partial}{\partial p} \delta(p-p'),
$$
where you do need to take due care about where the derivative is acting on, but which is basically identical to your expression.
(Note that it is often easier to encapsulate all the complexity by "calculating" the derivative as $\delta'(p-p')$, as mike stone's answer does, but if you do this then it is essential to keep in mind that the derivative of the Dirac delta is a distribution and needs to be treated as such, i.e., not as a function, but as a functional.)
Something equivalent happens for the matrix element of momentum between two position eigenstates:
$$
⟨x|\hat p|x'⟩ = -i\hbar \frac{\partial}{\partial x} \delta(x-x').
$$
In other words, the behaviour you're seeing there is generic.
