Expected value of the current density operator In Ullrich's TD-DFT book, the paramagnetic current-density operator is defined as
$$\hat{\mathbf{j}}(\mathbf{r})=\frac{1}{2i}\sum_{a=1}^{N}\left[\nabla_a\delta(\mathbf{r}-\mathbf{r}_a)+\delta(\mathbf{r}-\mathbf{r}_a)\nabla_a\right]$$
The expectation value of this operator is $\langle\hat{\mathbf{j}}(\mathbf{r})\rangle=\mathbf{j}(\mathbf{r},t)$. I know that the current density is 
$$\mathbf{j}(\mathbf{r},t)=\frac{1}{2i}\left(\Psi^*\nabla\Psi-\Psi\nabla\Psi^*\right),$$
and this is the result I think I have to try to get with $\langle\hat{\mathbf{j}}(\mathbf{r})\rangle$.
At first sight, I have to find some property of the delta function that makes negative one of the terms in the integral of $\hat{\mathbf{j}}(\mathbf{r})$. So, I have
$$
\langle\hat{\mathbf{j}}(\mathbf{r})\rangle=\frac{1}{2i}\sum_{a=1}^{N}\int\Psi^*(\mathbf{x}_1,\dots,\mathbf{x}_N,t)\left[\nabla_a\delta(\mathbf{r}-\mathbf{r}_a)+\delta(\mathbf{r}-\mathbf{r}_a)\nabla_a\right]\Psi(\mathbf{x}_1,\dots,\mathbf{x}_N,t)d\mathbf{x}_1\dots\mathbf{x}_N,$$
where $\mathbf{x}_i$ represents spacial and spin coordinates. I thought that this property could be useful:
$$\int f(x)\delta'(x-a)dx=-f'(a),$$
Thus, for simplicity $N=1$,
$$\begin{align}
\langle\hat{\mathbf{j}}(\mathbf{r})\rangle&=\frac{1}{2i}\left(\int\Psi^*(\mathbf{x}_1,t)\nabla\delta(\mathbf{r}-\mathbf{r}_a)\Psi(\mathbf{x}_1,t)d\mathbf{x}_1+\int\Psi^*(\mathbf{x}_1,t)\delta(\mathbf{r}-\mathbf{r}_a)\nabla\Psi(\mathbf{x}_1,t)d\mathbf{x}_1\right)\\
&=\frac{1}{2i}\left(-\Psi^*(\mathbf{x},t)\nabla\Psi(\mathbf{x},t)+\Psi^*(\mathbf{x},t)\nabla\Psi(\mathbf{x},t)\right)
\end{align}
$$
and this is equal to zero. Maybe I used the property in a questionable way or there is something I am missing. 
Also, my second issue is about obtaining a term with the gradient of $\Psi^*(\mathbf{x},t)$ and I have no clue of what property use (hermiticity?). Can you give me some hints to solve this? Further, this is useful for obtaining the continuity equation, so I need to resolve this problem.
 A: This answer is posted as a request from some users.
The expected value of the paramagnetic current density of an $N$ electron system is
$$
\begin{align}
\langle\hat{\mathbf{j}}(\mathbf{r})\rangle&=\frac{1}{2i}\sum_{a=1}^{N}\langle\Psi|\nabla_a\delta(\mathbf{r}-\mathbf{r}_a)+\delta(\mathbf{r}-\mathbf{r}_a)\nabla_a|\Psi\rangle \\
&=\frac{1}{2i}\sum_{a=1}^{N}\left[\int\Psi^{*}(\{\mathbf{x}_a\},t)\nabla_a\delta(\mathbf{r}-\mathbf{r}_a)\Psi(\{\mathbf{x}_a\},t)d\mathbf{x}_1\dots d\mathbf{x}_N\right.\\
&\qquad\qquad\quad\left.+\int\Psi^{*}(\{\mathbf{x}_a\},t)\delta(\mathbf{r}-\mathbf{r}_a)\nabla_a\Psi(\{\mathbf{x}_a\},t)d\mathbf{x}_1\dots d\mathbf{x}_N\right]\\
&=\frac{1}{2i}\sum_{a=1}^{N}\left[\int\Psi^{*}(\{\mathbf{x}_a\},t)\Psi(\{\mathbf{x}_a\},t)\nabla_a\delta(\mathbf{r}-\mathbf{r}_a)d\mathbf{x}_1\dots d\mathbf{x}_N \right.\\
&\qquad\qquad\quad\left.+2\int\Psi^{*}(\{\mathbf{x}_a\},t)\delta(\mathbf{r}-\mathbf{r}_a)\nabla_a\Psi(\{\mathbf{x}_a\},t)d\mathbf{x}_1\dots d\mathbf{x}_N \right]
\end{align}
$$
Using a property of the derivative of the Dirac delta (found in Zettili's Quantum Mechanics) on the first integral of the square brackets,
$$\int_{-\infty}^\infty f(x)\delta'(x-a)dx=-f'(a),$$
we get (using a $b\neq a$ that runs from $1$ to $N$)
$$
\begin{align}
\langle\hat{\mathbf{j}}(\mathbf{r})\rangle&=\frac{1}{2i}\sum_{a=1}^{N}\left[\int(-1)\Psi^{*}(\mathbf{r},\sigma_a,\{\mathbf{x}_{b\neq a}\},t)\nabla_\mathbf{r}\Psi(\mathbf{r},\sigma_a,\{\mathbf{x}_{b\neq a}\},t)d\sigma_a d\mathbf{x}_b\right.\\
&\qquad\qquad\quad\left.-\int\Psi(\mathbf{r},\sigma_a,\{\mathbf{x}_{b\neq a}\},t)\nabla_\mathbf{r}\Psi^{*}(\mathbf{r},\sigma_a,\{\mathbf{x}_{b\neq a}\},t)d\sigma_ad\mathbf{x}_b\right.\\
&\qquad\qquad\quad+2\left.\int\Psi^{*}(\mathbf{r},\sigma_a,\{\mathbf{x}_{b\neq a}\},t)\nabla_\mathbf{r}\Psi(\mathbf{r},\sigma_a,\{\mathbf{x}_{b\neq a}\},t)d\sigma_ad\mathbf{x}_b\right]\\
&=\frac{1}{2i}\sum_{a=1}^{N}\left[\int\Psi^{*}(\mathbf{r},\sigma_a,\{\mathbf{x}_{b\neq a}\},t)\nabla_\mathbf{r}\Psi(\mathbf{r},\sigma_a,\{\mathbf{x}_{b\neq a}\},t)d\sigma_ad\mathbf{x}_b\right.\\
&\qquad\qquad\quad\left.-\int\Psi(\mathbf{r},\sigma_a,\{\mathbf{x}_{b\neq a}\},t)\nabla_\mathbf{r}\Psi^{*}(\mathbf{r},\sigma_a,\{\mathbf{x}_{b\neq a}\},t)d\sigma_ad\mathbf{x}_b\right]\\
&=\frac{N}{2i}\left[\int\Psi^{*}(\mathbf{r},\sigma_a,\{\mathbf{x}_{b\neq a}\},t)\nabla_\mathbf{r}\Psi(\mathbf{r},\sigma_a,\{\mathbf{x}_{b\neq a}\},t)d\sigma_ad\mathbf{x}_b\right.\\
&\qquad\quad\left.-\int\Psi(\mathbf{r},\sigma_a,\{\mathbf{x}_{b\neq a}\},t)\nabla_\mathbf{r}\Psi^{*}(\mathbf{r},\sigma_a,\{\mathbf{x}_{b\neq a}\},t)d\sigma_ad\mathbf{x}_b\right]\\
&=N\int\text{Im}\left[\Psi^{*}(\mathbf{r},\sigma_a,\{\mathbf{x}_{b\neq a}\},t)\nabla_\mathbf{r}\Psi(\mathbf{r},\sigma_a,\{\mathbf{x}_{b\neq a}\},t)\right]d\sigma_1d\mathbf{x}_b
\end{align}
$$
Obtaining the required value (see this, page 5 and 6):
$$\langle\hat{\mathbf{j}}(\mathbf{r})\rangle=\mathbf{j}(\mathbf{r},t)$$
