Electronic (probability) current density is calculated by continuity equation: $$ \nabla\cdot \vec J=-\partial_t \rho \tag{1} $$ here $\vec J$ is current density operator and $\rho=e \psi_r^+\psi_r$ is density. One way to find $\vec J$ is to explicitly calculate time derivative of $\rho$ and use Schrodinger equation. This method is explained here. This method is easy to understand for me.
I am trying to calculate $\vec J$ using Heisenberg equation of motion that is $$ \nabla\cdot \vec J = -\frac{i}{\hbar}[H_0,e\psi_{r'}^+\psi_{r'}]\tag{2} $$ I am stuck at calculation of this commutation relation. I take $H_0$ as Hamiltonian of free particle: $$ H_0 = \frac{\hbar^2}{2m} \int d\vec r \nabla \psi_{r}^+ \nabla \psi_{r} \tag{3} $$ So: $$ [H_0,e\psi_{r'}^+\psi_{r'}]= [H_0,\psi_{r'}^+]\psi_{r'}+\psi_{r'}^+[H_0,\psi_{r'}]\tag{4} $$ And $$ [H_0,\psi_{r'}]=\frac{\hbar^2}{2m} \int d\vec r [\nabla \psi_{r}^+ \nabla \psi_{r},\psi_{r'}]= \frac{\hbar^2}{2m} \int d\vec r \nabla \psi_{r}^+[ \nabla \psi_{r},\psi_{r'}] +\frac{\hbar^2}{2m} \int d\vec r [\nabla \psi_{r}^+ ,\psi_{r'}]\nabla \psi_{r} \tag{5} $$ how to proceed? how to calculate $[\nabla \psi_{r}^+ ,\psi_{r'}]$?
Edit 1:
After a comment, i found that $[\nabla \psi_{r}^+ ,\psi_{r'}]=\nabla[ \psi_{r}^+ ,\psi_{r'}]$. After using commutation rules $[ \psi_{r} ,\psi_{r'}^+]=\delta(r-r')$, I get: $$ \nabla\cdot \vec J = \frac{-i\hbar e}{2m}\int d\vec r \{\psi_r^+ (\nabla_r[\psi_r^+,\psi_{r'}])\nabla \psi_r +\nabla \psi_r^+ (\nabla[\psi_{r},\psi_{r'}^+])\psi_r\}\\ = \frac{-i\hbar e}{2m}\int d\vec r \{\psi_r^+ \nabla_r\delta(r-r')\nabla \psi_r +\nabla \psi_r^+ \nabla\delta(r-r')\psi_r\}\\ = \frac{-i\hbar e}{2m} \{\psi_r^+ \nabla^2 \psi_r +\nabla \psi_r^+ (\nabla)\psi_r\} $$ The very term $\nabla\psi_r^+ (\nabla)\psi_r$, it should be $(\nabla^2\psi_r^+)\psi_r$. I wonder how the extra $\nabla$ moves back?
